Realistic model representation of ocean phytoplankton is important for simulating nutrient cycles and the biological carbon pump, which affects atmospheric carbon dioxide (pCO2) concentrations and, thus, climate. Until recently, most models assumed constant ratios (or stoichiometry) of phosphorous (P), nitrogen (N), silicon (Si), and carbon (C) in phytoplankton, despite observations indicating systematic variations. Here, we investigate the effects of variable stoichiometry on simulated nutrient distributions, plankton community compositions, and the C cycle in the preindustrial (PI) and glacial oceans. Using a biogeochemical model, a linearly increasing P:N relation to increasing PO4 is implemented for ordinary phytoplankton (PO), and a nonlinearly decreasing Si:N relation to increasing Fe is applied to diatoms (PDiat). C:N remains fixed. Variable P:N affects modeled community composition through enhanced PO4 availability, which increases N-fixers in the oligotrophic ocean, consistent with previous research. This increases the NO3 fertilization of PO, the NO3 inventory, and the total plankton biomass. Surface nutrients are not significantly altered. Conversely, variable Si:N shifts south the Southern Ocean’s meridional surface silicate gradient, which aligns better with observations, but depresses PDiat growth globally. In Last Glacial Maximum simulations, PO respond to more oligotrophic conditions by increasing their C:P. This strengthens the biologically mediated C storage such that dissolved organic (inorganic) C inventories increase by 34-40 (38-50) Pg C and 0.7-1.2 Pg yr-1 more particulate C is exported into the interior ocean. Thus, an additional 13-14 ppm of pCO2 difference from PI levels results, improving model agreement with glacial observations.