Within the framework of transverse momentum dependent factorization in combination with nonrelativistic quantum chromodynamics (QCD), we study charmonium and bottomonium production in hadronic collisions. We focus on quarkonium states with even charge conjugation, for which the color-singlet production mechanism is expected to be also dominant in the small transverse momentum region, qT2≪4Mc,b2. It is shown that the distributions of linearly polarized gluons inside unpolarized, longitudinally, and transversely polarized protons contribute to the cross sections for scalar and pseudoscalar quarkonia in a very distinctive, parity-dependent way, whereas their effects on higher angular momentum states are strongly suppressed. We derive analytical expressions for single and double spin asymmetries, which would allow for the direct extraction of the gluon transverse momentum dependent distributions, mirroring the phenomenological studies of the Drell-Yan processes aimed at the extraction of their quark counterparts. By adopting Gaussian models for the gluon transverse momentum dependent distributions, which fulfill without saturating everywhere their positivity bounds, we provide numerical predictions for the transverse single-spin asymmetries. These observables could be measured at LHCSpin, the fixed target experiment planned at the LHC. Published by the American Physical Society 2024