Abstract: We study in detail the impact of dynamical quarks on the gluon mass generation mechanism, in the Landau gauge, for the case of a small number of quark families. As in earlier considerations, we assume that the main bulk of the unquenching corrections to the gluon propagator originates from the fully dressed quark-loop diagram. The nonperturbative evaluation of this diagram provides the key relation that expresses the unquenched gluon propagator as a deviation from its quenched counterpart. This relation is subsequently coupled to the integral equation that controls the momentum evolution of the effective gluon mass, which contains a single adjustable parameter; this constitutes a major improvement compared to the analysis presented in Aguilar et al. [Phys. Rev. D 86, 014032 (2012)], where the behavior of the gluon propagator in the deep infrared was estimated through numerical extrapolation. The resulting nonlinear system is then treated numerically, yielding unique solutions for the modified gluon mass and the quenched gluon propagator, which fully confirms the picture put forth recently in several continuum and lattice studies. In particular, an infrared finite gluon propagator emerges, whose saturation point is considerably suppressed, due to a corresponding increase in the value of the gluon mass. This characteristic feature becomes more pronounced as the number of active quark families increases, and can be deduced from the infrared structure of the kernel entering in the gluon mass equation.