The understanding of long-time-series variations in air-sea CO2 flux in the Bering Sea is critical, as it is the passage area from the North Pacific Ocean water to the Arctic. Here, a data-driven remote sensing retrieval method is constructed based on a large amount of underway partial pressure of CO2 (pCO2) data in the Bering Sea. After several experiments, a Gaussian process regression model with input parameters of sea surface temperature, sea surface height, mixed-layer depth, chlorophyll a concentration, dry air mole fractions of CO2, and bathymetry was selected. After validation with independent data, the root mean square error of pCO2 was< 24 μatm (R2 = 0.94) with satisfactory performance. Then, we reconstructed the sea surface pCO2 in the Bering Sea from 2003 to 2019 and estimated the corresponding air-sea CO2 fluxes. Significant seasonal variations were identified, with higher sea surface pCO2 in winter/spring than in summer/autumn in both the basin and shelf area. Semiquantitative analysis reveals that the Bering Sea is a non-temperature-dominated area with a mean temperature effect on pCO2 of 12.7 μatm and a mean non-temperature effect of −51.8 μatm. From 2003 to 2019, atmospheric pCO2 increased at a rate of 2.1 μatm yr−1, while sea surface pCO2 in the basin increased rapidly (2.8 μatm yr−1); thus, the CO2 emissions from the basin increased. However, the carbon sink in the continental shelf still continuously increased. The whole Bering Sea exhibited an increasing carbon sink with the area integral of air-sea CO2 fluxes increasing from 6 to 19 TgC over 17 years. Meanwhile, the seasonal amplitudes in pCO2 in the shelf area also increased, approaching 14 μatm per decade. The reaction of the continuously added CO2 in continental seawater reduced the ocean CO2 system capacity. This is the first study to present long-time-series satellite data with high resolution in the Bering Sea, which is beneficial for studying the changes in ocean ecosystems and carbon sink capacity.