This is originally due to Ljunggren and Jacobsthal, 1952,  in "On the divisibility of the difference between two binomial coefficients" in Skand. Mat.-Kongr., Trondheim 1949, 42-54 (1952).

It is also attributed to G. S. Kazandzidis, who proved a somewhat stronger result. See his papers

 - "On a congruence and on a practical method for finding the highest power of a prime p which divides the binomial coefficient (AB)", Bull. Soc. Math. Grèce, N. Sér. 6, No. 2, 358-360 (1965).
 - "Congruences on the binomial coefficients", Bull. Soc. Math. Grèce, N. Ser. 9, No. 1, 1-12 (1968).
 - "On congruences in number-theory", Bull. Soc. Math. Grèce, N. Sér. 10, No. 1, 35-40 (1969).

Kazandzidis proved $$ \binom{pn}{pk} \equiv \binom{n}{k} \bmod p^{3+a}$$
where $a$ is the $p$-adic valuation of $nk(n-k)\binom{n}{k}$.

All these proofs are elementary. The modern treatment is a bit less so, but is more illuminating. Indeed, $p$-adic analysis can lead to a proof. See Chapter 7 of Alain M. Robert's book "A course in p-adic analysis" (GTM 198, Springer, 2000). The proof given there is based on the material in the paper "The Kazandzidis supercongruences. A simple proof and an application", Rend. Semin. Mat. Univ. Padova 94, 235-243 (1995), by Alain M. Robert and Maxime Zuber.

In the language of $p$-adic analysis, this congruence can be recast in terms of the $p$-adic Gamma function. In the same way that $\Gamma$ is a useful analytic continuation of the factorial faction to the complex plane, $\Gamma_p$, defined in $\mathbb{Z}_p$, is a continuous $p$-adic function that extends the usual Gamma function (well, not quite: to get continuity you must work with $\Gamma_p(n) := (-1)^n\prod_{1 \le m < n, \, p \nmid n}m$ -- the factorial does not extend to a continuous function on the $p$-adics). Then the relevant congruences are a direct consequence of the $p$-adic Taylor expansion of $\log \Gamma_p$. See the aforementioned Chapter 7 for full details.

Some remarks:

 - The congruence $\binom{ap^n}{bp^n} \equiv \binom{ap^{n-1}}{bp^{n-1}} \equiv p^{2n}$ has a combinatorial proof (as opposed to algebraic proofs mod $p^{3n}$). A reference for a combinatorial proof when the modulus is $p^n$ is Gian-Carlo Rota and Bruce Sagan's paper "Congruences derived from group action", Eur. J. Comb. 1, 67-76 (1980). I am not sure if this is the earliest reference of this kind of argument, though.
 - There are various $q$-analogues of this $p$-adic congruences.
 - The above papers include also the cases $p=3$ (and maybe $p=2$?) which are just a bit weaker.