The Integrated Public Use Microdata Series (IPUMS) based on the Demographic and Health Surveys (DHS) was the source of population data. DHS is a nationally representative survey conducted in developing countries. Since the data is publicly available and anonymous, no further ethical approval is necessary. More details about the DHS data (questionnaire design, sampling method, survey procedures, quality control, etc.) are recorded on the official website of DHS (https://www.dhsprogram.com/ ). The IPUMS DHS is created to facilitate the analysis of DHS data across time and space and has been utilized by over 100,000 researchers worldwide.
Until now, the IPUMS DHS contains standardized population and health data of 9 Asian countries and 32 African countries since the 1980s. In this study, we selected all DHS surveys during 1990–2020 in South Asia (3 countries) and sub-Saharan Africa (28 countries). Information on birthweight in the five years before the survey was obtained through questionnaires, which could have been derived from health cards (45.85%) or from the mother's recollections (52.64%). Supplementary Table 6 shows the countries included and the corresponding survey years. We obtained data from 1,839,959 women for the ages of 15–49 who had up to six recent delivery records, including information about their newborns such as infant sex, birth year, birth month, and birth weight. We also gathered demographic data on the mothers, such as age, education, body mass index (BMI), smoking status (yes or no), place of residence (urban or rural), family wealth index, and GPS information about the survey unit (village or residential cluster). In DHS, all geographical coordinates of sample units were processed and moved 2–10 kilometers in random distances and directions to ensure that individuals and families could not be identified.
The subjects without any delivery records, stillbirths, multiple births, births with abnormal birth weight (< 500g or > 6000g), and those with missing GPS or exposure information were excluded. We further restricted the births before January 2015 to match the availability of the Detection and Attribution Model Intercomparison Project (DAMIP) data, which would be used for attribution estimation of extreme temperature exposure during pregnancy on fetal growth in this study. A total of 451,252 infants from 369,018 mothers were included in the final analyses. See Supplementary Fig. 10 for the screening flow of studied subjects.
Data analysis
The association between extreme temperature and fetal growth, as well as the fetal growth burden attributable to climate change, were estimated by four stages.
Stage 1: Establishment of the exposure-response relationship by countries and regions. Birth weight and LBW (birth weight less than 2500 grams) were used as the main outcomes that indicate fetal growth in this study. To depict the exposure-response relationship between maternal pregnancy temperature exposure and birth weight/LBW for each country, we applied a generalized linear model. The formula is as follows:
$$g\left(Y\right)={\beta }_{0}+\text{n}\text{s}(temperature,2)+{\beta }_{1}*education+{\beta }_{2}*residence+{\beta }_{3}*wealth+{\beta }_{4}*\text{B}\text{M}\text{I}+{\beta }_{5}*sex+ns\left(age,3\right)+ns\left(month,4\right)+ns(humidity,3)$$
1
where \(g\)() is a link function, with ‘gaussian’ and ‘logit’ used for the outcome (Y), birth weight, and LBW, respectively. To test the potential non-linearity for the relationships10,11, pregnancy temperature was included in the model as a natural spline function with degrees of freedom (DF) of 2. A series of confounding factors including maternal education, residence (urban or rural), family wealth index, maternal BMI, and infant sex were added to the models for adjustment. At the same time, the maternal age (DF = 3), the month of delivery (DF = 4), and the relative humidity during pregnancy (DF = 3) were controlled by using the spline function.
Considering the small sample size of some countries and the similar exposure-response relationship pattern, we divided 31 countries into five geographical areas, Southern Asia, and four areas in sub-Saharan Africa (Eastern Africa, Western Africa, Central Africa, and Southern Africa). Then region-specific models were fitted for the exposure-response relationship between temperature and birth weight/LBW in specific regions. Here the country of the mother was additionally included for adjustment.
$$g\left(Y\right)={\beta }_{0}+\text{n}\text{s}\left(temperature,2\right)+{\beta }_{1}*education+{\beta }_{2}*residence+{\beta }_{3}*wealth+{\beta }_{4}*\text{B}\text{M}\text{I}+{\beta }_{5}*sex+ns\left(age,3\right)+ns\left(month,4\right)+ns\left(humidity,3\right)+{\beta }_{6}*\text{c}\text{o}\text{u}\text{n}\text{t}\text{r}\text{y}$$
2
Stage 2: Examination of the effects of extreme temperature exposure on fetal growth by countries and regions. For each country or region, the pregnancy temperature exposure was categorized into 10 layers by every 10-percentile interval, where < the 10th was defined as extreme cold, while > the 90th was extreme hot. From the established exposure-response relationships in stage 1, we obtained the finest temperature (i.e. the temperature corresponding to the lowest risk of reduced birth weight or LBW) for each country or region and the temperature layer where the finest temperature was located.
Then the generalized estimation equation (GEE) with an exchangeable working correlation was applied to estimate the country-specific association between extreme temperature and birth weight or LBW. The formula of GEE was similar to Model (1), except that the natural spline function of temperature in Model (1) was replaced with temperature categories and DHS geocoded residential cluster was added as the random effect. It should be noted that three types of exposure-response relationships in the above analysis were identified: extreme heat risk type, extreme cold risk type, and extreme heat and cold risk type. For the countries belonging to the extreme heat risk type, i.e. the risk of reduced birth weight or LBW was dominated by high-temperature exposure, we estimated only the effects of extreme heat (> the 90th ) by using the lowest temperature layer (< the 10th ) as the reference. On the contrary, only the effects of extreme cold (< the 10th ) were estimated for the countries of low-temperature risk type by referring to the highest temperature layer (> the 90th ). The effects of extreme heat and extreme cold were estimated for the countries where both high and low-temperature exposures exert considerable impacts on reduced birth weight or LBW. In this case, the temperature layer containing the finest temperature was treated as the reference. The effects were expressed with regression coefficients or ORs and their 95% CIs. We also used the same GEE to estimate the region-specific associations between extreme temperature and birth weight or LBW with additional adjustment of the country of mother.
Stage 3: Estimation of extreme temperature-related LBWs attributed to anthropogenic climate change within studied samples. Considering the high LBW rate in LMICs, we first need to calculate the relative risk (RR) using OR41.
$${RR}_{ie}={OR}_{ke}/[\left(1-{P}_{0i}\right)+{P}_{0i}*{OR}_{ke}]$$
3
The OR represents the estimated association between the exposure to extreme temperatures and LBW for each region in stage 2. P0 represents the LBW rate. The subscript character i represents a specific country, e represents extreme heat or cold, and k represents the region that the country belongs to.
In the studied population, we assessed the pregnancy exposure for each birth using both hist-nat and historical temperature simulated by six models of DAMIP. The average temperature of nine months before the delivery was calculated as the temperature exposure during pregnancy under counterfactual and factual scenarios, respectively. In each country, according to the pregnancy exposure, newborns were divided into three groups: extreme cold, extreme heat, and non-extreme temperature exposed groups using country-specific thresholds for extreme temperatures established during Stage 2. Then, the PARP, which indicates the proportion of LBW resulting from exposure to extreme temperatures, was calculated using the following formula 42:
$${PARP}_{ie}={P}_{ie}\text{*}({RR}_{ie}-1)/[{P}_{ie}\text{*}\left({RR}_{ie}-1\right)+1]$$
4
in which Pie represents the percentage of the newborns in country i exposed to e level. It should be noted again, for countries with extreme heat and cold risk types, the PARP attributed to extreme temperature should be the sum of PARPcold and PARPheat. For the countries of the extreme heat risk type, only extreme heat-related PARP was calculated because there was no RRcold available, and only extreme cold-related PARP was calculated for the countries belonging to the extreme cold risk type.
Here PARPs under both counterfactual and factual scenarios were calculated, and the difference of PARP under the factual scenario minus that under the counterfactual scenario was used to represent the extreme temperature-related LBWs attributable to anthropogenic climate change in various countries and regions. To intuitively show the extent to which human-induced climate change affects LBW related to extreme temperatures exposure, we further calculated the attributable proportion (AP).
$${PARP}_{anthropogenic}={PARP}_{\text{f}\text{a}\text{c}\text{t}\text{u}\text{a}\text{l}}-{PARP}_{counterfactual}$$
5
$${AP}_{anthropogenic}=\frac{{|PARP}_{anthropogenic}|}{\text{max}\left({PARP}_{\text{f}\text{a}\text{c}\text{t}\text{u}\text{a}\text{l}}, { PARP}_{counterfactual}\right)}*100\%$$
6
We calculated the PARP and AP by each country and region using the simulated temperature data from six DAMIP models separately. Six sets of PARP and AP were output and their mean was used as the main results. We also estimated the 95% CI for PARPs and APs. Specifically, we assumed that the regression coefficients for the relationship between extreme heat/cold and LBW follow a multivariate normal distribution, and used Monte Carlo to generate 1000 coefficient samples for each country and region. All these coefficients were converted into ORs for the calculaton of PARPs and APs. So that we were able to obtain a total of 6,000 sets of PARPs and APs, respectively. Then, the 2.5th to 97.5th quantiles of PARPs and APs were used as their corresponding 95% CIs24,25. This method enabled us to obtain the practical CIs after fully considering the uncertainty of the exposure-response function and the variations between different climate simulation models.
Stage 4: Estimation of extreme temperature-related LBWs attributable to anthropogenic climate change for the whole studied countries. We divided 31 studied countries into grids of 0.5 ° × 0.5 °, and attribution analysis was conducted on each grid point to capture the spatial variation. To achieve this, we collected 2000–2014 global population data (30 arc-second resolution) from the LandScan program (https://landscan.ornl.gov/). LandScan Global generates the finest resolution global population data representing the environment (average 24 hours) population since 2000, by combining geospatial science, remote sensing technology, and machine learning algorithms. We used the upscaling technology in ArcGIS 10.5 (Environmental Systems Research Institute, Redlands, California, USA) to make its resolution to be 0.5 ° × 0.5 °. In addition, the national crude birth rates (per 1,000 people) for 2000–2014 were obtained from the World Bank (https://data.worldbank.org/). We assume that the birth rate in all grids within a country is the same.
Since the population data were available after 2000, we quantified the extreme temperature-related LBWs attributable to anthropogenic climate change during the period 2000–2014. First, we needed to calculate the monthly number of baseline newborns per year for each mesh with the following formula:
$${baseBirth}_{jm}={{Y}_{j}*BirthRate}_{i}*{P}_{im}$$
7
where j is the grid and Yj represents the 24-hour averaged population of each grid in country i, BirthRatei represents the birth rate of country i in a specific year, and Pim represents the proportion of births in month m of country i, which was assumed to be the same as studied samples.
The temperature exposure during pregnancy for newborns in each month (January to December) and each grid were assessed using the same method as stage 3. By applying the thresholds for extreme heat and cold, we identified the type of pregnancy exposure for births in each month, and the proportion of the births with pregnancy exposure to extreme heat and cold during the year. Using formula (4), we calculated the extreme heat-related PARP and extreme cold-related PARP for each grid in the current year. Then the number of LBW attributable to anthropogenic climate change (aLBWs) could be estimated using the following formula.
$${\text{a}\text{L}\text{B}\text{W}\text{s}}_{j}={baseBirth}_{j}*{LBWrate}_{i}*{PARP}_{j}$$
8
Extreme temperature-related LBWs attributable to observed and anthropogenic climate change, as well as their 95% CIs were calculated in the same way as we used in stage 3. To further show the temporal trend, we summarized the LBWs attributable to anthropogenic climate change by two periods, 2000–2007, and 2008–2014 year.
Exploring the modulator for spatial variation in both association and attribution. From the perspective of physical geography, we speculated that elevation is an important driver for the spatial variation of the estimated association and attribution. Based on the elevation of included areas, we divided the studied samples into high-altitude (> 1000m) and middle- and low-altitude groups ( < = 1000m). The above-described analyses, including the establishment of the exposure-response relationship, and estimation of the association and attribution, were conducted in subgroups.
Sensitivity analyses. We performed several sensitivity analyses to evaluate the strength and reliability of the association between extreme temperature and birth weight/LBW. First, we performed a series of models, including an unadjusted model, a model only adjusting for relative humidity during pregnancy, a model with more covariates (relative humidity, maternal education, BMI, residence, family wealth index, and infant sex), and also a final model additionally included the delivery month and country, to test if the association between extreme temperature and birth weight/LBW was sensitive to covariate adjustment. Second, since the population data were collected retrospectively in the DHS survey, the estimation may be affected by the potential recall bias. To explore this, we restricted the samples within 3, 4, and 5 years before the survey to re-estimate the association. Meanwhile, we conducted the analyses on study subjects whose birth weight was derived from written records rather than the mother's self-report. Third, we used more cutoff values to define extreme temperatures, i.e., using the 5th or 25th percentiles to define extreme cold, and the 95th or 75th percentiles for extreme heat, respectively, to test if the effects of extreme temperature on fetal growth vary with definitions.