Life expectancy refers to the mean years of life remaining at a given age through the mortality rates observed at a given year [23]. This study considers the gap in gender-specific life expectancy (GGLE) as the dependent variable. Specifically, the value of GGLE is calculated as each country’s female life expectancy minus its male life expectancy (Table 1).

We calculate the population-weighted PM_2.5 (Table 1) at the country level as a proxy variable to quantify the magnitude of air pollution in the exposed population on the country scale by using the following formula: p w P M 2.5 = ∑ i = 1 n P M 2.5 i × p o p i p o p (1)

In Eq. 1, p w P M 2.5 represents the population-weighted PM_2.5 for a country; P M 2.5 i refers to the value of ground-level PM_2.5 in the ith grid. p o p i is the population count in the ith grid; pop is the total population at country level. The annual concentrations of ground-level PM_2.5 in the dataset are based on simulated aerosol optical depth [24].

The urbanization rate was found to be closely positively related to LE based on previous studies [25, 26]. Therefore, the urbanization rate is selected as one of the core independent variables in this study.

Several covariates were selected in our model, including poverty level, education level, nutrition level, and the death rate from smoking, which showed significant relations to LE based on existing studies. Poverty plays an important role in LE determinants [27, 28]. In this study, the poverty gap index was collected to reflect the level of poverty, which refers to the mean shortfall of consumption/income from the poverty line. Education level, especially for female education conditions, is an important factor in the LE [29, 30]. In this paper, we collected primary education level as the indicator to measure the education level, which refers to the female pupils as a percentage of total pupils at the primary level, including enrollments in public and private schools. Daily caloric supply refers to the average per capita caloric availability. This indicator indicates the caloric availability delivered to households but does not necessarily indicate the number of calories consumed. A prior study showed that food calories could increase LE [31]. The death rate from smoking (both sex and age-standardized) refers to the annual deaths attributed to smoking (per 100,000 people). Previous papers indicated the negative effects of smoking on LE [32–34].

In this study, the data on life expectancy for females and males in each country from 1960 to 2018 are obtained from the World Development Indicators (WDIs) database published by the World Bank [35]. Further, the annual concentrations of ground-level PM_2.5 (with dust and sea salt removed) data for 1998–2016 are collected from the Global Annual PM_2.5 Grids data with the spatial resolution of 0.01°×0.01° released by EARTHDATA (https://sedac.ciesin.columbia.edu/data/set/sdei-global-annual-gwr-pm2-5-modis-misr-seawifs-aod). We obtained the gridded population of the world at 1,000 m spatial resolution for the year 2005, according to the Socioeconomic Data and Applications Center of NASA (https://sedac.ciesin.columbia.edu/data/collection/gpw-v4). Data on the share of urban population (% of the total population) is from the World Urbanization Prospects: 2018 Revision released by the United Nations Population Division. The data of the poverty gap index is from the Poverty and Inequality Platform by the World Bank, which is based on primary household survey data (https://data.worldbank.org/indicator). Data on the death rate from smoking (both sex with age-standardized) is from the Global Burden of Disease Collaborative Network by the Institute for Health Metrics and Evaluation (IHME) (http://ghdx.healthdata.org/gbd-results-tool). Data on daily caloric supply is from the Our World in Data (OWID) based on UN FAO & historical sources (https://ourworldindata.org/calorie-supply-sources). Data on education are collected from the database of the UNESCO Institute for Statistics (http://uis.unesco.org/).

To model GGLE, a series of Bayesian spatiotemporal models are proposed (Supplementary Table S1), where for each country (i = 1, 2, ..., 134) and year step (t = 1960,1961, ..., 2018), the GGLE, y _it, was modeled by y i t ∼ N o r m a l   μ i t , σ e 2 . Here, μ i t denotes the expected value of y i t and σ e 2 measures the variance of y i t . The μ i t could be calculated as the following formula (model M0s): μ i t = α + s i + θ t + δ i t (2)

In Eq. 2, α is an intercept that measures the overall y i t during the study period (1960–2018); The spatial term s i , denotes the spatial random effect capturing the spatial dependency of y i t . We note that s i captures the overall spatial random effects common from 1960 to 2018, describes the difference between y i t in the i-th country or region relative to the global average level, α . When s i >0, it indicates that the y i t in country i is higher than overall y i t across the whole study period. The term θ t denotes a dynamic temporal trend that captures the overall temporal trend common to all countries. The δ i t is a space-time interaction random effect, representing a vector that varies through space and time [36]. The term δ i t allows spatial pattern changes from one time frame to another and temporal trend varies from one country to another.

Furthermore, a Bayesian spatiotemporal ecological regression was developed to investigate the impacts of air pollution and urbanization on GGLE and gender-specific life expectancy while adjusting for confounding factors and spatiotemporal variation. As an extension of model M0s, the variables of PM and UP (two interested variables) were considered, adjusting for confounding factors like poverty (P), education (E), calorie supply (C), and smoking (S). The new model was defined as the following formula (model M1s): y i t = α + s i + θ t + δ i t + β 1 ln ⁡ P M i t + β 2 ln ⁡ U P i t + β 3 ln ⁡ P i t + β 4 ln ⁡ E i t + β 5 ln ⁡ C i t + β 6 ln ⁡ S i t (3)

In Eq. 3, β 1 , β 2 , ..., β 6 refer to the corresponding regression coefficients of pwPM_2.5, urbanization, poverty, education, supply of calories, and smoking, respectively, measuring the effects of these six variables across 134 overall countries on y (GGLE, LEm, and LEf). A previous study found that the effects of air pollution and urbanization on LE and gender-specific LE may be spatially varied [37]. Model M1s is modified to allow each region to present its own regression coefficient to examine the spatially varied impacts of pwPM_2.5 and urbanization on GGLE and gender-specific LE. In this paper, a total of 134 countries are divided into six regions based on the definition of the World Health Organization (WHO) (i.e., Africa, Americas, Eastern Mediterranean, Europe, South-East Asia, and Western Pacific), to investigate the region-specific effects of pwPM_2.5 and urbanization on GGLE and gender-specific life expectancy, the modified regression model is shown as follow (model M2s): y i t = α + s i   θ t + δ i t + ∑ j = 1 n β 1 j ln ⁡ P M i t + ∑ j = 1 n β 2 j l n U P i t + β 3 ln ⁡ P i t + β 4 ln ⁡ E i t + β 5 ln ⁡ C i t + β 6 ln ⁡ S i t (4)

In Eq. 4, β 1 j and β 2 j represent the regression coefficients of influencing factors in a different region, j = 1, 2, …, n (n = 6). For detailed methods, please refer to the supplementary materials.

The average GGLE at the country level and the world is illustrated from 1960 to 2018 (Figure 1). Some key points can be drawn. Firstly, we could find the spatial difference across regions. Europe and North American countries show the highest peak values of GGLE, which shows a trend of rising first and then falling in GGLE. Secondly, other regions show an increasing trend during the study period, similar to the global trend. Furthermore, GGLE in South Asia (e.g., India) shows an increasing trend from negative to positive values, whereas other countries show positive values of GGLE during the study period. In addition, the GGLE in the Middle East and North Africa presents an anomaly peak during the 1980s, possibly related to the sharp decreased male LE values in Iran and Iraq.

Further, most countries show a trend from increase to decrease of GGLE with peak values, 41 countries have an increasing trend without a peak value, and only 12 countries/regions show a decreasing trend without a peak value during the study period. The difference is depicted in space, indicating that the countries with increasing or decreasing trends are mainly concentrated in South Asia, Middle Africa, the Middle East, and Latin America (Figure 1). Countries with peak values also present spatial differences. For example, Russia and countries in East Europe, Central Asia, Middle East show the higher level of peak values of GGLE, and countries in Western Europe, America, South Africa, Japan, South Korea, Australia, and New Zealand show the middle range of peak values of GGLE, while countries in North Africa, East-southeastern Asia show the lower peak values of GGLE (Figure 1).

The combined spatial effects of s i + θ t + δ i t at four different time points (1960, 1980, 2000, and 2018) for GGLE were estimated using the Bayesian spatiotemporal modeling, which depicts a spatial pattern and its changes during the study period 1960–2018 (Figure 2). Positive values of the combined spatial effects indicate that GGLE in country i at time point t is higher than the overall level, while negative values indicate that lower than the overall level. We found that the spatial pattern of GGLE slightly changed from 1960 to 2018. However, the whole spatial pattern is similar at these four-time points. For example, GGLE in Russia is consistently higher than the overall level. In contrast, most countries in Africa have lower GGLE than the overall level. China has a substantial change in GGLE from 1960 to 2018, which indicates a lower level in 1960, 1980, and 2000, but higher levels in 2018 compared to the overall level in the world. Other countries like the U.S., Japan, South Korea, Brazil, Argentina, and South Africa show that GGLE is consistently higher than the overall level.

The overall temporal trend of GGLE ( θ t ) is estimated using the Bayesian spatiotemporal model from 1960 to 2018. The overall temporal trend of GGLE presents a clear nonlinear temporal trend in the world (Figure 3). The GGLE increased rapidly from 1960 to 1985, and then grew slowly down but maintained an increasing trend. From 1995 to 2005, the GGLE decreased rapidly by a clear linear temporal trend. After 2005, the GGLE increased again. Moreover, we estimated the country-specific temporal trends ( θ t + δ i t ) of GGLE to explore whether there are huge differences across countries worldwide. Figure 3 shows temporal trends of GGLE in four selected countries from 1960 to 2018, including China, the U.S., Russia, and India. China slightly increases GGLE, while GGLE in the U.S. shows a nonlinear temporal trend. However, there are no clear significant temporal trends of GGLE in China and the U.S. during the study period. Figure 3 indicates that GGLE in Russia has a stronger upward trend from 1960 to 2005 but a clear downward trend from 2005 to 2018. India shows a significant increase in temporal trend. Our findings demonstrate apparent spatial differences in the local temporal trends across countries worldwide.

The global (Table 2) effects of air pollution and urbanization are investigated based on Bayesian regression modeling. The estimate of parameter φ suggested that the spatially structured component explained 7.7% of the variance in the overall spatial effects of GGLE (Table 2). However, the smaller DIC and WAIC scores derived from models with spatial random effects (Table 2) indicated that the spatial structure considered in the Bayesian spatiotemporal regression models could improve the data fitness. Several key findings could be observed. First, a significant positive relationship can be seen between pwPM_2.5 and GGLE, and no significant relations are found between pwPM_2.5 and LE in males, and females can be found. Therefore, considering spatial random effects on the global model, the GGLE increased by 0.00034% when a one percent increase in pwPM_2.5. Second, urbanization shows a significant positive relation to GGLE with spatial random effects, and significant positive relations are found between urbanization and LE in males and females as well (Table 2). Third, significant negative relations between poverty and GGLE can be found, while significant positive relations between calories supply, smoking, and GGLE can be identified.

Further, the six WHO regions exhibit geographic disparities in the effects of pwPM_2.5 and urbanization on GGLE and gender-specific life expectancy (Table 3). First, our results indicate that the negative effect of pwPM_2.5 on GGLE can be found in Africa, Eastern Mediterranean, and South-East Asia. A positive effect of pwPM_2.5 on GGLE is significant in America, Europe, and Western Pacific. No significant relations between pwPM_2.5 and LE in males and females are found in the six regions. Second, our results indicate that the negative effect of urbanization on GGLE can be found in Africa and Europe. In contrast, the positive effect of urbanization on GGLE is significant in America and Eastern Mediterranean. Significant positive relations between urbanization and LE in males and females are found in the six regions.