We develop the basic theory of eigenvalues of $p$-adic random matrices, analogous to the classical theory for random matrices over $\mathbb{R}$ and $\mathbb{C}$. Such eigenvalue statistics were proposed as a model for the zeroes of $p$-adic $L$-functions by Ellenberg-Jain-Venkatesh, who computed the limiting distribution of the number of eigenvalues in a unit disc. We compute the full joint distribution of the $n$ eigenvalues of an $n \times n$ matrix with Haar distribution, obtaining Coulomb gas type formulas as in the archimedean case, with Vandermonde terms leading to eigenvalue repulsion. From these Coulomb gas density functions we derive asymptotics of eigenvalue statistics as $n \to \infty$. These include exact computations, such as a closed form $$ρ(x,y) = 1 - θ_3(-\sqrt{p};||x-y||^2/p)$$ for the limiting pair correlation of eigenvalues in $\mathbb{Z}_p$, and similar results in quadratic extensions. Such formulas yield concrete numerical predictions on zeroes of $p$-adic $L$-functions. For eigenvalues in arbitrary extensions of $\mathbb{Q}_p$ we also give precise estimates on their pair-repulsion and expected number of eigenvalues in each extension. Finally, we compute the asymptotic probability that all eigenvalues lie in $\mathbb{Z}_p$. Our proofs combine results from several distinct areas: $p$-adic orbital integrals, roots of random $p$-adic polynomials, the Sawin-Wood moment method for random modules, and Markov chains associated with measures on integer partitions.