We use the Color Glass Condensate (CGC) framework to study the production of forward heavy quark-antiquark pairs in unpolarized proton-nucleus or proton-proton collisions in the small-x regime. In the limit of nearly back-to-back jets, the CGC result simplifies into the transverse-momentum dependent (TMD) factorization approach. For massless quarks, the TMD factorization formula involves three unpolarized gluon TMDs: the Weizs{a}cker-Williams gluon distribution, the adjoint-dipole gluon distribution, and an additional one. When quark masses are kept non-zero, three new gluon TMDs appear -- each partnered to one of the aforementioned distributions -- which describe the distribution of linearly-polarized gluons in the unpolarized small-x target. We show how these six gluon TMDs emerge from the CGC formulation and we determine their expressions in terms of Wilson line correlators. We calculate them analytically in the McLerran-Venugopalan model, and further evolve them towards smaller values of x using a numerical implementation of JIMWLK evolution.