I learned the second congruence as a version of Wolsteholme's theorem, and I would be a bit surprised if Kazandzidis was the first person to observe the equivalence between this form and any other form of Wolsteholme's result. As for the "reason" that this result is true, I wrote a proof for the Wikipedia page which is mostly, but not entirely, a direct counting argument and which you could call "the grounds" for the congruence. The conceptual content of the argument is as follows:
The result holds modulo $p^2$ because you can divide the $pm$-set into $m$ cycles of length $p$ and rotate them separately. You obtain $\binom{m}{n}$ equivalence classes of subsets of size 1, and the other equivalence classes have order $p^2$ or higher.
Then you can examine $\binom{-p}{p}$, interpreted as $\binom{p^4-p}{p}$. This binomial coefficient is algebraically equivalent to $\binom{2p}{p}/2$ up to sign. The orbit decomposition in the previous paragraph establishes a second relation with $\binom{2p}{p}$. The conclusion is that the mod $p^2$ contribution vanishes for one non-trivial pair $(m,n)$, and if it vanishes once, it vanishes always. This vanishing principle extends mod $p^4$ — your second binomial congruence holds mod $p^4$ — provided that it vanishes once "by accident". A prime for which this happens is called a Wolstenholme prime. Two such primes are known, 16843 and 2124679.
Another remark: There are two pieces of evidence that the $p^2$ congruence (Babbage's theorem) is entirely combinatorial, but the extension to $p^3$ (Wolstenholme) is essentially algebraic. First, that Wolstenholme's theorem doesn't hold for the primes $p=2,3$. Second, that Babbage's theorem has a $q$-analogue for Gaussian binomial coefficients, with the same orbit proof. But Wolsteholme's extension does not have a $q$-analogue as far as I know. In fact, you can tell that the known $q$-Babbage's theorem doesn't extend, because the difference polynomial doesn't have any more cyclotomic factors.
I stand corrected on a couple of points. First, on the math, there is a paper by George Andrews, "q-Analogs of the binomial coefficient congruences of Babbage, Wolstenholme, and Glaisher", that does give some version of a q-analogue of Wolstenholme's theorem. The proof given there is algebraic, so it does not shed light on possible combinatorial "explanations". Second, according to Andrews, the full binomial interpretation of Wolstenholme's result is due to Glaisher.
Here as an approach.
It is known that:
$P_{n}(x) = \sum_{k=0}^{n}\binom{n}{k}\binom{-n-1}{k}\left(\frac{1-x}{2}\right)^{k}$ (see for example http://en.wikipedia.org/wiki/Legendre_polynomials)
Therefore these expressions are equal
$$
\sum_{k=0}^{n}\binom{n}{k}\binom{n+k}{k}\sum_{j=0}^{k}\binom{k}{j}^{3}=\sum_{k=0}^{n}\binom{n}{k}\binom{-n-1}{k}\frac{1}{4 \pi^{2}}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\left(\frac{1-g(x,y)}{2}\right)^{k}dydx
$$
if
$$
\sum_{j=0}^{k}\binom{k}{j}^{3}=\frac{1}{4 \pi^{2}}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\left(\frac{g(x,y)-1}{2}\right)^{k}dydx
$$
One can see that
$$
\left(\frac{g(x,y)-1}{2}\right)^{k} = 4^{k} (\cos^{2}x + \cos x \sin y )^k
$$
Therefore we get
$$
\frac{1}{4 \pi^{2}}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\left(\frac{g(x,y)-1}{2}\right)^{k}dydx = \frac{4^k}{4\pi^{2}}\int_{0}^{2\pi}\int_{0}^{2\pi}\sum_{j=0}^{k}\binom{k}{j}\cos^{2k-j}x \sin^{j} y dxdy
$$
Since
\begin{align*}
\int_{0}^{2\pi} \cos^{m} x dx=\int_{0}^{2\pi} \sin^{m} x dx = \frac{2\pi} {2^m}\binom{m}{m/2} \quad \text{for} \quad m \quad \text{even, otherwise =0}
\end{align*}
We get
$$
\frac{1}{4 \pi^{2}}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\left(\frac{g(x,y)-1}{2}\right)^{k}dydx = \sum_{\ell=0}^{k/2}\binom{k}{2\ell}\binom{2\ell}{\ell}\binom{2(k-\ell)}{k-\ell}
$$
and now we want to show that
$$
\sum_{j=0}^{k}\binom{k}{j}^{3}=\sum_{\ell=0}^{k/2}\binom{k}{2\ell}\binom{2\ell}{\ell}\binom{2(k-\ell)}{k-\ell}
$$
The last equality follows from the identity:
$$
\sum_{j=0}^{k}\binom{k}{j}^{3}=\sum_{\ell=k/2}^{k}\binom{k}{\ell}^{2}\binom{2\ell}{k}
$$
See for example http://arxiv.org/pdf/math/0311195v1.pdf formula (2).
because
$$
\binom{k}{2\ell}\binom{2\ell}{\ell}\binom{2k-2\ell}{k-\ell}=\binom{k}{k-\ell}^{2}\binom{2k-2\ell}{k}
$$
Best Answer
Here are some computations that may be useful. Here below, $\big[ \,f\,\big]_a^b$ denotes $f(b)-f(a)$, and $n^{(k)}:=n(n-1)\dots(n-k+1)$.
1. (An integration by parts). Let $w$ be a solution of the third order linear ODE on the interval $[a,b]$: $$(x^4-34x^3+x^2)w''' + 3(2x^3-51x^2+x)w''+(7x^2-112x+1)w'+(x-5)w=0,$$ and put $$M(n):=\int_a^b x^nw(x)dx$$ for any integer $n\ge0$. Then, for any $n\ge 2$, $$n^3M(n)- (34n^3-51n^2+27n-5)M(n-1) + (n-1)^3M(n-2)= \Big[A(x)(x^n)'' + B(x)(x^n)' + C(x)x^n \Big]_a^b,$$ where $$ A:=pw \qquad B:= qw- (pw)' \qquad C:= (pw )'' - (qw)'+ rw ,$$ and $$p(x):=x^3-34x^2+x\qquad q(x):=3x^2-51x\qquad r(x):=x+10-x^{-1}.$$
Proof. We express the above polynomials of $n$ in the base $n^{(k)}$, then we absorb the latter terms as coefficients of derivatives of $x^n$, and finally we integrate by parts.
We have: $$n^3M(n)- (34n^3-51n^2+27n-5)M(n-1) + (n-1)^3M(n-2)=$$ $$\int_a^b\Big\{n^3x^nw -(34n^3-51n^2+27n-5)x^{n-1}w +(n-1)^3x^{n-2}w\Big\}dx=$$ $$\int_a^b\Big\{(n^ {(3)}+3n^ {(2)}+ n)x^nw -(34n^{(3)}+51n ^{(2)}+10n-5)x^{n-1}w +(n^{(3)}+n-1)x^{n-2}w\Big\}dx=$$ $$\int_a^b\Big\{(x^n)'''pw+(x^n)''qw+(x^n)'rw+x^n(5x^{-1} -x^{-2} )w\}dx=$$ $$ \int_a^bx^n \Big\{-(pw)'''+(qw)''-(rw)'+ (5x^{-1} -x^{-2} )w\}dx +\Big[ (x^n)''pw- (x^n)'(pw)'+ x^n(pw)'' +(x^n)'qw-x^n(qw)' +rw \Big]_a^b=$$ $$ -\int_a^bx^{n-1} \Big\{(x^4-34x^3+x^2)w''' + 3(2x^3-51x^2+x)w''+(7x^2-112x+1)w'+(x-5)w\Big\}dx+$$ $$+\Big[ (x^n)''pw+ (x^n)'\big(qw- (pw)'\big)+ x^n\big( (pw )'' - (qw)'+ rw \big) \Big]_a^b=$$ $$=\Big[A(x)(x^n)'' + B(x)(x^n)' + C(x)x^n \Big]_a^b.\qquad\square$$
2.(Consequence). Let $w$ a solution of the above linear equation on $(0,c)\setminus\{c_0\}$, with $\int_0^c w(x)dx=1$, and and assume it verifies the following linear boundary conditions, expressed in terms of the above coefficients $A,B,C$:
i) $A(x)=o(1),\quad B(x)=o(1),\quad C(x)=O(1)$, as $x\to0$;
ii) $A(x),\quad B(x),\quad C(x)$, are continuous at $x=c_0$ ,
iii) $A(x)=o(1) ,\quad B(x)=o(1) ,\quad C(x)=o(1)$, as $x\to c$.
Then the corresponding $M(n)$ are the Apéry sequence.
Indeed, computing the integral on $[0,c]$ for $M(n)$ as limit of integrals on $[\epsilon, c_0-\epsilon]\cup[c_0+\epsilon, c-\epsilon]$ as $\epsilon\to0$, and applying the above integration by parts formula, one gets that $M(n)$ satisfy the Apéry's recurrence, with $M(0)=1$ (note that $M(1)=5M(0)$ follows from the recurrence as well).
rmk. This also include the case where $w$ vanishes identically on $(0,c_0)$, and the condition is simply that $A, B, C$ should vanish both at $x=c_0$ and at $x=c$.
3. (Positive solutions of the third order ODE). Assume $u$ solves the second order linear ODE $$(x^3-34x^2+x)u''+(2x^2-51x+1)u'+\frac{1}{4}(x-10)u=0.$$ Then $w:=u^2$ solves $$(x^4-34x^3+x^2)w''' + 3(2x^3-51x^2+x)w''+(7x^2-112x+1)w'+(x-5)w=0.$$
Proof. Put $$P:=x^4-34x^3+x^2\qquad Q:=\frac{x^2}{2}-5x,$$ so the equation for $u$ (multiplied by $2x$) writes: $$2Pu''+P'u'+Qu=0,$$ and we have $$0=(2Pu''+P'u'+Qu)'u+ 3(2Pu''+P'u'+Qu)u'=$$ $$=(2Pu'''+2P'u''+P'u''+P''u'+Qu'+Q'u)u+3(2Pu''+P'u'+Qu)u'=$$ $$=P(2u'''u+6u''u')+3P'(u''u+u'^2)+(P''+4Q)u'u+Q'u^2=$$ $$=Pw''' +\frac{3}{2}P'w''+\Big(\frac{P''}{2}+2Q\Big)w'+Q'w,$$ which is the above third order equation for $w$.