Recent results form the James Webb Space Telescope (JWST) report space densities for bright and massive galaxies at z>7 that far exceed expectations of theoretical models of galaxy formation, prompting a revision of our understanding of the physical processes leading to the assembly of the first luminous structures. In this work we present predictions from a realization of the GAlaxy Evolution and Assembly (GAEA) model, which implements a prescription for a variable stellar initial mass function (IMF). This prescription is inspired by high-resolution numerical simulations that account for the role of cosmic rays (CR) as regulators of the star formation rate (SFR) in giant molecular clouds. In our approach, SFR density is assumed to be a proxy for the CR density, providing a link between the IMF shape and the predicted physical conditions of the star forming interstellar medium. Our results show that, in our model framework, assuming such a variable IMF reproduces several properties of the z>6 galaxy population, with no further modification of the feedback model, including their UV luminosity functions up to z~13. In order to compare model predictions with available estimates for the galaxy stellar mass function (GSMF), we reconstruct stellar masses from the model's synthetic photometry assuming a universal IMF, reflecting standard observational practice. Under this approach, we show that the model can reproduce the evolution of the GSMF up to the highest redshifts accessible. Our findings highlight the need to consider a variable IMF shape in the error budget associated with stellar mass estimates. We show that the evolution of both the slope and normalization of the gas-phase mass metallicity relation can be used as powerful discriminant between models of early galaxy formation assuming different IMF evolution.
