We establish a new connection between moments of {n times n} random matrices Xn and hypergeometric orthogonal polynomials. Specifically, we consider moments {mathbb{E}{rm 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. (J Math Phys 57:111901, 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 rightarrow 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.