The application of random matrix techniques in QCD and non-Abelian gauge theories in general has a long history e.g. in counting Feynman diagrams, going back to 't Hooft and others. In this talk I will focus on a different aspect that relates the two in the low energy spectrum of the QCD Dirac operator, as initiated by Shuryak and Verbaarschot. First, I will explain what is the approximation studied here where spectral statistics of random matrices apply, and where for example the technique of orthogonal polynomials can be useful in comparing to QCD lattice data. It is given by a particular finite volume low energy limit, the epsilon regime of chiral perturbation theory of Gasser and Leutwyler. I will mention how QCD parameters like quark masses, zero-modes, finite lattice spacing or chemical potential can be incorporated into the random matrix ensemble. In the last part I will discuss some recent work with my former student Tim Wurfel on the inclusion of finite temperature, that leads out of the standard classes of random matrices, but still remains analytically tractable. This talk is mainly based on the review arXiv:1603.06011 and the paper with Tim arXiv:2110.03617; part of the London TQFT Journal Club; it will be possible to follow this talk online (please register at https://london-tqft.vercel.app)