We develop a method to study quantum impurity models, small interacting quantum systems bilinearly coupled to an environment, in presence of an additional Markovian quantum bath, with a generic non-linear coupling to the impurity. We aim at computing the evolution operator of the reduced density matrix of the impurity, obtained after tracing out all the environmental degrees of freedom. First, we derive an exact real-time hybridization expansion for this quantity, which generalizes the result obtained in absence of the additional Markovian dissipation, and which could be amenable to stochastic sampling through diagrammatic Monte Carlo. Then, we obtain a Dyson equation for this quantity and we evaluate its self-energy with a resummation technique known as the Non-Crossing Approximation. We apply this novel approach to a simple fermionic impurity coupled to a zero temperature fermionic bath and in presence of Markovian pump, losses and dephasing.