  We show that when the QCD axion is directly coupled to quarks with $c_q/f \,
\partial_\mu a \, \bar{q} \gamma^\mu \gamma^5 q$, such as in DFSZ models, the
dominant production mechanism in the early universe at temperatures $1 \, {\rm
GeV}\lesssim T \lesssim 100 \,{\rm GeV}$ is obtained via $q \bar{q}
\leftrightarrow g a$ and $q g \leftrightarrow q a$, where $g$ are gluons.
Different heavy quarks $q_i$ can produce a thermal axion background that
decouples at a temperature $T_i$: (1) top quark at $T_t \lesssim 100 \, {\rm
GeV}$ for $f/c_t \lesssim 3\times 10^8 {\rm GeV}$; (2) bottom quark at $T_b
\lesssim m_b$, for $f/c_b\lesssim 8\times 10^{7} {\rm GeV}$; (3) charm quark at
$T_c \lesssim m_c $ for $f/c_c \lesssim 5\times 10^{7} {\rm GeV}$. Each of
these cases corresponds to a contribution to the effective number of
relativistic degrees of freedom, in the windows given by $0.027 \leq \Delta
N_\text{eff}\leq 0.031$, $0.037 \leq \Delta N_\text{eff} \leq 0.039$ and $0.039
\leq \Delta N_\text{eff}$, respectively. These contributions are larger than
the one obtained when thermalization happens only above the electroweak phase
transition, $\Delta N_{\rm eff}\lesssim 0.027$, and are within reach of future
CMB S4 experiments, thus opening an alternative window to detect the axion and
to test the early universe at such temperatures.
