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 (M2GP) model, which embeds this multi-stage structure within a novel non-stationary covariance function. We show that the M2GP 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 M2GP 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.