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.