Historical transboundary ozone health impact linked to affluence

Ozone pollution is a major transboundary threat to global health. Systematic improvement of mitigation strategy for transboundary ozone requires a socioeconomic understanding of historical lessons in countries at different affluence levels. Here, we explore the changes in transboundary ozone related premature deaths over 1951–2019 driven by anthropogenic emissions of four country groups categorized by income level. By integrating global emission datasets, a chemical transport model (GEOS-Chem), in situ ozone measurements worldwide and an ozone exposure-response model, we find that mortality caused by transboundary anthropogenic ozone increases by 27 times from 1951 to 2019, and on average contributes about 27% of global anthropogenic ozone related deaths. All groups exert and suffer from substantial transboundary ozone related mortality. The high-income and upper middle groups have each experienced an inverted U-shaped relationship between its affluence and per-million-people contribution to mortality caused by transboundary ozone, with the turning point around 23 000 USD and 6300 USD, respectively. The lower middle group has gradually matched the growth pathway of the upper middle group with a turning point less clear. Concerted efforts to ensure early turning points in less affluent countries will have considerable global health benefits.


Introduction
Ozone in the surface air is estimated to cause hundreds of thousands of premature deaths each year through worsening respiratory and cardiovascular diseases [1][2][3][4][5][6][7][8][9][10]. Formed from emitted precursor gases such as nitrogen oxides (NO x ) and nonmethane volatile organic compounds (NMVOCs) through highly nonlinear photochemistry, with a tropospheric lifetime of a few weeks, ozone can be transported long distances to other regions to exert transboundary impacts [11,12]. Despite decades of mitigation efforts especially in developed regions, ozone remains a major transboundary pollutant [11][12][13][14]. Attempts to systematically improve strategies for transboundary ozone control in the future must learn from common lessons in the past. Historically, changes in regional anthropogenic emissions of ozone precursors are determined by energy consumption [15][16][17], energy structure [18], technology and emission control levels [19,20], which are related to affluence level (indicated here by per capita gross domestic product, GDP). Yet it remains unclear how regional affluence is linked to its contribution to mortality associated with transboundary ozone.
Previous works have attempted to link regional environmental issues to its affluence level, usually based on local empirical data [17,[21][22][23][24]. A few studies have connected regional affluence levels to their anthropogenic ozone-related emissions [17,[22][23][24]. However, emission-based assessments do not take into account the strong nonlinear ozone chemistry [25,26]-with NMVOC emissions fixed, increasing NO x emissions enhances ozone formation in the NO x -limited chemical regime but decreases ozone in the NO x -saturated regime [27]. Other studies have linked regional affluence level to its ambient ozone pollution [28][29][30], with an implicit assumption that ozone is fully derived from local precursor emissions. These studies do not isolate the impacts of transboundary ozone originating outside that region's territory. For example, Asia is an important source region of ozone over North America, and vice versus [11,12,14,31,32]. The long-distance transboundary ozone means that ozone concentrations worldwide contributed by a region's emissions (as a source) differ from ozone concentrations over that region's territory (as a receptor). In addition to the complex ozone chemistry [27] and transboundary transport [11,12], the health response to ozone exposure may be nonlinear [3]. To date, a global-scale assessment of the potential connection between each region's transboundary ozone related health impact and its affluence level is lacking.
In this study, from a historical perspective, we quantify the impacts of regional anthropogenic emissions on global ozone-related premature mortality over 1951-2019, explore the role of transboundary pollution, and then reveal the linkage between regional transboundary impact and its affluence level. All countries are aggregated into four groups as emission source regions based on their per capita gross national income (GNI), as defined by the United Nations [33]. These groups include highincome (i.e. developed regions; with GNI per capita more than $12 235), upper middle income (e.g. China and Russia; with GNI per capita between $3956 and $12 235), lower middle income (e.g. India and Indonesia; with GNI per capita between $1006 and $3955) and low-income (GNI per capita less than $1005); see supplementary figure 1 of our previous study [34] for detailed regional disaggregation information. The classification of four groups is fixed, to explore the historical transboundary ozone related mortality impacts of country groups based on current worldwide economy. We only consider the impacts of emissions released in each group due to its production, and do not account for the emissions released in other groups associated with that group's consumption [35,36].

Method and datasets
This section describes our methodology, datasets and uncertainty estimates. Briefly, our study integrates historical emission data from the Community Emissions Data System (CEDS) [37] and the Multiresolution Emission Inventory for China (MEIC) [20,[38][39][40], the chemical transport model GEOS-Chem [41], in situ ozone measurements worldwide and an ozone exposure-response model [3]. We estimate the transboundary ozone related mortality caused by anthropogenic emissions in each income group for about every five year. We fix the meteorological conditions of GEOS-Chem at the 2014 level to focus on the sole impacts of anthropogenic emissions (section 2.2); a sensitivity test shows that using year-specific meteorology has a small effect on ozonerelated mortality results (supplementary figure 1). Simulated ozone concentrations by GEOS-Chem are evaluated and adjusted based on historical surface measurements (section 2.3). We employ Turner et al [3] ozone exposure-response model to estimate premature deaths of adults (⩾30 yr old) caused by respiratory and circulatory diseases associated with longterm exposure to ozone pollution. The model adopts a low concentration cutoff (LCC) of 26.7 ppbv as threshold below which ozone is assumed not to cause mortality [3]; choosing 0 ppbv as the LCC instead leads to a larger number of premature deaths and a strengthened relationship between regional affluence and its ozone morality impact (supplementary figure 2). Year-, country-and disease-specific baseline mortality rates since 1990 are calculated based on the Global Burden of Disease (GBD) 2019 database [42], and the rates at 1990 are applied to precious years [34]; a sensitivity test linearly extrapolating the rates to earlier years, following previous works [34,43], shows slightly different mortality results (by 2%-23% depending on the years, supplementary figure 3).

Historical population and GDP data
The historical country-based GDP and population data are obtained from the World Bank [44] for 264 regions, covering the period from 1960 to 2019. The GDP data are expressed in constant 2010 USD. When we calculate per capita GDP of each income group, we do not consider a few small countries of which GDP data in earlier years are missing.
To extend per capita GDP data to 1951, we employ exponential linear extrapolation. For each income group, we calculate the growth rate of per capita GDP from 1960 to 1970 (equation (1)), and then extend the group-based per capita GDP to 1951 (equation (2)) .
(2) G i,t p represents per capita GDP of country group i (high-income, upper middle income, lower middle income, or low-income) in the year t. The subscript p represents per capita. t1 denotes the year of 1960 and t2 the year of 1970. δ i represents the growth rate of per capita GDP of group i.

GEOS-Chem model simulations
A global chemical transport model GEOS-Chem version 11-01 [41] is employed in this study, to estimate the global ozone concentrations about every five years, including 1951, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2010, 2014 and 2019. All simulations cover 18 months, with the first six months used for model spin-up (e.g. the simulation for 2014 covers July 2013 through December 2014). The model setup follows our previous study [34]. Briefly, all simulations are driven by the MERRA2 assimilated meteorology provided by the NASA Global Modeling and Assimilation Office. The meteorological data are fixed in 2014 to quantify the sole impact of anthropogenic emissions. Supplementary figure 1 shows the mortality results using year-specific meteorology for 1981, which are similar to those using 2014 meteorology. The model is run with the full O x -NO x -VOC-CO-HO x gaseous chemistry, and aerosols are online calculated. The horizontal resolution is 2.5 • longitude × 2 • latitude grid with 47 vertical layers. Uptake of the hydroperoxyl radical (HO 2 ) on aerosols is accounted for with an uptake coefficient of 0.2, following previous works [14,45]. Vertical mixing of planetary boundary layer follows Lin and McElroy [46] Model convection uses the relaxed Arakawa-Schubert scheme [47]. Stratospheric ozone follows the Linoz scheme [48].
We use monthly gridded (0.5 • longitude × 0.5 • latitude) CEDS data from Hoesly et al [37] for gaseous (NO x , NMVOC, SO 2 , NH 3 ) and primary aerosol (BC and POA) pollutants globally over 1951-2014. CEDS emissions over China are replaced by emissions from MEIC [20,[38][39][40] since the year of 2000, because CEDS does not well represent the changes in Chinese emissions. Then we employ the country-and speciesspecific temporal changes over 2015-2019 in emission data from McDuffie et al [49] to derive the emissions over 1951-2019. NMVOC and NO x emissions of four groups are shown in supplementary figure 4. Supplementary material S1 and table S1 of our previous study [34] details more emission information, such as natural emission, international shipping and aviation emission inventories.
In each of the simulation years, we have six types of full-chemistry simulations to quantify surface ozone concentrations worldwide contributed by anthropogenic emissions of each country group. The control simulation (CTL) includes all anthropogenic emissions worldwide. Five sensitivity simulations zero out anthropogenic emissions worldwide and in four country groups group by group, including xANTH, xHIGH, xUPPER, xLOWER and xLOW. The difference between CTL and xHIGH represents the contribution of ozone by the high-income group, after the results are further adjusted to account for the nonlinear issue in ozone attribution (see below). The contributions of other groups are estimated in a similar way. International shipping emissions, aviation emissions and natural emissions are all included in xANTH simulations; we refer this part contribution as 'natural contribution' .
Considering the nonlinear relationship between ozone production and its precursors, we apply a linear weighted method to the model outputs on an hourly basis following a previous study [14]. For example, for the anthropogenic contribution of high-income group in each hour and grid cell, we first obtain the fraction of ozone contribution by this group to the total ozone summed based on the zero-out simulations: We then obtain the adjusted ozone contribution caused by the high-income group by applying f HIGH to the total surface ozone in CTL: Contributions of other groups are obtained similarly. The adjusted natural ozone is as follows:

GEOS-Chem evaluation and adjustment
We use ground-level ozone measurements worldwide to adjust for the effect of systematic model biases on our regional ozone attribution and health calculation. Hourly measurements in 2014 are taken from 1306 Environmental Protection Agency Air Quality System sites [50] for the US, from 216 sites of the National Air Pollution Surveillance program [51] for Canada, from 137 sites of European Monitoring and Evaluation Programme Chemical Coordination Centre network [52] for Europe, and from 945 sites of China National Environmental Monitoring Centre [53] for China. These data have been widely used for air quality assessment and other applications [54][55][56].
Since we use annual average daily maximum 8 h average (MDA8) ozone for health impact assessment (section 2.4), we further process the measurement data accordingly to facilitate comparisons with GEOS-Chem results. We convert the Chinese measurements from the original unit of µg m −3 at the standard condition (273.15 K and 1 atm) to mixing ratio (in ppbv). We further apply data quality control on all raw datasets, following a previous study [54]. We then convert the hourly data to MDA8, by calculating 8 h averages for 24 bins starting at 00:00 local time, and choose the highest of these 24 values. If there are less than six valid ozone values among each 8 h bin, the corresponding 8 h average value is considered missing. If there are two or more 8 h average values that are missing in a day, the corresponding MDA8 is considered missing. Lastly, we map all MDA8 data on a 2.5 • longitude × 2.0 • latitude grid for each day, with available data at all sites within a grid cell averaged to obtain the gridded MDA8 ozone.
We evaluate the annual average MDA8 ozone results of the CTL simulation. We first calculate the model MDA8 for each day and grid cell, and then select the model data only on days and in grid cells when the measured MDA8 data are available. As shown in the two columns on the left of supplementary figure 5, the CTL simulation overestimates the observed annual average MDA8 over China (by 9.8 ppbv on average), the US (by 8.9 ppbv), Europe (by 1.3 ppbv) and Canada (by 3.9 ppbv).
To remove the systematic model bias, we take a scaling approach in which a scaling factor (as the ratio of observed to CTL modeled annual mean MDA8 ozone) is derived for each range of population density (separated by yellow vertical lines in supplementary figure 6(a)). This population density-specific correction is used to facilitate the subsequent health impact calculations. It improves upon the approach by Shindell et al [2] that uses a single scaling factor for the whole globe. Here the scaling factors are derived as the average of the observed to modeled ratios in each population range, after combining all data in China, the US and Europe. Population data on a 0.1 • × 0.1 • grid are taken from GBD 2016 [57] health data and re-gridded to the model resolution (2.5 • longitude × 2.0 • latitude). Generally, the (observed to modeled) ratios are less than 1 and tend to decrease with increasing population density. The scatter in the ratios reflects local characteristics not captured by the model simulation.
Although not used for scaling the model results, the ratios derived based on data in each region are shown separately in supplementary figures 6(b)-(d) for comparison. Both the US-based and China-based ratios (supplementary figures 6(b) and (c)) are close to the ratios based on all data in the US, China and Europe together (supplementary figure 6(a)), although the Europe-based ratios tend to be larger reflecting the smaller model biases over Europe.
We apply the three-country based ratios (supplementary figure 6(a)) to the whole globe. After scaling, the modeled population-weighted global annual average MDA8 ozone (37.6 ppbv) matches the observed value (37.9 ppbv) (supplementary table 1). Supplementary figure 5 (the two columns on the right) further shows that the scaling greatly reduces the model overestimation at individual grid cells.
After scaling, the model bias is 1.3 ppbv averaged over China, 2.6 ppbv over the US, −4.2 ppbv over Europe, and −0.5 ppbv over Canada. The small bias for Canada indicates the robustness of our scaling approach, because the scaling factors are not derived based on Canadian measurements. To further test the robustness of our scaling approach over time, we employ the surface measurements over US and Europe in 2000. The results show our scaling approach effectively adjusts the modeled ozone to a small mean bias of 1.1 ppbv and −5.5 ppbv.
We then apply the scaling factors to all simulation years. For 1980 and earlier years when gridded GBD2016 population data are not available, we combine the gridded GBD2016 population data in 1990 with the historical country-based population time series from 'World Population Prospects' [58].
We further quantify the capability of GEOS-Chem in simulating the historical background ozone concentrations. For this purpose, we employ long-term ozone measurements at 12 background sites from tropospheric ozone assessment report (TOAR) [59]. These sites are categorized as 'background' sites, have records before 1980 and have data for at least five years. Sites that are not included here do not satisfy both criteria. Supplementary table 2 shows the geographical information of these sites. For these 12 sites, maximum daily 8 h average data were averaged to obtain monthly mean values. When we calculate the annual means, the corresponding simulated monthly values are set as missing if the observed monthly ones are missing. Given the scarcity of historical measurement data, for each site we averaged all yearly mean measurement data within ten year of each simulation year (1951, 1960, 1970, 1980, etc) to match the respective model simulation. Model values are selected at the grid cell and height of each site. The scatter plots in supplementary figure 7 show that GEOS-Chem reproduces the observed ozone both before and after the population-based adjustment, with a low bias and a slope close to 1 for each simulation year.

Premature mortality by ozone exposure
We use the pollution exposure-response relationship by Turner et al [3] to calculate the impact of annual average MDA8 ozone on premature mortality. This relationship is an update from Jerrett et al [60], and it has been widely used [1,2,5,61]. For premature mortality caused by each of respiratory and circulatory (including diabetes) diseases: Here, M g represents the number of premature deaths caused by each disease due to exposure to annual average MDA8 ozone in grid cell g. The global total mortality is summed over all grid cells. b g denotes the disease-specific gridded baseline mortality rate for people ⩾30 yr old, i.e. the fraction of total deaths caused by each disease. P g denotes the number of people ⩾30 yr old exposed to ozone at grid cell g. The country-, year-and age-specific population data are from the 'World Population Prospects' [58], and the gridded total population data are from GBD 2016 database [57]. The country-based age structure is applied to each grid cell within that country. The country-, disease-, age-and year-specific premature death data are from GBD 2019 data resources [42]. The baseline mortality rate of a country is calculated based on its population and mortality data, and is applied to all grid cells within that country. Since the GBD 2016 database provides mortality data since 1990, we apply the baseline mortality rates in 1990 to prior years.
β is the effect coefficient to present the association between annual average MDA8 ozone exposure and disease-specific premature mortality. It represents the natural logarithm of the Hazard Ratio (HR) for a 10 ppbv increase in ozone exposure. According to Turner et al [3], the circulatory (including diabetes) HR has a value of 1.03 (95% CI: 1.01-1.05), and the respiratory HR has a value of 1.12 (95% CI: 1.08-1.12).
∆O g is the difference between the GEOS-Chem simulated and bias-corrected annual average MDA8 ozone concentration and an LCC (low-concentration cutoff). Below the LCC, ozone is assumed to have no effect on mortality. We use choose an LCC value of 26.7 ppbv in the main analysis [3], following previous studies [1,2].
A previous study [62] shows that at there is no significant difference in ozone-related health impacts at lower end of observed exposures compared to those at higher exposures. Therefore, we also include health impact estimation without an exposure threshold, that is, choosing 0 ppbv as the LCC instead. The results are shown in supplementary figure 2. Ozone pollution caused by anthropogenic emissions drops below the LCC decades ago, which leads to a sharp decrease of ozone-related premature deaths, such as the sharp change of mortality impacts of upper middle group with its per capita GDP around $1000 as shown in figure 4(a). Choosing 0 ppbv as the LCC narrows the change in the per-million-people contributions to mortality impact of upper middle group when its per capita GDP is around $1000 as shown in supplementary figure 2. Therefore, not only the number of premature deaths get larger, but also the relationship between regional affluence and its ozone mortality impacts gets strengthened.
For each year, we first calculate the total mortality (M g ) due to annual average MDA8 ozone from all anthropogenic and natural sources (after ozone bias corrected; equation (6). We then calculate the number of deaths contributed by each income group by multiplying M g by the fraction of ozone contributed by that group. Equation (9) provides an example for the high-income group. Mortality contributions of other groups are obtained similarly. Such a direct proportion approach has been used in many previous studies [34,35,63,64] O g = HIGH g + UPPER g + LOWER g + LOW g Supplementary figure 8 compares our global total mortality results to other studies. The spread in mortality outcomes reflects the large uncertainty in the current knowledge of the ozone mortality impact. Our results in 2010 are close to the recent study by Shindell et al [2] (with ozone bias corrected and with Turner et al [3] exposure-response function) on a global total basis.

Decomposition of changes in ozone-related mortality
Following our previous study [34], historical changes in anthropogenic ozone-related mortality is decomposed into three factors, including baseline mortality rate (b), ozone-related term (T = 1 − e −β∆O ) and exposed population (P), according to the calculation formula in section 2.4.
The change in anthropogenic ozone-related mortality between two simulation years (∆M) could be expressed as: Contributions of one factor with the rest two kept constant are denoted by three terms on the right-hand side. There are six possible decompositions of these three terms, and the averaged decomposition results is used in the main text. See more detailed information in our previous study [34].

Uncertainty and limitation
This study is subject to uncertainties from a few sources. First, the historical CEDS [37,49] and MEIC [20,[38][39][40] inventories used here contain errors ranging from 10% to 170% depending on the pollutant and region, which is further discussed in the limitation discussion below. For our ozone simulations and mortality impact estimates, there errors are embedded in the derivation of the GEOS-Chem model error (σ1).
Second, we use the bias between model simulation and measurements to represent the error of GEOS-Chem, following previous studies [34,63,64]. Supplementary table 2 shows the mean relative model bias after bias adjustment over mainland China, US, Europe, Canada and globe (combining the four regions together). We use the relative bias over Europe (−14.1%, the worst case) to represent the overall model error. The error is referred to as σ1 (one standard deviation).
Third, the ozone exposure-response model by Turner et al [3] is subject to errors in linking pollution exposure, specific diseases and premature mortality. We further discuss this as one limitation below. As shown in section 2.4, we use the 95% CI of the HR for respiratory and circulatory diseases to characterize errors in the exposure-response model. The overall error is about 30% and referred to as σ2 (one standard deviation).
The overall uncertainty in the historical premature mortality calculation is estimated as the sum in quadrature of σ1 and σ2. The majority of these error terms do not change with time, thus the temporal changes (in a relative sense) in mortality impacts are less affected by these errors.
Besides the uncertainty sources mentioned before, this study is subjected to several limitations. First, the anthropogenic emission inventories employed in this study tend to contain larger uncertainties in the early years and over some less affluent regions. Therefore, we use other two widely used global anthropogenic emission inventories to explore the linkage between regional affluence level and its ozone precursors' emissions. As shown in supplementary figure 9, the patterns are similar to those in figures 3(c) and (d).
Second, simulated ozone concentrations are all obtained using GEOS-Chem. Although we have employed global surface ozone observations and long-term background ozone observations to validate and adjust the simulated results, further mutimodel analyses would be helpful to estimate the intrinsic biases embedded in GEOS-Chem. Although we conduct consistent quality control processes to all raw observation datasets in this study, the monitoring networks over different regions have different internal quality control standards and internal procedures [65,66]. Besides, the surface ozone observations are lacking over less affluent regions [66] (such as rural areas, and globally such as Africa), which could potentially introduce classification error when quantifying ozone exposures [67]. Therefore, we call for improvement in observation metadata reporting, more consistent monitoring procedures in global monitoring networks and more studies focusing on observations and emission datasets over less affluent regions in future to better quantify the uncertainty.
Third, the application of ozone health model is subjected to errors in the internal validity and the generalizability of the health model [67][68][69]. A full evaluation of this error is technically prohibitive [67]. Following previous studies [1,2,34], we estimate the uncertainty ranges, and further conduct sensitivity tests to better address its limitation. Besides the ozone health model from Turner et al [3], we also employ another risk function from Jerrett et al [60] to estimate the ozone-related premature deaths with similar results to a previous study [1] (supplementary figure  8). We chose 0 ppbv as LCC instead of 26.7 ppbv in another test, which leads to a large number of ozonerelated premature deaths and a strengthened relationship between regional affluence and its contribution to mortality caused by transboundary ozone (supplementary figure 2). The baseline mortality rates, as the input datasets, are fixed at 1990 for earlier years, due to data limitation. To test this corresponding effect, we employ linear-extrapolated countrybased baseline mortality rates to earlier years and find slightly different mortality results (supplementary figure 3). Further ozone-related cohort studies based on various races, age groups and other socio-economic conditions in future will help to reveal more detailed in transboundary ozone related mortality.  contribute about 71% of the total ozone-related mortality, is not elaborated further in this study.

Aggravating ozone health impact
To better understand the drivers of the increase in anthropogenic ozone mortality impact, we attribute the mortality growth to three factors using a decomposition analysis (section 2.5). The three driving factors considered are ozone concentration, population, and baseline mortality rate. Globally, the baseline mortality rate declines continuously as a result of improved medical and living conditions, based on available data since 1990. However, the increasing ozone concentration and population overcompensate for the beneficial effect of declining baseline mortality rate and enhance the anthropogenic ozone health burden. The worsening ozone exposure alone causes a 16-fold increase in the total mortality increase from 1951 to 2019, and the population change causes a 24fold increase over the same period.
The mortality contribution by ozone from each of the four groups has changed substantially over the years. The high-income group contributes as much as 77% of the global total premature mortality in 1951. In absolute terms, ozone from the high-income group leads to 6 From 1980 to 1990, the ozone concentration contributed by the high-income group declines with its decreasing precursor emissions, but this effect is offset by the growth in global population exposed to ozone pollution originating from the high-income group. Although the mortality caused by the highincome group declines from 1990 and 2010, there is slight growth (by about 1%) from 2010 to 2019. This resumption of mortality growth is primarily because of worsened ozone pollution originating from Europe through the nonlinear chemistry driven by its declining NO x emissions (supplementary figure  4(a)). The recent worsening European ozone pollution is also found by previous studies based on surface measurements [56,[70][71][72]. The complex effects of emissions on ozone concentrations and associated mortality highlight the challenges to mitigate ozone pollution.
For the other three groups, their impacts on global ozone mortality have continued to increase since 1951, as a result of general increases in attributable ozone concentration (due to increasing emissions, supplementary figure 4) and population overcompensating the effect of declining baseline mortality rates. As an exception, the decreases in ozone exposure caused by the upper middle and lower middle income groups between 1990 and 2000 are due to significant drop in emissions in many countries [37,49] which used to be parts of the Soviet Union. By 2019, the premature deaths caused by the upper middle, lower middle and low-income groups reach 94. ] thousand, respectively. In other words, by 2019 the lower middle group becomes the top contributor to annual global premature mortality (53%), followed by the upper middle (34%), high-income (10%) and lowincome group (3%). Figure 1(b) shows inter-decadal changes in annual premature deaths attributable to transboundary anthropogenic ozone (caused by four source groups but occurring outside their territories) and contributions from the three driving factors (ozone, population and baseline mortality rate) from 1951 to 2019. Results for every five years are shown in supplementary figure 10(b) The worsening ozone pollution alone leads to a 10fold increase, the rising population causes a 21-fold increase, whereas the declining baseline mortality rate slightly offsets these increases. Between 1951 and 2019, on average, premature deaths caused by transboundary anthropogenic ozone account for about 27% of annual global anthropogenic ozone related deaths.

Substantial mortality caused by transboundary anthropogenic ozone
Ozone originating from the high-income group contributes as much as 84% of annual global mortality caused by transboundary anthropogenic ozone in 1951 and more than 50% before 1990. In absolute terms, the annual contribution of high-income group peaks in 1990 with 17.5 [95% CI: 9.0-26] thousand deaths. The contributions from other three groups increase continuously from 1951 to 2019. In recent decades, the percent contributions of highincome and low-income groups to mortality caused by transboundary ozone are larger than their percentage contributions to total (local + transboundary) mortality ( figure 1(b) versus 1a). In 2019, the highincome group contributes about 25.7% of deaths caused by transboundary ozone ( figure 1(b)) but only 10% of total ozone related deaths ( figure 1(a)), and the respective numbers for the low-income group are 10.5% and 3%. The contributions of upper middle and lower middle groups are about 34.4% and 29.5% in 2019, lower than their contributions to total deaths. Figure 2 compares the transboundary impact caused by and exerted upon each group, by contrasting the ratio of deaths outside a region's territory and global deaths caused by that region (x-axis, as a source) against the ratio of deaths in a region caused by foreign pollution to that region's total deaths (y-axis, as receptor). The fractional transboundary impacts both caused by and exerted upon the highincome group are increasing over the past decades to reach 0.6 and 0.3 in 2019, with the former (as a source) always greater than the latter (as a receptor). The fractional transboundary impacts caused by and exerted upon the upper middle and lower middle groups decline with time, as a result of their growing local fractional contributions. The fractional impacts for the low-income group remain at high values of 0.6-0.9 throughout the years, reflecting stable, strong tie of its ozone problem with other groups. Figure 3 shows the transboundary ozone related mortality inter-impacts between the four groups, in terms of the impact of every million people in a source region, for every decade.  (over 1960-1975) or even higher than (in 1951) the self-impacts of upper middle and lower middle groups. Figure 3 also shows that the mortality impacts attributable to transboundary ozone caused by every  . This is due to the regional differences in population exposed to ozone and in baseline mortality rate. The same phenomenon exists between the upper middle and the lower middle groups. Figure 4 links the transboundary impact by every million people in a given group (as a source region) to its affluence level (i.e. per capita GDP), by putting together data in all groups and years. This allows to evaluate the evolution of regional transboundary ozone impact in association with affluence from a macroscopic, socioeconomic perspective. We examine the per-million-people contributions of each source region to its transboundary ozone related fractional mortality impact (FMI, the ratio of mortality number to total population; figure 4(a)), to transboundary population-weighted MDA8 ozone ( figure 4(b)), and to its own NO x (figure 4(c)) and NMVOC (figure 4(d)) emissions. Analysis of transboundary FMI caused by every million people in each source region removes the influence of historical population growth. Results for local FMI are shown in supplementary figure 12 for comparison.

Linkage between mortality attributable to transboundary ozone and affluence
For the high-income group, a clear inverted U-shaped relationship exists between its affluence level and per-million-people transboundary FMI ( figure 4(a)). The transboundary FMI increases with affluence until a turning point, after which the transboundary FMI decreases substantially with further increasing affluence. The per capita GDP at the turning point is about 23 000 USD. Similar inverted U-shaped patterns are evident for this group in per-million-people contributions to transboundary ozone, NO x emissions and NMVOC emissions (figures 4(b)-(d)).
For the upper middle group, its per-millionpeople NO x emissions peak at an affluence level of about 6300 USD, which is much lower than the turning point of high-income group. The lower turning point is partly due to emission control measurements in China [20,37,49], which belongs to this group. Prior to this turning point, the upper middle group follows the growth pattern of per-million-people NO x emissions of the high-income group. Per-millionpeople NMVOC emissions of the upper middle group deviate from the path of high-income group and start to decline at an affluence level around 7500 USD.
For the upper middle group, its per-millionpeople transboundary FMI is very small when the affluence level is below 1000 USD, but grows substantially with affluence to quickly match the path of high-income group until a turning point similar to those for precursor emissions. The small mortality impact at low-affluence situations is because ozone concentrations are below the LCC (26.7 ppbv) in many downwind regions of the upper middle group. There is emerging evidence that the health impact of ozone may not have such an exposure threshold [62]. Removing this threshold would substantially increase the mortality impact at low-affluence situations (supplementary figure 2).
Per-million-people NO x and NMVOC emissions and transboundary ozone impact of the lower middle group have followed similar paths as the upper middle group, albeit with some differences at low affluence situations around 1000 USD. For the per-millionpeople transboundary FMI, the gap between the paths of lower middle and upper middle groups narrows down along with the increasing affluence level.
Among the four source regions, the low-income group has undergone the weakest growth in affluence in both absolute and relative terms. For this group, there is substantial scattering in its per-millionpeople emissions, transboundary ozone and mortality impact. The relatively small changes in affluence render an affluence-pollution relationship more difficult to be established.

Conclusion
This study shows considerable transboundary ozone related premature mortality over the course of 1951-2019 and its evolution with regional affluence. Bidirectional transboundary pollution between less and more affluent regions are substantial. The highincome group has undergone an inverted U-shaped pattern between its affluence and its mortality impact caused by transboundary ozone, after removing the effect of historical population growth in both source and receptor regions. The upper middle group has experienced a turning point at a lower affluence level, whereas turning points are not obvious in the lower middle and low-income groups.
As less affluent countries try to achieve sustained economic growth, they continue to suffer from increasing pollutant emissions in general. Starting at a low-affluence turning point to invert this emission trend will have local and nonlocal health benefits via atmospheric transport. This task is possible, because the scientific knowledge and financial and technical capabilities for ozone control nowadays are much better than decades ago [73]; the turning point of the upper middle group at an affluence level much lower than the high-income group is strong evidence. Although the short-term economic costs of pollution mitigation can be substantial for less affluent regions, such burden can be compensated by its health benefits [2,5,74]. Stringent mitigation policies may also provide new opportunities for technological innovation [19,75,76], efficiency improvement [77,78], and labor demand enhancement [79,80]. Measures to cut pollutant emissions often reduce carbon dioxide simultaneously with substantial cobenefits for climate change mitigation.
Besides domestic actions by less affluent countries to reduce emissions, international financial and technological aids from more affluent countries would lead to a win-win health outcome. This is important especially given that large fractions of present emissions in less affluent countries are associated with export of goods to supply consumption worldwide [35,36] in a globalized economy. To this end, our study offers a socioeconomic insight into historical inter-regional pollution linkage to facilitate concerted mitigation action.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.