In an era where scientific experimentation is often costly, multi-fidelity emulation provides a powerful tool for predictive scientific computing. While there has been notable work on multi-fidelity modeling, existing models do not incorporate an important multi-stage property of multi-fidelity simulators, where multiple fidelity parameters control for accuracy at different experimental stages. Such multi-stage simulators are widely encountered in complex nuclear physics and astrophysics problems. We thus propose a new Multi-stage Multi-fidelity Gaussian Process (M$^2$GP) model, which embeds this multi-stage structure within a novel non-stationary covariance function. We show that the M$^2$GP model can capture prior knowledge on the numerical convergence of multi-stage simulators, which allows for cost-efficient emulation of multi-fidelity systems. We demonstrate the improved predictive performance of the M$^2$GP model over state-of-the-art methods in a suite of numerical experiments and two applications, the first for emulation of cantilever beam deflection and the second for emulating the evolution of the quark-gluon plasma, which was theorized to have filled the Universe shortly after the Big Bang.