  We establish a new connection between moments of $n \times n$ random matrices
$X_n$ and hypergeometric orthogonal polynomials. Specifically, we consider
moments $\mathbb{E}\mathrm{Tr} X_n^{-s}$ as a function of the complex variable
$s \in \mathbb{C}$, whose analytic structure we describe completely. We
discover several remarkable features, including a reflection symmetry (or
functional equation), zeros on a critical line in the complex plane, and
orthogonality relations. An application of the theory resolves part of an
integrality conjecture of Cunden et al. [F. D. Cunden, F. Mezzadri, N. J. Simm
and P. Vivo, J. Math. Phys. 57 (2016)] on the time-delay matrix of chaotic
cavities. In each of the classical ensembles of random matrix theory (Gaussian,
Laguerre, Jacobi) we characterise the moments in terms of the Askey scheme of
hypergeometric orthogonal polynomials. We also calculate the leading order
$n\to\infty$ asymptotics of the moments and discuss their symmetries and
zeroes. We discuss aspects of these phenomena beyond the random matrix setting,
including the Mellin transform of products and Wronskians of pairs of classical
orthogonal polynomials. When the random matrix model has orthogonal or
symplectic symmetry, we obtain a new duality formula relating their moments to
hypergeometric orthogonal polynomials.
