  The coevolution of galaxies and their metal content serves as an important
test for galaxy feedback models. We analyze the distribution and evolution of
metals within the IllustrisTNG simulation suite with a focus on the gas-phase
mass-metallicity relation (MZR). We find that the IllustrisTNG model broadly
reproduces the slope and normalization evolution of the MZR across the redshift
range $0<z<2$ and mass range $10^9 < M_*/\mathrm{M}_\odot < 10^{10.5}$. We make
predictions for the high redshift ($2<z<10$) metal content of galaxies which is
described by a gradual decline in the normalization of the metallicity with an
average high redshift ($z>2$) evolution fit by $\mathrm{d\;log(Z)}/\mathrm{dz}
\approx - 0.064$. Our simulations indicate that the metal retention efficiency
of the interstellar medium (ISM) is low: a majority of gas-phase metals ($\sim$
85 per cent at $z=0$) live outside of the ISM, either in an extended gas disk,
the circumgalactic medium, or outside the halo. Nevertheless, the redshift
evolution in the simulated MZR normalization is driven by the higher gas
fractions of high redshift galaxies, not by changes to the metal retention
efficiency. The scatter in the simulated MZR contains a clear correlation with
the gas-mass or star formation rate of the system, in agreement with the
observed fundamental metallicity relation. The scatter in the MZR is driven by
a competition between periods of enrichment- and accretion-dominated
metallicity evolution. We expect that while the normalization of the MZR
declines with redshift, the strength of the correlation between metallicity and
gas-mass at fixed stellar mass is not a strong function of redshift. Our
results indicate that the "regulator" style models are best suited for
simultaneously explaining the shape, redshift evolution, and existence of
correlated scatter with gas fraction about the MZR.
