We perform for the first time a dedicated analysis of cosmological constraints on Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) QCD axion models. Such constructions are especially interesting in light of the recent XENON rrexcess and of hints from stellar cooling. In DFSZ models, for m(a) greater than or similar to 0.1 eV, scatterings of pions and muons can produce a sizable cosmic background of thermal axions, that behave similarly to massive neutrinos. However, the pion coupling depends on the alignment between the vacuum expectation value (vevs) of two Higgs doublets, and can be significantly suppressed or enhanced with respect to the Kim-Shifman-Vainshtein-Zakharov scenario (KSVZ). Using the latest Planck and BAO data, we find m(a) <= 0.2 eV at 95% C.L., when the axion coupling to pions c(a pi) is maximal. Constraints on m(a), instead, can be significantly relaxed when c(a pi) is small. In particular, we point out that in the so-called DFSZ-II model, where the axion coupling to leptons does not vanish simultaneously with c(a pi), production via muons gives m(a) < 0.6 eV at 95% C.L., whereas in the DFSZ-I model bounds on m a can be fully lifted. We then combine cosmological data with recent hints of a DFSZ axion coupled to electrons from the XENON1T experiment, finding in this case that the axion mass is constrained to be in the window 0.07 eV less than or similar to m(a) less than or similar to 1.8(0.3) eV for the DFSZ-I (DFSZ-II) model. A similar analysis with stellar cooling hints gives 3 meV less than or similar to m(a) less than or similar to 0.2 eV for DFSZ-II, while no constraint arises in the DFSZ-I case. Forthcoming cosmic microwave background stage 4 experiments will be able to further test such scenarios; for instance the XENON1T window should be fully probed at 2 sigma for a DFSZ-I axion.