Up to now, the nonparametric analysis of multidimensional continuous-time Markov processes has focussed strongly on specific model choices, mostly related to symmetry of the semigroup. While this approach allows to study the performance of estimators for the characteristics of the process in the minimax sense, it restricts the applicability of results to a rather constrained set of stochastic processes and in particular hardly allows incorporating jump structures. As a consequence, for many models of applied and theoretical interest, no statement can be made about the robustness of typical statistical procedures beyond the beautiful, but limited framework available in the literature. To close this gap, we identify $β$-mixing of the process and heat kernel bounds on the transition density as a suitable combination to obtain $\sup$-norm and $L^2$ kernel invariant density estimation rates matching the case of reversible multidimenisonal diffusion processes and outperforming density estimation based on discrete i.i.d. or weakly dependent data. Moreover, we demonstrate how up to $\log$-terms, optimal $\sup$-norm adaptive invariant density estimation can be achieved within our general framework based on tight uniform moment bounds and deviation inequalities for empirical processes associated to additive functionals of Markov processes. The underlying assumptions are verifiable with classical tools from stability theory of continuous time Markov processes and PDE techniques, which opens the door to evaluate statistical performance for a vast amount of Markov models. We highlight this point by showing how multidimensional jump SDEs with Lévy driven jump part under different coefficient assumptions can be seamlessly integrated into our framework, thus establishing novel adaptive $\sup$-norm estimation rates for this class of processes.