key: cord-0965478-fp5jwrmm
authors: Ji, Xue-Yue; Huang, Li-Yuan; Song, Jia; Fei, Chun-Nan; Liu, Jun; Liu, He
title: Short-term effects of meteorological factors, air pollution, and sunspot on childhood hand, foot, and mouth disease in Tianjin, China: a new time series regression, 2014–2018
date: 2020-06-23
journal: Environ Sci Pollut Res Int
DOI: 10.1007/s11356-020-09794-x
sha: 123589feb77f17e27c7b93326aaae484bff34424
doc_id: 965478
cord_uid: fp5jwrmm

This study is aimed at defining the relationship between a set of environmental factors and childhood HFMD and then at estimating the related effect. The 16 environmental factors included meteorological, air pollution, and sunspot. A traditional TSR modified by using susceptible-infectious-recovery models and distribution lag nonlinear model was applied to estimate the short-term effects of daily environmental factors on children HFMD occurrence in 2014–2018 with adjustment of potential confounding factors. A total of 70,027 children aged 0–15 years with HFMD were enrolled. No significant effect was observed for daily sunspot numbers and average visibility. We found positive effects of the ambient average temperature, with an approximately m-shaped curve of the overall cumulative relationship, peaking at 25.6 °C with a relative risk (RR) of 1.45 (95% confidence intervals 1.21–1.73). The largest RR value of hot effect was achieved on the current day and then decreased by 2 days (total group, male group, and scatter group) or 1 day (female group and nursery group), and the effect lasted about 6 to 8 days from the lag 4 or lag 6 day. A greater association of temperature with HFMD for the female group and the scattered group was observed. This study suggests that ambient average temperature might be a risk factor for children HFMD in Tianjin. Further studies are warranted to confirm these findings.

Hand, foot, and mouth disease (HFMD), first reported in 1957 in New Zealand, is a worldwide common viral infectious disease caused mainly by coxsackievirus A16 (Cox A16) and enterovirus 71 (EV 71). HFMD is named after its general clinical characteristics, mild symptoms of febrile illness and rashes involving the skin of the hands, feet, mouth, and occasionally the buttocks and genitalia (Huang et al. 2016; Zhuang et al. 2015) . Infants and children younger than 5 years were particularly susceptible to HFMD (Zhu et al. 2012; Zhu et al. 2016) , the number of cases in children older than 5 years continually increased from 2001 (Zhuang et al. 2015) . Some patients can quickly develop potentially fatal nervous system and systemic complications, especially in cases associated with EV 71 (Xing et al. 2014) . Fortunately, the China Food and Drug Administration (CFDA) approved the first inactivated EV 71 whole virus vaccine for the prevention of severe HFMD in 2015 (Mao et al. 2016) . HFMD is highly contagious and spreads through direct contact or fecal-oral transmission through nasopharyngeal secretions such as saliva or nasal mucus (Huang et al. 2018) . In mainland China, since it has been made statutorily notifiable in May 2008, an average of about 2 million HFMD cases was reported every year, which has been one of the leading causes of child death and has raised substantial public health concerns in the Asia-Pacific region .

It is well known that weather plays a vital role in the transmission of many infectious diseases (Xiao et al. 2017 ). The seasonality of the incidence of HFMD has led to extensive research by epidemiologists in an attempt to quantify the association with meteorological variables (Huang et al. 2018; Sumi et al. 2017; Yu et al. 2019; Zhao et al. 2017 ). The results obtained from different regions on the effects of various meteorological factors were not entirely consistent, and the temperature-HFMD relationship was the most controversial. Individual studies have shown that the relationship between temperature and HFMD was negative (Samphutthanon et al. 2013) . Still, most studies have also confirmed the relationship was positive, whether linear (Gou et al. 2018 , Onozuka and Hashizume 2011 , Zhang et al. 2016a or nonlinear (Hii et al. 2011; Huang et al. 2016; Xiao et al. 2017; Xu et al. 2015; Yu et al. 2019; Zhu et al. 2016) . The non-linear relationship was reflected in the temperature-HFMD incidence curve, some of which were characterized by J-shape (Hii et al. 2011) , some of which were inverted V shape (Xu et al. 2015 , Yu et al. 2019 , Zhu et al. 2016 , and some of which are irregular (Huang et al. 2016; Xiao et al. 2017) .

In addition to temperature, other meteorological factors such as humidity, wind speed, and sunshine in some studies were used as influencing factors (Yu et al. 2019) , while others were considered confounding factors to control (Huang et al. 2016; Xu et al. 2015) . The inconsistency in results can be attributed in part to the diversity of methodologies and data sources. Still, it is also possible to suggest that temperature-HFMD relationships are modified by some other factors, such as air pollutants (Huang et al. 2018; Yu et al. 2019) or sunspot (Towers 2017 ) and location-specific variables, that contribute to the influence of temperature on the incidence of HFMD (Xiao et al. 2017) . Of particular concern is particulate matter, such as PM 2.5 or PM 10 ; that particulate matter having an aerodynamic diameter of 2.5 μm or 10 μm can enter the human respiratory system, even penetrate lung cells, and enter the circulatory system, resulting in serious health problem, especially in children. Yu et al. found that low PM 2.5 (< 17 μg/m 3 ) and high ozone (> 186 μg/m 3 ) exerted a particularly protective effect on HFMD incidence (Yu et al. 2019) . Huang et al. found that the highest relative risk (RR) with HFMD for PM 10 was 80 mg/m 3, and the lag effect was 5 days (Huang et al. 2018) . In contrast, it has also been suggested that the PM 10 level showed no relationship to the incidence of HFMD (Huang et al. 2016) . Therefore, the impact of exposure to environmental factors on the incidence of HFMD has not been fully determined.

The methodologies, used by most of the above studies, were time series regression (TSR), auto-regressive integrated moving average (ARIMA) models, and distributed lag nonlinear model (DLNM), which were originally developed to estimate the associations of weather and air pollution with all-cause mortality (Gasparrini 2016 , Thurston and Ito 2001 , Yu et al. 2011 , Zhang et al. 2016b or morbidity (Vicedo-Cabrera et al. 2019 ) in places where this is overwhelmingly due to non-infectious diseases ) (e.g., cardiovascular diseases (Masselot et al. 2018) ). Using the above methodologies to study the association between environmental factors and infectious diseases, not only the characteristics of infectious diseases, such as pathogenesis, immune population, and strong autocorrelation, but also the adaptation of mathematical modeling must be considered .

Distributed lag non-linear model (DLNM), generalized by Armstrong (Armstrong 2006) , is a class of models with different options for the functions applied to model nonlinearity and distributed lag effects. The susceptible-infectiousrecovery (SIR) model is based on immune to re-infection and other known mechanisms for the dynamics of immunity and transmission among the population . The purpose of this study is to modify traditional TSR practice that includes using SIR models and DLNM to estimate the short-term effects of environmental factors on the incidence of HFMD in Tianjin, a coastal city in China.

Tianjin (longitude 116°430′ to 118°040′; latitude 38°340′ to 40°150′) is the largest coastal city in northern China (Fig. 1) . As one of the four municipalities under the direct administration of the central government of China, Tianjin lays southeast of Beijing (the capital city) and borders Hebei Province, located along the coast of Bohai Gulf, comprising an area of 11,917.3 km and a total population of 15,621,200 as of 2016 estimation.

Tianjin has four divided seasons, monsoon-influenced climate-dry and windy springs, hot and rainy summers (due to monsoons), cool and pleasant autumns, and dry and cold winters (due to vast Siberian anticyclones). July tends to be the hottest month, while January is the coldest.

Daily cases of HFMD in Tianjin, China, on 1 January 2014 and 31 December 2018 were obtained from Tianjin Centers for Disease Control and Prevention based on the date of onset and current residence. According to the requirements of the "Standard management of infectious disease information report" issued by the National Health and Family Planning Commission of China, the infectious diseases to be reported, including HFMD, should be reported to the national infectious disease surveillance system via the Internet by a unified format, a prescribed time limit, and a specific process. The clinical diagnosis of HFMD patients was based on the "Guidelines for the diagnosis of hand, foot and mouth disease" issued by the Ministry of Health of China (National Health Commission of the People's Republic of China 2018): fever, or without fever, with hand, foot, mouth, or buttock rash. Cases without rash should not be diagnosed as HFMD in clinical practice.

Daily data on average temperature (T,°C), maximum temperature (TM,°C), minimum temperature (Tm,°C), atmospheric pressure at sea level (SLP, hPa), average relative humidity (H, %), total rainfall and/or snowmelt (PP, mm), average visibility (V, km), average wind speed (W, km/h), and maximum speed of wind (WM, Km/h) were crawled from the tutiempo.net (https://en.tutiempo.net). Daily concentrations of particular matters up to 10 or 2.5 μm in size (PM 10 or PM 2.5 , μg/m 3 ), sulfur dioxide (SO 2 , μg/m 3 ), carbonic oxide (CO, mg/m 3 ), nitrogen dioxide (NO 2 , μg/m 3 ) , and ozone (O 3 , μg/m 3 ) were crawled from aqistudy.cn (https://www.aqistudy.cn). Daily data on sunspot was downloaded from the website http://www.sidc.be. Tianjin has only started monitoring PM 2.5 since 2014, so our research data began from that year.

Statistical analysis was carried out in three steps. The first step was to organize the data, including data cleansing and integration of HFMD with environmental factor data. The second step was a descriptive analysis of all data, with a Spearman correlation analysis daily of HFMD cases and environmental factors. The third step was based on TSR, combined with the principle of DLNM and SIR model to evaluate the associations of environmental factors and HFMD.

The quasi-Poisson regression model that allows overdispersion of data was built as follows.

It is a traditional TSR analysis that seeks to identify how an observed time-varying variable of interest such as temperature x t explains a series of the observed HFMD cases on day t, where f(t) is a smooth function for avoiding confounding by long-term trends and seasonality, α is the model intercept, and β are regression coefficients. Other measured risk factors are denoted as z p,t .

The SIR model is based on the mechanisms for the dynamics of immunity and transmission among the population. The model simplified from Koelle et al. (Koelle et al. 2005; Koelle and Pascual 2004) can be written as follows:

Equation (2) is the result of integrating Eqs.

(3) and (4) and taking logarithms and making a Taylor series approximation (Koelle and Pascual 2004) .

where S t is the susceptible population size, estimated from subtracting the sum of fractions of past HFMD cases where fractions of immune κ i are assumed to steadily decline with the intervening time step t − i. N t is the total population size, θ t is virus transmissibility at time t. δ and γ are parameters related to the type of mixing between individuals. m is the entire duration of immunity (in time steps); ε t is multiplicative noise . Since many previous studies have shown that the meteorological indicator-HFMD relationship may be non-linear and considering the incubation period of infectious diseases, there may be a delay in the correlation between them, so we incorporate the DLNM into the TSR.

Eq.(5) is the final model, combined with DLNM and SIR, where μ t ≡ Log(Y t ); Y t is a series of the observed HFMD cases on day t, the day of observation from 1 to 1826, assumed to arise from a distribution belonging to the log-link family. l is the lag days. x t, l is a matrix obtained by applying cross-basis function in the dlnm (Gasparrini 2011) to examine variables that were closely related to the incidence of HFMD, where ns() is a natural cubic spline. z t represents the several other variables that should be controlled due to their modifying effect on HFMD incidence. After model selection and sensitivity analysis, only the average temperature was significantly associated with HFMD incidence; thus, other environmental factors were removed from the model. The cross-basis function between the average temperature and HFMD was to use a natural cubic spline with 4 degrees of freedom (df) for the exposure-response relationship and natural cubic splines with 4 df for the lag-response relationship. Spline knots were placed at equal spaces in each variable range and at equal intervals in the log scale of lags using the default setting of dlnm. The maximum lag time was set to 16 days. Time is the indicator variable, using the natural cubic splines with 10 df per year to remove the long-term trends and seasonality. DOW t is day of the week on the day. Holiday is a binary variable that is "1" if day t was a holiday (including the summer and winter vacations for schools and kindergartens and national public holidays), and σ is the coefficient. Log(Y t-1 ) was described in Eq. (2), which was the logarithm of lagged outcome counts. It suggested that for infectious diseases, rather than including lagged residuals as a term in the model, it matches the likely mechanism better to include the log(Y t-1 ) .

The median of average temperature was as the reference value for estimating RR. Stratified analysis was conducted by gender and care style. The care style group was divided into live scattered (mostly 0-3 years old), nursery (mostly 3-6 years old), and student (mostly 6-15 years old) according to daily activities.

The determination of all model parameters was based on previous studies (Huang et al. 2016; Huang et al. 2018; Yu et al. 2019) , with the Akaike information criterion for quasi-Poisson (Q-AIC) as the main indicator and comprehensive consideration by model sensitivity analysis. Sensitive analyses were performed by changing the maximum lag days (4-16 days) for environmental factors. The range of df for seasonal and long-term trends was varied from 1 to 15 and 3 to 6 for environmental factors and lag space. The R code includes a function for determining the number of knots in cross-basis and a function for calculating QAIC, which is the code exposed by Gasparrini et al. (2017) .

All statistical analyses were carried out in R 3.6.1 (R Core Team 2019) with the main packages glmulti 1.0.7.1 (Calcagno and Mazancourt 2010) for model selection and dlnm 2.3.9 (Gasparrini 2011) for exploring and evaluating the delayed effects. R itself and all packages used are available from CRAN at http://CRAN.R-project.org/. A P value of less than 0.05 was considered to be statistically significant. Table 1 shows the daily and total cases of HFMD by gender and care style. From 1 January 2014 to 31 December 2018, 70,027 HFMD cases aged 0-15 years were reported in Tianjin, of which 57.3% were scattered children, 34.8% were nursery children with centralized care, and 7.9% were students. The gender ratio was 1.46:1. The highest number of daily HFMD cases was 217, with a median of 23. Table 2 illustrates the descriptive statistics for daily data on environmental factors during the study period. The average daily values of main meteorological factors, namely, average temperature, atmospheric pressure at sea level, average relative humidity, average visibility, and average wind speed, were 15.7°C (range − 14.2~35°C), 1017 hPa (range 995.1 1043.9 hPa), 53% (range 10~99%), 15.4 km (range 0.33 0 km), and 9.8 km/h (range 3.5~25.2 km/h). The average daily concentrations of main air pollution factors, PM 2.5 , PM 10 , and O 3 , were 54 μg/m 3 (range 0~383.4 μg/m 3 ), 92.7 μg/m 3 (range 0~984 μg/m 3 ), and 86 μg/m 3 (range 32 88 μg/m 3 ). Table 3 reveals that the three daily temperature measures were strongly correlated with daily HFMD cases (r > 0.7). In addition to sunspot numbers and average visibility, the correlation between the daily cases of HFMD and other environmental factors was statistically significant (P < 0.001). Figure 2 shows the monthly distribution of HFMD from 2014 to 2018 in Tianjin, China, indicating that the number of HFMD in June and July (about 100 cases per month) was higher than that in other months. Figure 3 presents the time series distribution of total cases and environmental factors from 2014 to 2018 in Tianjin. Total cases, three temperature factors, air pressure, relative humidity, and ozone showed a distinct seasonality. We found that between 2014 and 2016, the peak of HFMD incidence occurred before June, while the daily average temperature peak was the opposite. In 2017 and 2018, the peak of HFMD incidence and temperature occurred after June, and the peak in 2017 was significantly lower than that in other years.

The number of cases in the student group was small (7.9% of the total), and the cross-basis of temperature-cases-lag was not statistically significant in the model. Therefore, we did not estimate the temperature-HFMD effect and lag effect of the student group. Figure 4 displays the general pattern of RR, as a function of temperature and lag, by showing three-dimensional plots of RR along temperature and 16 lag days. Overall, the estimated effects of temperature on HFMD incidence were non-linear, with higher relative risks at hotter temperatures.

The estimated associations represented in the RR scale, between temperature and HFMD in Tianjin, are illustrated in Fig. 5 . All panels display the overall cumulative exposureresponse curves, interpreted as the risk cumulated over the entire lag period of 0-16 days. Except for the nursery group in Fig.5e , the curve patterns of other groups (Fig. 5a-d) presented an m shape. The first shorter peak all appeared at about 1°C, with no statistical significance. The second appeared at about 25°C, with a higher peak and statistical significance. Also shown in Fig. 5a is a histogram of the daily average temperature, with the reference temperature (15.7°C) and the cut-off values for defining extreme cold (− 6.3°C) and Figure 6 estimates the risks of different percentile average temperatures relative to the median temperature (15.7°C) for total and group-specific HFMD cases along with the lags. Regardless of the group-specific or the total group, the low temperatures (1st and 25th) have smaller RRs than the high temperatures (75th and 99th). No significant associations between low temperatures and HFMD incidence were observed at different lag days. The exposure-response curves exhibited by higher temperature and HFMD incidence were the same in each group. The largest RR value of hot effect(75th) was achieved on the current day and then decreased by 2 days (total group, male group, and scatter group) or 1 day (female group and nursery group), and the effect lasted about 6 to 8 days from the 4th or 6th lag days. The temperature of each group in the peak column of Fig. 6 was the temperature at which the overall relative risks were the largest during the 0-16 lag days, as shown in Fig. 5 . Table 4 summarizes the cumulative effects of the different average temperatures at different lag days. For example, 25.6°C, for the total group, has the highest RR value, 1.45 (1.21-1.73), over the maximum lag days (0-16 days), while 26.8°C was the temperature that had the highest RR over all combinations of lag days, which here refers to the RR value of 1.51 (1.29-1.77) at the 0-13 lag days. 

In this paper, we conducted a time series analysis to examine the short-term effects of temperature with gender and care style specific on HFMD among children (0-15 years) in Tianjin, China, during 2014-2018. As some previous studies have shown, data can be better explained by the non-linear relationship and lag effects between temperature and HFMD. It further confirms that temperature can affect the transmission of HFMD, and the non-linear behavior of the temperature-HFMD relationship means that the effect of temperature on HFMD may involve complex mechanisms (Xiao et al. 2017 ). An interesting finding in our study is that the exposureresponse curve is approximately m-shaped, which is inconsistent with previous studies. Most of the previous studies presented a J shape (Hii et al. 2011) or inverted V shape (Yu et al. 2019; Zhu et al. 2015; Zhu et al. 2016) . A study in Beijing, a city adjacent to Tianjin, also showed that the temperature-HFMD curve was inverted an V shape, with the largest RRs at 25.0~27.5°C (Xu et al. 2015 ) over 13 days. Huang et al. found that in Ningbo, a coastal city similar to Tianjin, the curve of temperature-HFMD was m-shaped, although it was not obvious but also bimodal, and the highest peak appears at 31°C (Huang et al. 2016) . Although both of our studies in Ningbo showed a double-peak curve, the first shorter peak was not statistically significant. In this respect, it is consistent with most studies that the temperature-HFMD curve exhibits an inverted V shape, especially in the nursery group of this study. The RRs increased with the increment of temperature from the median temperature of 15.7°C, and it reached the peak (RR 1.45) at 25.6°C. It then began to decrease the groupspecific RRs that followed the same trends as the total RRs, with the peak at 25.6°C (male, RR 1.46), 26.1°C (female, RR 1.60), 25.9°C (scattered, RR 1.60), and 24.5°C (nursery, RR 1.31), respectively. There is evidence of a positive correlation between temperature and host activity or excretion of enterovirus (Belanger et al. 2009 ), so accompanying warming may increase the chances of susceptible individuals being exposed to infectious sources or contaminated environments (Xiao et al. 2017) . However, it is well known that temperature and ultraviolet radiation are the two main factors leading to the inactivation of enterovirus (Bertrand et al. 2012) . Thus, extremely hot temperatures may shorten the survival time of enteroviruses in the environment and thus reduce their chances of transmission back to the hosts. It could explain the risk reduction at very hot temperatures.

We found a greater association of temperature with HFMD for the female group and the scattered group. To date, few studies have examined the gender and care style effects of temperature on HFMD. Fortunately, the conclusions of limited research were consistent with our findings (Huang et al. 2018; Xu et al. 2016; Zhu et al. 2015) . Although the male-tofemale sex ratio was 1.46:1 in this study, it showed that female children were more sensitive to HFMD with temperature, consistent with the studies by Huang et al. (2018) and Xu et al. (2015) . The inconsistency of sensitivity to temperature in male and female children might be due to the different gender resistance, exposure opportunities, and the number of peer infections. There are fewer cases among female children's companions, which significantly highlights the effect of Fig. 4 Relative risks of daily HFMD by daily average temperature along 16 lag days for all groups in Tianjin, adjusting for DOW, holidays, autocorrelation, seasonal trend, and long trend temperature changes. Compared with nursery children, the incidence of HFMD in scattered children was more strongly related to ambient temperature, probably because of younger age (0-5 years old), no good health habits, and poor body resistance.

Currently, most research was focused on temperature and HFMD. Still, the Spearman correlation results in our study showed that in addition to the average visibility and sunspots unrelated to HFMD, the correlation between the other 14 environmental factors and the daily cases of HFMD was statistically significant. However, through repeated verification, after considering the lag effect, controlling the long-term trend and the seasonal impact, only the average temperature was included in the final model. Yu et al. proposed that average temperature, wind speed, precipitation, PM 2.5 , and O 3 were significantly correlated with the incidence of HFMD, and their effects varied with different lag days and subgroups (Yu et al. 2019) . Taking precipitation and average relative humidity as examples in meteorological factors, considering that the virus may live longer in a humid or watery environment and can tolerate various factors leading to virus inactivation (such as temperature, salinity, pH) (Salo and Cliver 1976) , some studies have shown that it is correlated with the incidence of HFMD.

Huang et al. showed a positive correlation between relative humidity and the incidence of HFMD, observing RR peak at a lag of 6 days (Huang et al. 2016) . Studies conducted in Hefei and Shenzhen have shown that extremely high humidity levels and wind speeds increase the risk of HFMD, possibly due to the favorable attachment of the virus to the surfaces of various objects (such as toys) with high humidity and the transmission of particles containing enterovirus through the Fig. 5 The overall relative risks of average temperature(°C) for total and group-specific HFMD cases over 16 days wind (Wong et al. 2010) . But most studies have controlled for these meteorological factors as confounding factors in the model (Xiao et al. 2017; Xu et al. 2015; Zhu et al. 2015) . Xiao et al. found that the temperature-HFMD relationships varied from region to region. Rainfall may play an important role in modifying the temperature-HFMD relationship (Xiao et al. 2017) . The Tianjin city we studied is a city with less rain. During the 5 years from 2014 to 2018, 79.2% (1446/1826) of the days were without rain, and 15.8% were light rain (< 9.9 mm/24 h). Relative humidity has a certain periodicity from the distribution pattern, and the day-to-day variation was not large. Therefore, the TSR cannot explain the incidence of HFMD from the day-to-day changes of the exposure (e.g., relative humidity, rainfall).

Interestingly, our study further confirmed that temperature, wind speed, visibility, rainfall, and O 3 were positively correlated with the incidence of HFDM. Conversely, PM 2.5 , PM 10 , SO 2 , CO, and NO 2 were negatively correlated. Taking PM 10 Fig. 6 The relative risks of average temperature (°C) for total and group-specific HFMD cases among different lags and PM 2.5 in the air pollution factors as an example, it is a hot issue in current environmental epidemiology research, mainly focusing on research on chronic diseases (Maji et al. 2017) and mortality (Choi et al. 2018; Renzi et al. 2018) . Whether it can affect the incidence of infectious diseases, especially HFMD, has also attracted the attention of relevant scholars. Yin et al. found that PM 10 -HFMD association was non-linear, with a peak value of 101-218 μg/m 3 , and male children were more sensitive to PM 10 effect (Yin et al. 2019 ). According to Huang et al.' s studies, there is no statistically significant relationship between PM 10 and total HFMD (Huang et al. 2016 ). Still, a certain relationship can be observed between PM 10 and female HFMD cases (Huang et al. 2018) . Zhong et al. reported that improving the air quality, especially decreasing PM 2.5 and PM 10 , could decrease the risk of HFMD outbreaks (Zhong et al. 2018) . These discrepancies can be attributed to a variety of reasons, such as methodological heterogeneity (e.g., data metrics and statistical analysis methods), as well as different climates, geographic conditions, economic factors, and air pollution control between cities. The Air Pollution Action Plan released in 2013 should be China's most influential environmental policy in the past 5 years. It is aimed at significantly improving air quality by setting PM 2.5 targets and requires a significant reduction between 2013 and 2017-of 25% in Tianjin. It can be seen from our study that the average PM 2.5 concentration in the whole year decreased from 1.43 (1.14-1.79) e a X = 0-13 lag day; b X = 0-12 lag day; c X = 0-15 lag day; d X = 0-15 lag day; e X = 0-12 lag day *The temperature that had the highest RR value over the maximum lag days (0-16 days) + The temperature that had the highest RR value over X lag days 86.5 μg/m 3 in 2014 to 50.5 μg/m 3 in 2018, a decrease of 41%, exceeding the target. The human intervention characteristics of air pollution factors may also be unrelated to the incidence of HFMD, which needs further research to confirm. Our study may be the first to explore the relationship between sunspot numbers and the incidence of HFMD. Sunspot activity at the maximum or minimum stage can significantly affect the earth's climate, causing extreme weather events such as heat, drought, and severe cold (Gachari et al. 2014) . A study by a colleague of the authors suggested that the twin peaks of the sunspot cycle in 2002 and 2012 were associated with the outbreak of severe acute respiratory syndrome (SARS) in China and Middle East respiratory syndrome (MERS) occurred in the Middle East (Qu and Wickramasinghe 2017) . The sunspot number was considered an excellent factor to predict the scope and the extent of severe diseases for humans during the maximal sunspot number (Kim 2019) . However, the sunspot number was not associated with the incidence of HFMD from our study. It may indicate that sunspot numbers are more suitable for the prediction of pandemic or outbreaks of new infectious diseases worldwide.

This study has three major strengths. First, methodologically, to our knowledge, this was the first study to examine the association between environmental factors and HFMD by combining SIR models with DLNM-based traditional TSR. Second, in terms of breadth of variables, we had explored 16 related factors, including nine weather factors, six air pollution factors, and one sunspot indicator. Third, in terms of stratification, we did not routinely stratify by age. Instead, we considered the characteristics of individual activities and grouped them according to different ways of care.

Our study should consider several limitations. The first was data quality issues, resulting in potential bias caused by under-reporting and repeated cases. In most cases, the use of clinical diagnostic methods rather than laboratory diagnostics was another source of reporting bias. The weaknesses of surveillance data were difficult to avoid in all cities. Because time series analysis was based on their previous observations, the impact of under-reporting bias should be limited (Xiao et al. 2017) . The second was that the small number of HFMD cases in student children (7.9% of the total) restricted us from analyzing the temperature-HFMD effects in the student groups. The third was that studies have shown that adjusting seasonality variation alone to control autocorrelation is not sufficient , so we included autoregressive item log(Y t−1 ) in the model. Also, the number of immunity to HFMD causes variation in the susceptible population that should be considered, but this information does not always exist.

In summary, our study provides an estimate of the temperature-HFMD relationship based on large datasets in Tianjin, China, including 5 years of daily HFMD cases, meteorological, air pollution, and sunspot. Only the temperature showed a non-linear relationship with the incidence of HFMD with the highest relative risk at the temperature range of 24.5°C (over 16 lag days in the nursery group) to 27.1°C (over 13 lag days in the male group). A greater association of temperature with HFMD for the female group and the scattered group was observed. Our results should help deepen the understanding of environmental factors, especially the temperature on the influence of HFMD. Also, this evidence could support the related public health decisions, such as establishing weather-based disease warning systems. Future studies should also take into account factors such as the immune population, spatiotemporal trends, temperature change between neighboring days, and socioeconomic.

Models for the relationship between ambient temperature and daily mortality

Influence of weather conditions and season on physical activity in adolescents

The impact of temperature on the inactivation of enteric viruses in food and water: a review

Glmulti: an R package for easy automated model selection with (generalized) linear models

Temporal variability of short term effects of PM10 on mortality in Seoul

Sunspot numbers: implications on Eastern African rainfall

Distributed lag linear and non-linear models in R: the package dlnm

Modelling lagged associations in environmental time series data

A penalized framework for distributed lag non-linear models

Different responses of weather factors on hand, foot and mouth disease in three different climate areas of Gansu

Short term effects of weather on hand, foot and mouth disease

Effects of meteorological parameters and PM10 on the incidence of hand, foot, and mouth disease in children in China

Impact of PM10 and meteorological factors on the incidence of hand, foot, and mouth disease in female children in Ningbo, China: a spatiotemporal and time-series study

A systematic review of methodology: time series regression analysis for environmental factors and infectious diseases

Time series regression model for infectious disease and weather

Spanish Flu, SARS, MERS-CoV by CO2 emission and maximal sunspot number

Disentangling extrinsic from intrinsic factors in disease dynamics: a nonlinear time series approach with an application to cholera

Refractory periods and climate forcing in cholera dynamics

Burden of disease attributed to ambient PM2.5 and PM10 exposure in 190 cities in China

EV71 vaccine, a new tool to control outbreaks of hand, foot and mouth disease (HFMD)

Aggregating the response in time series regression models, applied to weather-related cardiovascular mortality

Guidelines for hand, foot and mouth disease diagnosis and treatment

The influence of temperature and humidity on the incidence of hand, foot, and mouth disease in Japan

SARS, MERS and the sunspot cycle

R: a language and environment for statistical computing

Shortterm effects of desert and non-desert PM10 on mortality in Sicily

Effect of acid pH, salts, and temperature on the infectivity and physical integrity of enteroviruses

Spatiotemporal distribution and hotspots of hand, foot and mouth disease (HFMD) in northern Thailand

Association between meteorological factors and reported cases of hand, foot, and mouth disease from 2000 to 2015 in Japan

Epidemiological studies of acute ozone exposures and mortality

Sunspot activity and influenza pandemics: a statistical assessment of the purported association

Hands-on tutorial on a modeling framework for projections of climate change impacts on health

Human enterovirus 71 and hand, foot and mouth disease

The effect of meteorological factors on adolescent hand, foot, and mouth disease and associated effect modifiers

The exposure-response relationship between temperature and childhood hand, foot and mouth disease: a multicity study from mainland China

Hand, foot, and mouth disease in China, 2008-12: an epidemiological study

Non-linear association between exposure to ambient temperature and children's handfoot-and-mouth disease in Beijing

Impact of temperature variability on childhood hand, foot and mouth disease in Huainan

Analysis of the effect of PM10 on hand, foot and mouth disease in a basin terrain city

Short-term effects of meteorological factors and air pollution on childhood hand-foot-mouth disease in Guilin

The shortterm effect of ambient temperature on mortality in Wuhan, China: a time-series study using a distributed lag non-linear model

Short-term effects of meteorological factors on hand, foot and mouth disease among children in Shenzhen, China: non-linearity, threshold and interaction

Age patterns and transmission characteristics of hand, foot and mouth disease in China

Impact of weather factors on hand, foot and mouth disease, and its role in short-term incidence trend forecast in Huainan City

Forecasting hand, foot, and mouth disease in Shenzhen based on daily level clinical data and multiple environmental factors

Retrospective study of the incidence of HFMD and seroepidemiology of antibodies against EV71 and CoxA16 in prenatal women and their infants

The impact of ambient temperature on childhood HFMD incidence in inland and coastal area: a two-city study in Shandong Province, China

Assessment of temperature effect on childhood hand, foot and mouth disease incidence (0-5 years) and associated effect modifiers: a 17 cities study in Shandong Province, China

Epidemiological research on hand, foot, and mouth disease in mainland China

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Acknowledgments The authors thank Tianjin Centers for Disease Control and Prevention for providing data.

Ethical approval The Institutional Ethics Review Board of Tianjin Centers for Disease Control and Prevention approved this study.Competing interests The authors declare that they have no conflict of interest.Abbreviations HFMD, hand, foot, and mouth disease; TSR, time series regression; Cox A16, coxsackievirus A16; EV 71, enterovirus 71; CFDA, China Food and Drug Administration; PM 2.5 , particulate matter having an aerodynamic diameter of 2.5 μm; PM 10 , particulate matter having an aerodynamic diameter of 10 μm; RR, relative risk; ARIMA, auto-regressive integrated moving average; DLNM, distributed lag non-linear model; SIR, susceptible infectious recovery; T, average temperature; TM, maximum temperature; Tm, minimum temperature; SLP, atmospheric pressure at sea level; H, average relative humidity; PP, total rainfall and/or snowmelt; VV, average visibility; W, average wind speed; WM, maximum speed of wind; SO 2 , sulfur dioxide; CO , carbonic oxide; NO 2 , nitrogen dioxide; O 3 , ozone; df, degrees of freedom; CI, confidence intervals; Q-AIC, quasi-Akaike information criterion; SARS, severe acute respiratory syndrome; MERS, Middle East respiratory syndrome