Zhoushan City is one of the 27 cities in the central area of the Yangtze River Delta in the northeast of Zhejiang Province. It is a coastal city with a subtropical monsoon climate, a warm winter, cool summer, mild and humid, with sufficient sunlight. By the end of 2018, the annual average PM_2.5 concentration in Zhoushan was 22 μg/m³; additionally, the proportion of days with good daily air quality (AQI) was 94.8%. We collected the average daily records of environmental pollutants (μg/m³) from the Zhoushan Environmental Monitoring Center from January 2014 to December 2020, including the average daily levels of PM_2.5, sulfur dioxide (SO₂), nitrogen dioxide (NO₂), carbon monoxide (CO), and O₃. Additionally, the daily average temperature, relative humidity, and atmospheric pressure of Zhoushan during the same period were obtained from the China Meteorological Data Sharing Service System (http://zshbj.zhoushan.gov.cn/index.html). Finally, the overall air quality data of Zhoushan City was obtained from the comprehensive correction of data collected from Putuo East Port, Dinghai Tain, Feng, and Lincheng New Area.

The population of Zhoushan has stably increased from 1.145 million to 1.176 million between 2014 and 2020 (19). This study selected the Zhoushan Hospital as the research site since it is the only tertiary, first-class general hospital in Zhoushan City and is preferred by residents due to its technology, professionalism, service, and advanced facilities and equipment.

Daily pediatric outpatient records from January 2014 to December 2020 were obtained from the Zhoushan Center for Disease Control and Prevention. Outpatient records for respiratory events were coded using J00–J99 and R04–R09 of the International Classification of Diseases, 10th edition. This study included children aged 0–18 years, regardless of history of hospitalization. All data were anonymized to protect participants' privacy. This observational study was approved by the Institutional Review Board of the Shanghai International Medical Center (SIMC-IRB No: 20210513).

A generalized additive model with natural splines was constructed to link data by date to explore the association between ambient air pollutants and daily pediatric respiratory outpatients. Additionally, the associations between daily pediatric respiratory outpatients, daily concentrations of ambient air pollutants (PM_2.5, PM₁₀, SO₂, O₃, NO₂, and CO), and meteorological factors (temperature, humidity, and atmospheric pressure) were analyzed. A basic model was constructed without air pollutants. Natural spline functions of time and meteorological parameters were incorporated into the model to adjust for confounding effects, with the day of the week (DOW) as a dummy variable. To prevent model multilinearity, we used Spearman's test to examine the correlation between daily concentrations of ambient air pollutants (20).

The model was examined using Akaike's information criterion (AIC). Additionally, partial autocorrelation functions were used to examine the autocorrelations of the residuals. Synthesizing the results from AIC, PACF, 3 df for the relative humidity, and 6 df for the mean temperature were added in our work as indicator variables. The model is described as follows.

Where: E(Y_t) is the number of expected pediatric outpatients at day t; β_i, the association between the log-relative rate of outpatient visits and air pollutant unit increase; X_i, the concentration of air pollutants on day t; DOW, the dummy variable; ns (time/temperature/humidity, 6/6/3) is the natural spline function for time, temperature, and humidity with 6/6/3 df.

Relative risks and their confidence intervals (CIs) were estimated to quantify the influence of each air pollutant on pediatric outpatients. All results are presented as relative risk (RR) and 95% CI of daily pediatric respiratory outpatients, which were associated with every unit increase in ambient air pollutant concentration.

Data were stratified according to seasonal patterns, namely, the cool season (November to March), the hot season (June to August), and the transition season (April, May, September, and October). In addition, we built two-pollutant models to evaluate the stability of PM_2.5 effects after adjustment for co-pollutants. Lag analysis was conducted to observe the aspect of time in revealing some delayed effects of air pollutants on pediatric respiratory outpatients. The following lag structures were used: single-day lag (0–5) and multi-day lag (01,03,05). Considering the lagged and usually non-linear relationship, we applied a distributed lag non-linear model (DLNM) to explore the influence of PM_2.5 on pediatric respiratory outpatients (21).

All statistical tests were performed using R software (version 4.0.4) with “mgcv” (https://github.com/cran/mgcv) and “DLNM” (https://github.com/gasparrini/dlnm) packages. Statistical significance was set at P < 0.05. The results are presented as the percentage change of pediatric respiratory outpatient visits per μg/m³ increase in PM_2.5/day.

This is one of the few studies exploring the relationship between the number of pediatric respiratory outpatient visits and air pollution in developing country. Previous studies, most in developed country, have discussed the impact of air pollutants, such as NO₂ and PM_2.5, on the development of asthma in children (23–25). Other studies have found that changes in air quality affect the development of respiratory diseases in children, with an incidence rate ratio of ~0.8 for NO₂ and PM_2.5 (26), consistent with other studies in north China (13, 14, 27, 28). However, the ocean currents, climate, and economic development in south coastal cities have led to variations in the type and components of pollutants in different places (29). Another US study found a difference in PM_2.5 levels between indoors and outdoors (29); however, we found that PM_2.5 levels were five times elevated. Therefore, three important points are emphasized in this study. First, this study fills the gap of studies focus on the relationship between air pollution and children respiratory diseases in southern Chinese cities. Second, this study emphasized the impact of air pollutants in coastal cities, represented by Zhoushan, and the children with respiratory diseases. Third, the outpatient data of children with respiratory diseases were used. This was a comprehensive overview of all possible variations in childhood respiratory diseases. We used ICD codes to classify diseases to avoid missing relevant diagnoses.

Long-term and short-term exposure to fine particulate matter, such as PM_2.5 and PM₁₀, may indirectly affect the central nervous system and gastrointestinal tract (30). The small size of the particles allows easy penetration of the body's immune barriers, causing acute and chronic diseases and even premature death (31–33). A synergy exists between these major air pollutants, which travel from the nasal cavity to the trachea or the lungs, where they interact with each other. The mucosal layer, epithelial tissue, or respiratory tract cells coordinate with immune cells such as dendritic cells, macrophages, and monocytes to protect the body from pollutants. However, air pollutants can interrupt the body's homeostasis and release cytokines, which can damage respiratory health (34, 35). The air pollutant-induced excessive release of cellular inflammatory factors and immune cells may cause tissue damage and lung inflammation. Additionally, the immune system may be targeted, leading to immune dysfunction due to damage to immune cells such as natural killer cells (36–38). In addition, organic compounds in air pollutants can release electrons to form superoxides, damaging the body's antioxidants and causing inflammation (39, 40). Other studies have shown that prolonged exposure to low doses of PM_2.5 causes mast cell overexpression, which leads to increased inflammatory cytokine and histamine production (41, 42). Overall, the number of outpatient visits was the most sensitive indicator of impaired respiratory function. In Yangcheng, China, a time-series analysis found a 1.69% increase in the number of outpatient visits for respiratory diseases caused for every 10 μg/m³ increase in PM_2.5−10 (43), which was consistent with that of eastern and western cities in China (44, 45). These findings prove the relationship between air pollutants and the number of outpatient visits.

Additionally, air pollution and temperature exhibit a synergistic effect. Studies have found that air pollutants in drought conditions can exacerbate chronic respiratory diseases and cause acute reactions (46, 47). Our study found differences in pollutant concentrations at different temperatures. Notably, carbon monoxide concentration was highest in the cold season but lowest in the hot season. This may be attributed to the status of Zhoushan city as a coastal city, where residents burn coal for outdoor heating in the cold season. There was no relationship observed between the mean annual air temperature and the annual trend of pollutant concentration, suggesting the presence of many unknown confounding factors between temperature change and pollutant concentration.

In this study, we found that short-term exposure (within 1 month) to PM_2.5 was significantly associated with increased pediatric respiratory outpatient visits, especially in the lag mode of 1–7 days. In addition, the lag patterns of 1–3 and 1–30 days showed significantly rising levels. Jiang et al. found a significant association between coarse particulate matter and a 5-day lag in respiratory and cardiovascular outpatient visits (43), Other studies have shown similar results, such as PM_2.5 being associated with a a-day lag in asthma visits and PM_2.5−10 being associated with 1-day lag mortality (48, 49). We observed that an increase in PM_2.5 concentration by 10 μm/m³ resulted in ~4 and 10% in pediatric respiratory clinic visits for a lag of 1–3 days and 1–7 days, respectively. This is much higher than the 1% to 3% observed in previous studies and may be attributed to the growth and development period of the children's respiratory system. During this period, children's lungs are more vulnerable to air pollutants compared to adults. The lag time we observed was between 1 and 7 days, which is similar to previous studies. Additionally, the lag model of pollutant synergy revealed that the combination of PM_2.5 and O₃ reached a maximum value of 0.7% at lag03, which reached 1.6% in the Yancheng study. This may be attributed to the relatively low level of O₃ and the weak coordination effect of PM_2.5 caused by Zhoushan city being a coastal city with a lower industrial level and vehicle number than Yancheng.

Lower concentrations of PM_2.5 caused a gentle decrease in slope (Figure 3); in contrast, an increasing slope was observed with increased PM_2.5 concentration. This finding may be due to the small size of PM_2.5; furthermore, it is mainly deposited in the lower respiratory tract, and a high concentration is required to reach the lungs. However, further randomized clinical trials are needed to validate these results.

The reason we choosed DLNM, which is a non-parametric model, based on following two points. Firstly, we discovered that the distribution of scatter plot is a non-linear wavy curve; Secondly, the DLNM takes account of both lag effects and the exposure-reaction non-linear relationship, and initially applied in epidemiological studies in 2006 (50). Multiple studies have suggested that the effects of air pollution and meteorological factors on outcomes have lag effects, while traditional linear models, which only study effects on certain time points without considering lag effects, are likely to produce high linearity, biasing the results (51, 52). Therefore, DLNM is a better choice to study the lag effects of exposure, especially in the analysis of air pollution. This study also use GAM with smooth spline function to fit variables in order to reduce the model risk caused by linear settings. Which means predictor variable in the model will be divided into multiple parts, and then fits each part separately through a polynomial function.

This study has several limitations. First, exposure error bias greatly affects time series analysis, limiting the extrapolation of this study's inferences. Second, we used outpatient data; therefore, misclassification bias was inevitable. Third, the cities included were located in a specific region, which may affect the generalizability of the results. Future studies using larger sample sizes over a wider range of regions are needed to validate the findings of this study.