Long-term ozone exposure associated cause-specific mortality risks with adjusted metrics by cohort studies: A systematic review and meta-analysis

2.1 Search strategy

We searched 3 research databases (MEDLINE, Embase, and Web of Science) from 1 September, 2015 till 1 February, 2022 to finish our systematic review and meta-analysis, updated from 2 previous reviews on long-term O₃ exposure-associated mortality.^{15, 16} Search terms also referred to these 2 previous systematic reviews with modifications to enhance the inclusion of potential relevant studies, as we combined the keywords relevant to the cause-specific mortalities (i.e. “mortality”, “death”, “premature death”, “all-cause”, “non-accidental”, “cardiopulmonary”, “respiratory”, “chronic obstructive pulmonary disease”, “pneumonia”, cardiovascular, “lung cancer”, “cerebrovascular”, “stroke”, “ischaemic heart disease”, “congestive heart failure”), the pollutant of research interest (i.e. “ozone”), and qualified epidemiological study types (i.e. “long-term”, “cohort study”, “prospective”, “retrospective”, “longitudinal study”). The detailed search strategies were listed in Table S1. Health outcomes considered in the systematic review were: mortality from (1) all causes (AC, ICD9: 001-799, ICD10: A00-R99); (2) all respiratory diseases (RESP, ICD9: 460-519, ICD10: J00-J98); (3) chronic obstructive pulmonary diseases and allied conditions (COPD, ICD9: 490-496, ICD10: J19-J46); (4) all cardiovascular diseases (CVD, ICD9: 390-459, ICD10:I00-I98); (5) all cerebrovascular diseases (CEVD, ICD9: 430-438, ICD10: I60-I69); (6) ischaemic heart disease (IHD, ICD9: 410-414, ICD10: I20-I25); (7) congestive heart failure (CHF, ICD9: 428, ICD10: I50); (8) ischaemic stroke (ICD9: 434, ICD10: I61-I64); (9) pneumonia (ICD9: 480-487, ICD10: J12-J18); and (10) lung cancer (LC, ICD9: 162, ICD10: C33-C34).

2.2 Study eligibility criteria

As an updated systematic review, literatures identified in the previous 2 reviews underwent examination together with the newly retrieved ones. Studies were included during screening following the criteria as: (1) the epidemiological research should be conducted based on cohorts; (2) the exposure should include O₃ as an individual risk factor; (3) the health outcomes should be all-cause or cause-specific deaths at individual level; (4) studies provided hazard ratio (HR), risk ratio (RR) or odds ratio (OR) and their 95% confidence intervals (CIs) clearly and reported by every increase unit (e.g. 10-ppbV) of exposure concentrations, assuming linear risk relationship with adjusting key confounders; (6) the study should be published as an original research article in scholarly peer-reviewed journals in English. For articles from the same cohort, only one single study covering the widest populations and the longest follow-up period was reserved for meta-analysis, unless the subgroups of participants and study follow-up periods are clearly stated to be of mild overlapping; We followed the Preferred Reporting Items for Systematic Review and Meta-Analyses (PRISMA) guidelines to process the included studies on ambient O₃ exposure induced mortality.

2.3 Study selection and scrutinisation

All searched literatures were archived in Clarivate™ Analytics Endnote X9.3.1 reference manager software. Two literature review investigators (HZS and CL) conducted title and abstract pre-screening independently for all web-searched records and reviewed the full text for the pre-screened studies. Disagreements were resolved by discussions with a third reviewer (PY).

2.4 Data extraction

Details from each screened-out literatures were extracted and labelled for the purpose of meta-analysis, including (1) the authors with publication year as study labels of reference; (2) basic descriptive information of the study cohort embracing the cohort name, country, follow-up periods, numbers of cases and total participants, population genders and ethnics, exposure metrics, health outcomes, and major confounders; (3) the risk association effects preferably quantified in HR (and also RR/OR as substitute choices) per unit incremental exposure with 95% confidence interval (CI).

2.5 Study quality assessment

All screened-out studies underwent quality evaluation using the Quality Assessment Tool of Observational Cohort and Cross-Sectional Studies developed by National Institute of Health (NIH) (https://www.nhlbi.nih.gov/health-topics/study-quality-assessment-tools), aiming to ensure the studies considered for meta-analysis are adequately reliable. The assessments were cross-validated by two authors (HZS and CL) independently, with the third author (PY) supervising any disagreements. Table S2 listed 14 assessment items assigned with 1 score for each, and the tallied scores were translated into a literature-specific rating of quality. Studies scoring full-mark 14 were categorised to be “Good”, while 10-13 as “Fair” and <10 as “Poor”.

Besides applying the quality assessment tool to determine which reviewed studies should be included for meta-analysis, checking the epidemiological evidence quality from the included literatures for each cause of mortality was finished through the Grading of Recommendations Assessment, Development, and Evaluation (GRADE) system^{21, 22} to yield rating bands ranging across “high”, “moderate”, “low”, and “very low”. This grading system by default rated “high-quality” for cohort studies as the starting point of evaluation, and the rate would be downgraded by five limitations as the existence of (1) risk of bias examined by the Quality Assessment Tool (Table S2), (2) imprecision (i.e. studies did not report the central risk estimations with confidence intervals), (3) inconsistency (i.e. the directions of the estimated risks were controversial across studies), (4) indirectness (i.e. studies did not include the desired population, exposure, or health outcomes), and (5) publication bias (i.e. researchers tended to publish studies with positive results); and upgraded by three strengths as reporting (1) exposure-response trend, (2) residual confounding (i.e. adjusting the confounders highlighted the risks), and (3) strong associations. Publication biases were graphically presented by funnel plots,²³ and statistically tested by trim-and-fill method.²⁴ The review was registered in PROSPERO (CRD42021270637).

2.6 Exposure adjustment

2.7 Meta-analysis

We collectively named relative risks (RR) for HR/RR/OR throughout our meta-analysis. All literature-reported RRs were converted into adjusted incremental risk ratios with a 10-ppbV O₃-exposure increase by target metric (i.e. 6mDMA8 in this study), following eq 6 as shown below: where ln is the natural logarithm, RR_Originalis the originally reported risk estimates scaled into 10-ppbV incremental exposure, and k_{Original→Adjusted} is the conversion factor for metric unification. Multi-study pooled risks with 95% confidence interval (CI) were calculated from the adjusted RRs by Hunter-Schmidt random effect meta-regression estimator to correct the potential errors and biases caused from the diversity of study population and methodologies.³⁵

We applied the Higgins I ²to quantify the heterogeneity across studies. The Higgins statistics I ² is defined as where Q is the Cochran’s non-parametric heterogeneity statistic assessing whether there are any cross-study differences in risks based on x ² distribution and df is the corresponding degrees of freedom.³⁶ Low I ² values indicate no important heterogeneity observed and high I ² values, especially >75%, indicate considerable heterogeneity.

Subgroup analyses were conducted by grouping the selected studies upon the gender, regions, O₃ exposure metrics, and methodological reliability of individual exposure assignment; together with the adjustment of ethnicity, body mass index (BMI), smoking history, lifestyle features, and exposure to PM_2.5 and NO₂. Subgroups should contain at least 3 studies. Leave-one-out sensitivity analyses were also accomplished to test the robustness of synthesised overall risks by meta-analysis. All meta-analyses were performed in R 4.1.1 with packages meta, metafor, and metainf.

The most widely recognised approach to construct the integrated exposure-response (IER)³⁷ relationships required sufficient epidemiological studies to comprehensively sample the population exposure levels. However, studies on long-term O₃ exposure health effects were relatively limited, under which circumstance we made methodological modifications to make better use of the variabilities in exposure levels by statistically imputing the exposure distributions for each study from the provided statistics (e.g. mean, standard deviation, and percentiles) for curve fitting as elaborated in Supplementary Text S1. Supplementary Text S2 described the detailed procedures of exposure distribution imputations with a demonstration provided in S3, through which high uncertainties were still observed in the fitted IER curves due to insufficient epidemiological evidences.

3.1 Study characters

From the 3 databases during September 2015 till February 2022, a total of 339 studies (77 from MEDLINE, 102 from Embase, and 160 from Web of Science) were searched; and together with 34 additional literatures added manually from the 2 previous systematic reviews,^{15, 16} 373 studies underwent duplication censoring, deleting 101 duplicated studies. After detailed scrutinisation for 272 studies, a total of 25 studies concerning long-term O₃ exposure and multi-cause mortalities were finally enrolled for quality evaluation, meta-analysis and further discussions (Figure 1).^38-62 Table 1 summarised the basic information of the 25 included studies sorted by the publication year and surname of the first author.

• Download figure
• Open in new tab

Figure 1 Schematic flowchart of study assessment and selection processes for literature review and meta-analysis.

3.2 Metrics and exposure assignments

Our updated systematic review stressed more on the exposure metrics and methodologies to obtain O₃ exposure, as summarised in Table 2. Abbey et al. (1999),³⁸ Jerrett et al. (2013)⁴⁶ and Lipsett et al. (2011)⁴³ did not state the metric they used clearly, but based on comparisons between the reported surface O₃ concentrations and TOAR observational archives, we reasonably assumed ADA24 for the first study, and ADMA8 for the rest two. Details of the metric matching were given in Supplementary Text (S4). Lipfert et al. (2006)³⁹ used the highest 95^th percentile by hourly resolved O₃ concentrations as the peak exposure metric, which was only used in this one singular study, and hence approximated to DMA1. Krewski et al. (2009)⁴¹ and Smith et al. (2009)⁴² were both studies on ACS CPS II, and thus the same exposure assignment methodologies and metrics were assumed as Jerrett et al. (2009).⁴⁰ Likewise, Cakmak et al. (2018)⁵³ and Weichenthal et al. (2017)⁵² were assumed to inherit Crouse et al. (2015)⁴⁸ as all these 3 studies were on CANCHEC. Warm season was defined as 6 months from April to September in terms of the northern hemisphere by default, but we made no exceptions to 3 studies as Zanobetti et al. (2011)⁴⁴ using May to September, and Crouse et al. (2015)⁴⁸ and Paul et al. (2020)⁵⁷ using May to October, due to limited number of studies searched.

Across all included studies, multiple methods were applied to obtain gap-free surface O₃ concentrations for individual-level exposure assignment. The most basic way was the nearest neighbour matching between participant residential locations and in situ observation sites, which were more frequently used in earlier studies.^{39, 40} A comparatively more complicated way was statistical spatial interpolation, by inverse distance weighting⁴⁶ or ordinary kriging⁴¹. Full spatial coverage products, such as satellite-based remote-sensing⁵¹ and chemistry transport models⁵⁶, were used in some studies by supervised-learning-based data fusion techniques including but not limited to universal kriging embedded land use regression,⁴⁷ Bayesian hierarchical model⁵⁰ and ensemble learning⁵¹ to enhance the spatial extrapolation accuracy, which were evaluated to be of higher credibility than the basic ones described previously. All basic interpolation methods using merely the observations were rated as “Low”; applying chemical transport model simulations without calibration from the observations as “Moderate”; linearly coupling the observations with simulations as “Good”; and multi-source data assimilation by means of more sophisticated approaches as “High”. To sum up, 8 studies were rated “High”, 5 were “Good”, 2 were “Moderate”, and 10 were “Low”. Methodological progresses with time were evident as manifested in Table 2, prefiguring an explosion of population-based environmental health studies in the age of big data.

Based on the TOAR and CNEMC in situ observations, the cross-metric non-intercept linear conversion factors were estimated with regression accuracies given in Figure 2. Synthesised from the recent relevant studies, the 6mDMA8 metric was more recommended to highlight the peak exposure; and therefore, we chose to convert the originally reported RRs uniformly into the 6mDMA8 scale as standard. The O₃ exposure levels by the original and unified metric were listed in Supplementary Text S1. Demonstrations for the conversion interpretation and procedures were presented in Supplementary Text S5, respectively.

• Download figure
• Open in new tab

Figure 2 Cross-metric linear relationships and conversion accuracies.

The cross-metric linear relationships were scaled by Pearson’s correlation coefficients. The cross-metric conversion factors with 95% confidence intervals (95% CI) were estimated by non-intercept linear regression models, accompanied with fitting accuracies quantified by coefficient of determination (R²) and root mean square error (RMSE) in ppbV. The conversion factors were defined as multiples from the original metric by column into the target harmonised metric by row, e.g. ADMA8 = 1.671 ADA24, R² = 0.9736, RMSE = 7.78 ppbV. Note that by non-intercept linear regression, the values of R² should no longer be equal to the squared Pearson’s linear correlation coefficients. As the cross-metric conversion coefficients were estimated statistically, indirect conversions were not recommended, since regression noises restricted the validity of equation k_A→B = k_A→C. k_C→B.

3.3 Meta-analysis results

We conducted meta-analyses for long-term O₃ exposure-associated into 10 categories of mortalities as (1) all causes (AC), (2) all respiratory diseases (RESP), (3) chronic obstructive pulmonary diseases and allied conditions (COPD), (4) all cardiovascular diseases (CVD), (5) all cerebrovascular diseases (CEVD), (6) ischaemic heart disease (IHD), (7) congestive heart failure (CHF), and (8) lung cancer (LC), with the exposure metrics adjusted into 6mDMA8.

3.3.1 All-cause mortality

A total of 23 studies were included into O₃ exposure-associated all-cause mortality meta-analysis, pooling the overall risk into RR = 1.014 (95% CI: 1.009–1.019, I²: 97.8%) with every 10-ppbV incremental exposure by 6mDMA8 as presented in Figure 3. Sub-group meta-analysis by originally reported metrics concluded the significances of risks vary across metrics, as high-concentration highlighted metrics like 6mDMA8 were of the highest positive risk (RR = 1.022, 95% CI: 1.014–1.031) while the smoothed metric ADA24 reported negative association (RR = 0.980, 95% CI: 0.960–1.001), as shown in Figure S1. Grouped by study regions, significant discrepancies of the risk pattern was found (Figure S2), as the studies in North America revealed positive associations as RR = 1.019 (95% CI: 1.014–1.024), while researches on European populations showed reversed risks as RR = 0.910 (95% CI: 0.827–1.001), though not significant. The cross-region divergences did not necessarily indicate differences in population vulnerability, as (1) less and younger study population, (2) shorter follow-up durations, and (3) use of smoothed exposure metrics for studies in Europe could all potentially obscure the potential risk associations. Subgroup analysis manifested that high inter-study heterogeneities originated from metric inconsistency, methodological reliability of individual exposure assignment, and population variabilities, as encapsulated in Table S3. The funnel plot was visually symmetrical (Figure S3), and studies reporting risks below the pooled value were even slightly more, indicating no detected severe potential publication biases.

• Download figure
• Open in new tab

Figure 3 Pooled estimates of all-cause mortality risk associated with every 10-ppbV incremental O₃ exposure by 6mDMA8 metric.

Size of the shaded squares in the forest plot represents the weight of each study estimated by random-effect model.

No significant inter-gender differences were observed based on the limited studies reporting gender-specific risk association strengths. Further subgroup analyses were unfeasible due to the lack of reporting in the literature. Alternatively, grouped RRs were estimated based on whether the original researches had adjusted the confounding effects from ethnicity, body mass index (BMI), smoking history, lifestyle features, exposure levels of PM_2.5 and NO₂, and no inter-group divergences were observed (Table S3).

3.3.2 Respiratory mortality

Meta-analysis for O₃ exposure-associated all respiratory mortality includes 16 studies, pooling which gave the overall RR = 1.025 (95% CI: 1.010–1.040, I²: 83.9%) for every 10-ppbV incremental O₃ exposure by 6mDMA8 (Figure 4). Based on sub-group meta-analysis for different metrics (Figure S4), peak metrics showed more significant increasing risks than ADA24, where most of the heterogeneities were from (I² = 87.0%). Cross-metric divergences were generally higher than intra-metric discrepancies. Studies on North America populations showed better homogeneity in positive risks (RR = 1.029, 95% CI: 1.011–1.047, I² = 71.1%, Figure S5) than the European cohorts, pooling from which the overall risks were congruously insignificant (RR = 0.941, 95% CI: 0.856– 1.036, I² = 91.2%). For O₃-COPD mortality association, the pooled RR was 1.056 (95% CI: 1.029–1.084, I² = 94.5%) for 10-ppbV incremental O₃ exposure by 6mDMA8 from 7 studies. No apparent positive publication biases were detected for both respiratory and COPD mortalities from the funnel plot (Figure S3).

• Download figure
• Open in new tab

Figure 4 Pooled estimates of respiratory diseases and COPD mortality risks associated with every 10-ppbV incremental O₃ exposure by 6mDMA8 metric.

3.3.3 Cardiovascular mortality

A total of 15 studies were included to pool the overall O₃ exposure-induced CVD mortality risks as RR = 1.019 (95% CI: 1.004–1.035, I² = 97.7%) for each 10-ppbV additional O₃ exposure by 6mDMA8 (Figure 5). Given the fact that the lower bound of uncertainty interval was so close to the null hypothesis (i.e. RR = 1), the positive risk association found in this review could be controversial, and thus would require more studies to support or refute the finding. Heterogeneities (I² > 79.2%) were observed through all 3 metric-grouped studies as presented in Figure S6. Positive risk associations were found on 10 North American cohorts (RR = 1.036, 95% CI: 1.019–1.053) while oppositely for 5 European cohorts (RR = 0.934, 95% CI: 0.865–1.008), as shown in Figure S7. There was no need to be concerned with the publication bias, and no more inter-group divergences were spotted except for grouping by exposure assignment methodological credibility (Table S3). The pooled risk for the congestive heart failure-induced mortality from 4 studies was RR = 1.074 (95% CI: 1.054–1.093, I² = 85.8%) with every 10-ppbV incremental O₃ exposure by 6mDMA8.

• Download figure
• Open in new tab

Figure 5 Pooled estimates of cardiovascular diseases and congestive heart failure mortality risk associated with every 10-ppbV incremental O₃ exposure by 6mDMA8 metric.

3.3.4 Other mortality causes

The other cause-specific mortality risks attributable to long-term O₃ exposure were not statistically significant (Figure 6), as IHD mortality risk pooled from 10 studies was RR = 1.012 (95% CI: 0.987–1.039, I² = 98.7%), CEVD mortality risk pooled from 6 studies was RR = 0.993 (95% CI: 0.979–1.008, I² = 80.6%), and LC mortality risk pooled from 13 studies was RR = 0.966 (95% CI: 0.926–1.007, I² = 84.2%). For all 8 studied mortality causes, we also provided pooled risks by 3 more widely used metrics (i.e. 6mDA24, ADMA8, and ADA24) besides 6mDMA8, as listed in Table 3 for reference.

View this table:

• View inline
• View popup
• Download powerpoint

Table 3 Pooled RRs for long-term 10-ppbV incremental O₃-exposure attributable multi-cause mortalities by 4 most widely used metrics and crude risks without metric harmonisation.

• Download figure
• Open in new tab

Figure 6 Pooled estimates of ischaemic heart disease, cerebrovascular diseases, lung cancer, ischaemic stroke, and pneumonia mortality risks associated with every 10-ppbV incremental O₃ exposure by 6mDMA8 metric.

3.4 Study assessment

All 25 studies included into our final meta-analysis were rated to be above “Fair” (14 as “Fair” and 11 as “Good”) by the Quality Assessment Tool for Observational Cohort Studies, as listed in Table S4. All studies well met 10 out of 14 assessment items, while 9 studies did not sufficiently clarify the participant exclusion criteria; 2 re-analysis study reports did not clearly state the O₃ exposures;^{41, 42} 2 studies were of relatively insufficient follow-up durations (e.g. less than 5 years) to observe the outcomes resulting from long-term exposure;^{45, 63} and 10 studies were of methodological deficiencies in individual exposure assignment,^{38-46, 55} most of which were conducted before 2013 when data assimilation techniques were not maturely developed to fuse observations with other full spatial coverage products such as satellite-based remote sensing and atmospheric mechanistic simulations. The satisfactory assessment results overall indicated inconspicuous risks of bias, laying the reliable foundation for meta-analyses.

Tables S5 displayed GRADE epidemiological evidence assessment results for each mortality cause from all involved corresponding studies. In brief, the overall judgements for all-cause, respiratory, cardiovascular, ischaemic heart disease, congestive heart failure, and lung cancer mortality risks were “High”, while the rating for the rest 2 cause-specific mortality risks (COPD and cerebrovascular diseases) were both “Moderate”.

Inconsistence of the risk directions (i.e. positive or negative associations) was the most common reason for downgrading, except for the CHF-induced mortality. There were 6 studies having reported the O₃-mortality exposure-response trends to support the additional risks, as an assessment upgrading item for the pooled RRs of all-cause, respiratory and cardiovascular mortality. Cakmak et al. (2018) spotted higher RRs after adjusting the confounder compared to the crude values,⁵³ which gave prominence to the positive risk associations and thus correspondingly upgraded the rating for all-cause, ischaemic heart disease, and lung cancer mortalities. No substantial positive publication biases were found based on the collected evidences.

3.5 Sensitivity analysis

Leave-one-out sensitivity analyses showed stable risk estimates as summarised in Table S6, except for the lung cancer mortality risks after eliminating Kazemiparkouhi et al. (2019), the only study reporting positive risk association,⁵⁵ while the rest 11 studies concluded insignificant risks or even protective effects. Since the metric harmonisation in our study was an innovative attempt, we provided both metric-adjusted and unadjusted crude results for reference as presented in Table 3. The crude results were pooled from the originally reported relative risk values only unified into per 10-ppbV incremental exposure, without being transformed into any metrics for congruity. Along with the meta-analyses on all qualified studies, the relative risks were also pooled by keeping only one latest study with the largest population for each separate cohort, as summarised in Table S7. Under this circumstance, the pooled unit incremental mortality risks with every 10-ppbV incremental O₃ exposure by 6mMDA8 metric were RR = 1.008 (95% CI: 1.006– 1.009, I² = 82.6%) for all causes, RR = 1.034 (95% CI: 1.017–1.050, I² = 81.7%) for all respiratory diseases, RR = 1.060 (95% CI: 1.040–1.080, I² = 90.2%) for COPD, RR = 1.032 (95% CI: 1.010–1.055, I² = 98.2%) for all cardiovascular diseases, RR = 1.008 (95% CI: 0.973–1.045, I² = 99.2%) for ischaemic heart disease, and RR = 0.966 (95% CI: 0.931–1.002, I² = 83.8%) for lung cancer. Studies for mortality risks of cerebrovascular diseases and congestive heart failure were respectively conducted on different cohorts, and hence such supplementary analysis was unnecessary.

4.1 Improvements as an updated review

This work improves on 2 previous high-quality reviews^{15, 16} by covering up-to-date peer-reviewed studies, and expanding the O₃-exposure associated causes of mortality into wider range of categories. It is the first systematic review of the association between long-term O₃ exposure and cause-specific mortality highlighting the issue of inconsistent use of exposure metrics to our best knowledge. Since tropospheric O₃ is a photochemical pollutant which largely depends on solar radiation, the surface O₃ concentrations can vary drastically between day and night, as well as warmer and cooler seasons. We pointed out that a 10-ppbV increase in annual daily 24-hour average concentration (ADA24) is more constrained in magnitude than a 10-ppbV increase in warm-season daily 8-hour maximum average concentration (6mDMA8) owing to the wider variability in the range of the latter metric. Taking the observations by TOAR and CNEMC in situ monitoring networks during 1990-2019 as an example, the surface O₃ concentrations were 27.6 ± 6.1 (IQR: 24.1–31.0) ppbV by ADA24, while correspondingly 53.1 ± 10.6 (IQR: 47.7–61.4) ppbV by 6mDMA8, which indicated a 10-ppbV change fell below the IQR by the 6mDMA8, but could exceed the IQR using the ADA24 metric. This was why we believe adjusting the exposure metrics was necessary for O₃ exposure-attributable health risk meta-analysis.

We also put forward a feasible approach to mutually convert the O₃ exposure concentrations and corresponding risk strengths in various metrics by non-intercept linear projections following the operational suggestions from EPA,³² but update the linear conversion factors based on global in situ surface O₃ observations during 1990-2019. The methodological innovation took advantage of multi-dimensional information from the original studies, which could inspire further observation collections and researches for corroborations and improvements.

4.2 Metric relevant issues

Although linear coefficients were applied onto the cross-metric conversions, irreducible noises still existed given the high root mean squared errors (RMSE) as shown in Figure 2, which exposed the limitation of risk strength adjustment into the same exposure metric by simple linear conversion, as the actual cross-metric relationships could be way more complicated. However, there was no other way but using the linear conversion coefficients as surrogates to unify the RRs by different metric reported in original studies, and thus to avoid uncertainties brought by the conversion of metrics, using a promissory consistent exposure metric or estimating the unit excess RRs in multiple metrics would be highly advocated in future long-term O₃-exposure epidemiology studies.

Such linear conversion of risk associations could be validated by Kazemiparkouhi et al. (2020),⁵⁵ where multiple metrics were applied to estimate the mortality risks. For COPD mortality, the RR was 1.072 (95% CI: 1.067–1.077) by 6mDMA1 for every 10-ppbV additional exposure; and after converting into 6mDMA8 metric using the linear coefficient 0.831 (Figure 2), the estimated RR was 1.087 (95% CI: 1.081–1.093), close to the literature reported 1.084 (95% CI: 1.079–1.089),⁵⁵ which justified our linear conversion method. Converting Cross-metric linear conversions would not change the risk association direction, but using different exposure metrics when estimating the O₃-exposure attributable mortality risks could potentially cause discrepancies. For an instance, Kazemiparkouhi et al. (2020) concluded excess hazards of long-term O₃ exposure on all-cause mortality using 6mDMA1 and 6mDMA8 as quantitative metrics, but 6mDA24 led to a specious prevention effect (RR = 0.990, 95% CI: 0.988–0.991), which should be attributed to the existence of a theoretical exposure safety level for O₃ below which no negative health effects should occur. Under this circumstance, lower-level metrics (e.g. ADA24) by averaging the peak O₃ exposures might obscure the effective doses above the threshold, and also reduce the signal-to-noise ratios, so that were of lower credibility in recognising hazardous population exposures than higher-level metrics (e.g. 6mDMA8).

Data mining techniques were able to realise high-accuracy predictions of surface O₃ concentrations, but errors were never avoidable. Carey et al. (2013) used a basic IDW spatial interpolation approach to obtain the surface O₃ concentrations where the R² were 0.24–0.76,⁴⁵ while years later Di et al. (2017) applied an ensemble learning approach, achieving R² = 0.80, RMSE = 2.91 ppbV.⁵¹ Carey et al. (2013) reported the IQR of O₃ exposure concentrations as 3.0 ppbV, which was comparable to the RMSE of Di et al. (2017).⁵¹ Besides, lower R² could be accompanied with higher prediction errors, which might have concealed the highest and lowest quartiles, and led to failures in distinguishing the population-level exposures. This concern had been reflected in our subgroup meta-analysis by exposure metrics, that lower-level metrics were more inclined to report insignificant risks, which also cast sceptics on the reliability of studies covering narrow exposure variabilities. We therefore are in favour of the Lancet suggestions to use peak metrics to quantify the long-term O₃ exposure such as 6mDMA8, and also speak highly of the state-of-the-art data techniques to reduce errors in O₃ concentration prediction, so as to make a distinction between the high- and low-exposure populations.

4.3 Pathogenesis supports

Atkinson et al. (2016) concluded insignificant pooled risks for long-term O₃ exposure associated all-cause and respiratory mortality,¹⁶ which contradicted our results. It should mainly be ascribed to the heterogeneity between the more recent studies and earlier ones. The majority of studies collected in Atkinson et al. (2016) applied primitive statistical methods (i.e. nearest neighbourhood matching, IDW and ordinary kriging interpolation) for individual exposure assignment, which might have weakened the individual-level exposure distinguishment. In addition, some studies using ADA24 as the exposure metric could have also obscured the significance of associations.^{43, 45, 46} In contrast, studies after 2016 more frequently applied advanced numerical simulation models and data assimilation techniques to increase the precision of population exposure assessment; and most of them used 6mDMA8 metric to foreground the high exposures.^{50, 52, 53, 55-57} These recent studies stuck out the significant O₃-mortality associations.

From another aspect, pathological mechanisms had been at least partially ascertained by laboratorial experiments. The inhaled O₃ could constrict the muscles in the airways leading to shortness of breath, and damage the lining with inflammation.⁶⁴ Long-term O₃ exposure could increase the oxidative stress in the cardiovascular system,⁶⁵ and cause progressive thickening of the carotid arteries to restrict cerebral blood supply.⁶⁶ Additionally, strong associations had already been found between short-term O₃ exposure and a variety of cardiopulmonary symptoms as reported by a number of observational epidemiological studies,⁶ which also supported the long-term exposure effects, as it was unreasonable to presume no incremental risks by long-term exposure given the verified significant short-term effects. We thus were inclined to approve of the opinion that long-term O₃ exposure would increase mortality risks in agreement with GBD report.²⁰

To alleviate the population health loss resulting from O₃ exposure, the U.S. EPA appealed for optimisations in real-time accessibility of air quality index, with which residents could be able to avoid unnecessary high pollution exposure (https://www.epa.gov/ground-level-ozone-pollution/health-effects-ozone-pollution). Appropriate diets and supplements including carotenoids, vitamin D and vitamin E were recognised to be preventive against air pollution induced respiratory damages, which was a practical protective measure for the vulnerable.⁶⁷

4.4 Concentration-response relationship

Few studies had examined the concentration-response curves between long-term O₃ exposure and mortality, and thus the threshold exposure level (also known as theoretical minimum risk exposure level, TMREL) below which no adverse health effects would be assumed to occur was still controversial. For all-cause mortality, Di et al. (2017) reported a safe exposure level as 30 ppbV by 6mDA24 metric (approximately as 49.9 ppbV by 6mDMA8),⁵¹ while Shi et al. (2022) suggested a lower level as 40 ppbV by 6mDMA8, both estimated from the Medicare beneficiary cohort.⁵⁸ For respiratory mortality, Jerrett et al. (2009) tested the concentration-response relationships and estimated the threshold level as 60 ppbV by 6mDMA1 (49.9 ppbV by 6mDMA8),⁴⁰ while Lim et al. (2019) failed to identify a significant threshold level.⁵⁶ For cardiovascular mortality, Lim et al. (2019) showed no apparent health hazards below 45 ppbV by 6mDMA8,⁵⁶ and Paul et al. (2020) prescribed a threshold level around 35 ppbV by 6mDMA8 metric for diabetic patients.⁵⁷ These evidence-based threshold exposure levels were all no higher than the current standards, as 70 ppbV for daily maximum 8-hour exposure under NAAQS (The National Ambient Air Quality Standards regulated by the U.S. EPA)⁶⁸ and 50 ppbV by warm-season DMA8 under WHO global air quality guidelines.⁶⁹ However, whether the standard guidelines should be revised to be more strictly would require more further studies.

To synthesise epidemiological evidences, Burnett et al. (2014) developed an integrated exposure-response (IER) function-based curve-fitting method to pool the risk associations from multiple studies.³⁷ We attempted to construct the IER for long-term O₃ exposure associated mortalities in this review, with statistically reproduced exposure levels to enhance the curve fitting, as illustrated in Supplementary Texts S1-S3. The exposure imputing had revealed high reliability, but the high uncertainties of the IER curves could still not be addressed, which should be attributed to the limited effective epidemiological evidences. Empirically, this approach would require sufficient studies to cover a wide range of exposure levels, which had been frequently adopted for particulate matter exposure researches,^70-72 but seldomly used for O₃-health studies.²⁰ The main probable reason might be that the population long-term O₃ exposure levels would not be as comparably distinguishable to the particulate matters. In addition, a reasonable prescribed TMREL would be necessary to establish the IER curves,³⁷ and hence the indeterminacy of the threshold level could exacerbate the uncertainties in the estimated concentration-response relationships. Therefore, more relevant studies on long-term O₃ exposure associated risks are urgently appealed for, based on which discussions, optimisations, or corrections on our enhanced exposure-response curve-fitting methodologies, are encouraged.

4.5 Hierarchical classification of diseases

The causes of mortalities analysed in our study followed hierarchical subordinate relationships, as the all-cause mortality consisted of cardiovascular diseases, respiratory diseases, cancer and other causes; chronic obstructive pulmonary disease belonged to respiratory category; and ischaemic heart disease, stroke, congestive heart failure and other cerebrovascular diseases all subordinated to cardiovascular symptoms. On this occasion, estimating all O₃-exposure induced mortalities could follow a bottom-up scheme by adding up subgroups of diseases. However, for the historical O₃-associated mortalities, GBD attributed all O₃-associated mortalities onto COPD-induced premature deaths,²⁰ which we thought were of spaces for further explorations. Long-term O₃ exposure had shown significant association with excess cardiovascular mortalities, and thus we should update the mortality estimations in further studies by including CVD altogether into consideration.

4.6 Application in mortality estimations

The widest applications of the estimated risk association strengths were to project how many people would be affected by long-term ambient O₃ exposure. For example, Malley et al. (2017) estimated 1.23 (95% UI: 0.85–1.62) million respiratory deaths attributable to O₃ exposure in 2010,⁷³ using the risk strength by Turner et al. (2016) as HR = 1.12 (95% UI: 1.08–1.16).⁵⁰ This estimation was much higher than the 2019 GBD report: 0.31 (95% UI: 0.15–0.49) million, as had been highlighted in another recent study,²⁵ which should be attributed to the use of high HR value among all included studies. We had also found some other studies using one singular HR value for population risk estimations,^{17, 74-77} but we would still encourage further relevant studies to consider multi-study pooled RRs, which could effectively reduce the potential biases from a single study. The adaptability of the pooled RRs could be verified from the coverage of exposure levels, as the 25 studies identified in our review had embraced a wide range of exposure concentrations (Supplementary Text S1) to encompass the global surface O₃ variability.²⁵ On the other hand, the leave-one-out sensitivity analyses (Table S6) had revealed the robustness of the meta-analysis results when including sufficient numbers of studies, which was a circumstantial reflection for the representativeness of the synthesised risk association strengths. The annual GBD reports were also based on the generalisability presumption of the synthesised epidemiological evidences, but cohort-based researches in the unstudied regions are always appealed for to provide more convincing discoveries.

4.7 Limitations

Although the total number of studied participants for risk pooling was adequately high to ensure the statistical power, the cohort-based O₃-health studies were factually rare according to our literature search, and thus long-term follow-up studies are urgently encouraged. Additionally, current literatures seldomly reported grouped RRs, which made meta-analyses by sub-categories (like gender, age, socio-economic status, smoking and alcohol history, etc.) unfeasible. Scarcity of credible evidences also restricted the effects of conventional approaches to construct exposure-response curves, and our methodological innovation would require further relevant studies for substantiation. The cross-metric linear conversion factors were estimated relying on observations from available sites, which however might not be sufficiently representative of the global residential areas, as observational sites in India, Africa, and Latin America were still sparse. With ever-increasing popularisation of the in situ monitoring networks, the cross-metric conversion factors might need calibration with more comprehensive observations, so that the pooled RRs should also be updated accordingly.

4.8 Further study suggestions

We suggest that further environmental epidemiology studies, especially long-term O₃ exposure related researches, clearly report i) the methodologies to obtain ambient O₃ concentrations, the spatiotemporal resolution, and prediction accuracy of the database; ii) the exposure metrics used for risk estimation; and iii) the statistical distribution of the O₃ exposure concentrations. The data-oriented methodologies to accomplish full spatial coverage ambient air O₃ concentrations for individual-level exposure assignment should be transparent as the construction credibility of air pollution concentration database should also be rigorously assessed, which were the foundation of epidemiological follow-up studies. We would advocate the report of exposure metrics in future O₃-health studies so as to avoid confusions when comparing the risks with literature and conducting meta-regression; and according to the recent consensus, warm-season average (6mDMA8) shall be preferred as epidemiological study metrics.¹⁹ We recommend future studies estimate risks with multiple O₃ metrics for reference; and describing the statistical distribution of the O₃ exposure levels is another suggested element to assess the reliability of risk estimation models, which can also be useful in exposure-response tendency exploration. We also propose future cohort studies estimate subgroup-specific RRs which can be conducive to identify the vulnerable populations.

Our review highlights a deficiency existing in current environmental health research literatures, that studies on long-term O₃ exposure health effects are still rather rare compared to particulate matter-based studies.⁷⁸ Also, the meta-analysis results might be geographically-biased, since large-scale O₃ exposure health risk studies till 2022 did not cover Asia, Africa or Latin America regions. However, there are some thriving cohorts such as the Multi-Country Multi-City (MCC) Collaborative Research Network covering over 22 countries or regions,⁷⁹ and the China Kadoorie Biobank (CKB) enrolling over 0.5 million people,⁸⁰ working on environmental exposure projects. We are optimistic that more research will come out to fill in the literature gap of multi-region population-based studies.