We establish a new connection between moments of n×n random matrices Xn and hypergeometric orthogonal polynomials. Specifically, we consider moments ETrX−sn as a function of the complex variable s∈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. 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→∞ 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.