Standard analytical chemical evolution modelling of galaxies has been assuming the stellar initial mass function (IMF) to be invariant and fully sampled allowing fractions of massive stars to contribute even in dwarf galaxies with very low star formation rates (SFRs). Recent observations show the integrated galactic initial mass function (IGIMF) of stars, i.e. the galaxy-wide IMF, to become systematically top-heavy with increasing SFR. This has been predicted by the IGIMF theory, which is here used to develop the analytical theory of the chemical evolution of galaxies. This theory is non-linear and requires the iterative solution of implicit integral equations due to the dependence of the IGIMF on the metallicity and on the SFR. It is shown that the mass-metallicity relation of galaxies emerges naturally, although at low masses the theoretical predictions overestimate the observations by 0.3--0.4 dex. A good agreement with the observation can be obtained only if gas flows are taken into account. In particular, we are able to reproduce the mass--metallicity relation observed by Lee et al. (2006) with modest amounts of infall and with an outflow rate which decreases as a function of the galactic mass. The outflow rates required to fit the data are considerably smaller than required in models with invariant IMFs.