  In this article we use the Schwinger-Dyson equations to compute the
nonperturbative modifications caused to the infrared finite gluon propagator
(in the Landau gauge) by the inclusion of a small number of quark families. Our
basic operating assumption is that the main bulk of the effect stems from the
"one-loop dressed" quark loop contributing to the full gluon self-energy. This
quark loop is then calculated, using as basic ingredients the full quark
propagator and quark-gluon vertex; for the quark propagator we use the solution
obtained from the quark gap equation, while for the vertex we employ suitable
Ans\"atze, which guarantee the transversality of the answer. The resulting
effect is included as a correction to the quenched gluon propagator, obtained
in recent lattice simulations. Our main finding is that the unquenched
propagator displays a considerable suppression in the intermediate momentum
region, which becomes more pronounced as we increase the number of active quark
families. The influence of the quarks on the saturation point of the propagator
cannot be reliably computed within the present scheme; the general tendency
appears to be to decrease it, suggesting a corresponding increase in the
effective gluon mass. The renormalization properties of our results, and the
uncertainties induced by the unspecified transverse part of the quark-gluon
vertex, are discussed. Finally, the dressing function of the gluon propagator
is compared with the available unquenched lattice data, showing rather good
agreement.
