The air-sea CO2 fluxes in Eastern Boundary Upwelling Systems (EBUS) vary strongly in time and space with some of the highest flux densities globally. The processes controlling this variability have not yet been investigated consistently across all four major EBUS, i.e., the California (CalCS), Humboldt (HumCS), Canary (CanCS), and Benguela (BenCS) Current Systems. In this study, we diagnose the physical and biological mechanisms that contribute to historical (1920–2015) CO2 flux variability in these regions using simulation results from the Community Earth System Model Large Ensemble (CESM-LENS), a global coupled climate model ensemble that is forced by historical and RCP8.5 radiative forcing. Differences between simulations can be attributed entirely to internal climate variability. We find that the deviations from the ensemble mean, i.e., the anomalous CO2 fluxes, in the CalCS and CanCS are strongly affected by modes of variability associated with atmospheric subtropical gyres: the North Pacific Gyre Oscillation (NPGO) and the North Atlantic Oscillation (NAO), respectively. The CalCS (CanCS) has anomalous uptake (outgassing) of CO2 during the positive phase of the NPGO (NAO). The HumCS is mainly affected by El Niño Southern Oscillation (ENSO), with anomalous uptake of CO2 during an El Niño event. Variations in dissolved inorganic carbon (DIC) and sea surface temperature (SST) are the major contributors to these anomalous CO2 fluxes, and are generally driven by changes to gyre circulation, upwelling, the mixed layer depth, and biological processes. A better understanding of the sensitivity of EBUS CO2 fluxes to modes of climate variability may improve our ability to predict the ocean–atmosphere carbon cycle in EBUS, which are particularly susceptible to ocean acidification.