Numerical Stochastic Perturbation Theory (NSPT) allows for perturbative computations in quantum field theory. We present an implementation of NSPT that yields results for high orders in the perturbative expansion of lattice gauge theories coupled to fermions. The zero-momentum mode is removed by imposing twisted boundary conditions; in turn, twisted boundary conditions require us to introduce a smell degree of freedom in order to include fermions in the fundamental representation. As a first application, we compute the critical mass of two flavours of Wilson fermions up to order
in a
gauge theory. We also implement, for the first time, staggered fermions in NSPT. The residual chiral symmetry of staggered fermions protects the theory from an additive mass renormalisation. We compute the perturbative expansion of the plaquette with two flavours of massless staggered fermions up to order
in a
gauge theory, and investigate the renormalon behaviour of such series. We are able to subtract the power divergence in the Operator Product Expansion (OPE) for the plaquette and estimate the gluon condensate in massless QCD. Our results confirm that NSPT provides a viable way to probe systematically the asymptotic behaviour of perturbative series in QCD and, eventually, gauge theories with fermions in higher representations.