The scale evolution of parton distributions is determined by universal splitting functions. As a milestone towards the computation of these functions to four-loop order in QCD, we compute all contributions to the pure-singlet quark-quark splitting functions that involve two closed fermion loops. The splitting functions are extracted from the pole terms of off-shell operator matrix elements, and the workflow for their calculation is outlined. We reproduce known results for the non-singlet four-loop splitting functions and validate our new pure-singlet results against fixed Mellin moments.