Axions, if realized in nature, can be copiously produced in the early universe via thermal processes, contributing to the mass-energy density of thermal hot relics. In light of the most recent cosmological observations, we analyze two different thermal processes within a realistic mixed hot-dark-matter scenario which includes also massive neutrinos. Considering the axion-gluon thermalization channel we derive our most constraining bounds on the hot relic masses $m_a < 7.46$ eV and $sum m_ u< 0.114$ eV both at 95 per cent CL; while studying the axion-pion scattering, without assuming any specific model for the axion-pion interactions and remaining in the range of validity of the chiral perturbation theory, our most constraining bounds are improved to $m_a<0.91$ eV and $sum m_ u< 0.105$ eV, both at 95 per cent CL. Interestingly, in both cases, the total neutrino mass lies very close to the inverted neutrino mass ordering prediction. If future terrestrial double beta decay and/or long baseline neutrino experiments find that the nature mass ordering is the inverted one, this could rule out a wide region in the currently allowed thermal axion window. Our results therefore strongly support multi-messenger searches of axions and neutrino properties, together with joint analyses of their expected sensitivities.