The air–sea CO2 fluxes in eastern boundary upwelling systems (EBUSs) 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 EBUSs, i.e., the California (CalCS), Humboldt (HumCS), Canary (CanCS), and Benguela (BenCS) Current systems. In this study, we diagnose the climatic modes of the air–sea CO2 flux variability in these regions between 1920 and 2015, 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 (unforced) climate variability, whose contribution can be diagnosed by subtracting the ensemble mean from each simulation. We find that in the CalCS and CanCS, the resulting anomalous CO2 fluxes are strongly affected by large-scale extratropical modes of variability, i.e., the North Pacific Gyre Oscillation (NPGO) and the North Atlantic Oscillation (NAO), respectively. The CalCS has anomalous uptake of CO2 during the positive phase of the NPGO, while the CanCS has anomalous outgassing of CO2 during the positive phase of the NAO. In contrast, 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 large-scale 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 is key in improving our ability to predict the future evolution of the atmospheric CO2 source and sink characteristics of the four EBUSs.