The entanglement entropy is a unique probe to reveal universal features of strongly interacting many-body systems. In two or more dimensions these features are subtle, and detecting them numerically requires extreme precision, a notoriously difficult task. This is especially challenging in models of interacting fermions, where many such universal features have yet to be observed. In this Letter we tackle this challenge by introducing a new method to compute the Rényi entanglement entropy in auxiliary-field quantum Monte Carlo simulations, where we treat the entangling region itself as a stochastic variable. We demonstrate the efficiency of this method by extracting, for the first time, universal subleading logarithmic terms in a two-dimensional model of interacting fermions, focusing on the half-filled honeycomb Hubbard model at T=0. We detect the universal corner contribution due to gapless fermions throughout the Dirac semi-metal phase and at the Gross-Neveu-Yukawa critical point, where the latter shows a pronounced enhancement depending on the type of entangling cut. Finally, we observe the universal Goldstone mode contribution in the antiferromagnetic Mott insulating phase.