Determining the shape of the stellar initial mass function (IMF) and whether it is constant or varies in space and time is the Holy Grail of modern astrophysics, with profound implications for all theories of star and galaxy formation. On a theoretical ground, the extreme conditions for star formation (SF) encountered in the most powerful starbursts in the Universe are expected to favour the formation of massive stars. Direct methods of IMF determination, however, cannot probe such systems, because of the severe dust obscuration affecting their starlight. The next best option is to observe CNO bearing molecules in the interstellar medium at millimetre/ submillimetre wavelengths, which, in principle, provides the best indirect evidence for IMF variations. In this contribution, we present our recent findings on this issue. First, we reassess the roles of different types of stars in the production of CNO isotopes. Then, we calibrate a proprietary chemical evolution code using Milky Way data from the literature, and extend it to discuss extragalactic data. We show that, though significant uncertainties still hamper our knowledge of the evolution of CNO isotopes in galaxies, compelling evidence for an IMF skewed towards high-mass stars can be found for galaxy-wide starbursts. In particular, we analyse a sample of submillimetre galaxies observed by us with the Atacama Large Millimetre Array at the peak of the SF activity of the Universe, for which we measure 13C/18O⋍1. This isotope ratio is especially sensitive to IMF variations, and is little affected by observational uncertainties. At the end, ongoing developments of our work are briefly outlined.