We consider Gaussian and bootstrap approximations for the supremum of additive functionals of aperiodic Harris recurrent Markov chains. The supremum is taken over a function class that may depend on the sample size, which allows for non-Donsker settings; that is, the empirical process need not have a weak limit in the space of bounded functions. We first establish a non-asymptotic Gaussian approximation error, which holds at rates comparable to those for sums of high-dimensional independent or one-dependent vectors. Key to our derivation is the Nummelin splitting technique, which enables us to decompose the chain into either independent or one-dependent random blocks. Additionally, building upon the Nummelin splitting, we propose a Gaussian multiplier bootstrap for practical inference and establish its finite-sample guarantees in the strongly aperiodic case. Finally, we apply our bootstrap to construct a uniform confidence band for an invariant density within a certain class of diffusion processes.