We study monic polynomials Qn(x) generated by a high order three-term recursion xQn(x) = Qn+1(x) + an−pQn−p(x) with arbitrary p ≥ 1 and an > 0 for all n. The recursion is encoded by a two-diagonal Hessenberg operator H. One of our main results is that, for periodic coefficients an and under certain conditions, the Qn are multiple orthogonal polynomials with respect to a Nikishin system of orthogonality measures supported on starlike sets in the complex plane. This improves a recent result of Aptekarev-Kalyagin-Saff where a formal connection with Nikishin systems was obtained in the case whenAn important tool in this paper is the study of 'Riemann-Hilbert minors', or equivalently, the 'generalized eigenvalues' of the Hessenberg matrix H. We prove interlacing relations for the generalized eigenvalues by using totally positive matrices. In the case of asymptotically periodic coefficients an, we find weak and ratio asymptotics for the Riemann-Hilbert minors and we obtain a connection with a vector equilibrium problem. We anticipate that in the future, the study of Riemann-Hilbert minors may prove useful for more general classes of multiple orthogonal polynomials.