Climate-Driven Trends in Pathogen Spectrum of HFMD

Background: Pathogen spectrum is vital during infectious disease epidemics. However, the role of pathogen spectrum during hand, foot, and mouth disease (HFMD) epidemics and the major contributing factors remains unknown. Methods: HFMD cases from 2016 to 2018 in Guangzhou, China were collected. The relationship between HFMD cases and pathogen spectrum, as well as the association between pathogen spectrum and climate factors were studied using general linear models. We transformed the pathogen spectrum to the isometric log-ratio (ILR) components included in the model. Additionlly, vaccination rates were adjusted in the climate driven models. Results: We observed seasonal trends in HFMD cases, pathogen spectrum, and climate factors. The model regressing case numbers on pathogen spectrum revealed negative associations with both the ILRs of CAV16 (RR = 0.725, P<0.001) and EV71 (RR = 0.421, P<0.001). The model regressing pathogen spectrum on driven factors showed that the trends for EV71 proportions were inversely related to vaccination rate (%, β = -0.152, P=0.098) and temperature (°C , β = -0.065, P=0.004). Additionally, the trends for CVA16 proportions were inversely related to vaccination rate (%, β = -0.461, P=0.004) and temperature (°C , β = -0.068, P=0.031). The overall trends for pathogen spectrum showed that EV71 decreased signicantly, while the trends for CVA16 increased annually. Conclusions: Our ndings suggest a potential pathway for climate factors - pathogen spectrum - HFMD cases. Our study is practical and useful for targeted prevention and control, as well as provides environmental-based evidence.


Introduction
Hand, foot, and mouth disease (HFMD) is a common and potentially fatal infectious disease in children under 5 years of age (U.S. Centers for Disease Control and Prevention, 2020). Deaths from HFMD occur infrequently, with an annual average of 191 deaths over the past ve years in China (Department of Disease Prevention and Control, 2020). In China, HFMD has caused substantial disease burden, and remained at the top of the list for Group C infectious diseases in terms of morbidity and mortality since the May 2008 outbreak in Anhui (Department of Disease Prevention and Control, 2020). In preschools, the presence of HFMD cases can lead to suspension of classes (Ministry of Education of the People's Republic of China, 2006). The parents/guardians of infected children are often nancially impacted due to income lost resulting from not working to care for their child. Therefore, HFMD remains a major societal concern.
Variation in pathogen spectrum is potentially an important contributor for the consistently high number of HFMD cases in China. A few serotypes of enteroviruses including coxsackie virus A16 (CAV16), enterovirus 71 (EV71), CAV4, 5, 9, and 10, and CBV 2 and 5 is potentially the cause of HFMD (Bruu, 2003).
In this study, pathogen spectrum was de ned as a variety of viruses are responsible for causing HFMD.
The pathogens that cause HFMD are RNA viruses which have a high mutation rate (Duffy, 2018). The EV71 monovalent vaccine, which was launched in 2016 (Mao et al., 2016), does not seem to have made a signi cant impact on the overall HFMD epidemic that may be due to the herd immunity not being achieved yet (Du et al., 2020). In addition, different virus types have different activity and transmission characteristics, which can cause cases of infection on different scales. Therefore, it is of major public health importance to understand the characteristics of the pathogen spectrum for HFMD prevention and control.
Climate factors may affect the HFMD epidemic by in uencing the spectrum of pathogens. The relationship between HFMD and temperature or relative humidity have been widely reported in previous studies (Du et  This study will examine the relationship between HFMD cases and pathogen spectrum. We also investigated the association between HFMD pathogen spectrum and climate factors to develop a greater understanding on the potential mechanism of climate factors -pathogen spectrum -HFMD cases.

Study area
Guangzhou city, one of China's three largest cities, is the capital of Guangdong province in South China with a population of 14,904,400 at the end of 2018 (Wikipedia, 2020). The area of Guangzhou city is 7,434.4, km 2 on both sides from 112°57′ to 114°03′ E longitude and 22° 26′ to 23° 56′ N latitude in southcentral Guangdong (Wikipedia, 2020). Guangzhou city has a humid subtropical climate, with a short mild, and relatively dry winter, while a long hot, and very wet summer (Wikipedia, 2020).

Data collection
HFMD is a noti able disease in China, therefore all probable and con rmed HFMD cases must be reported online to the surveillance system (China Information System for Disease Control and Prevention) within 24 hours of diagnosis using a standardized form. We collected data on HFMD cases from the surveillance system managed by China Center for Disease Control (CDC) in Guanghzou from January 2016 to December 2018. The CDC collects data for all HFMD diagnoses by law, including birth date, onset date, and virus types. The virus types of the earliest 5 cases (at least) per month from each sentinel site are determined using uorescence quantitative PCR. All HFMD cases in Guangzhou were included in the current analysis. Three monthly time-series for the proportions from three reported virus types were calculated as the pathogen spectrum of HFMD.
Vaccination data were obtained from the surveillance system managed by the Guangzhou CDC. It is important to note that there is currently only a vaccine against EV71. The CDC collects date of birth, vaccination date, and cumulative vaccine doses for all inoculations administered in Guangzhou. We computed the monthly vaccination rate by dividing monthly number of vaccinations by the population under 5 years of age in each year. HFMD pathogen spectrum and vaccination rate were linked by the calendar month.
Climate data including air temperature (°C) and relative humidity (%) were collected from the China Meteorological Data Service Center (CMDC, http://data.cma.cn/en). There are two national weather stations located in Guangzhou. Each station records the daily average air tempeature (°C) and relative humidity (%). Daily data from two stations were aggregated into monthly average data for Guangzhou and linked with HFMD pathogen spectrum by the calendar month.

Statistical analysis
To further analyze the pathogen spectrum of HFMD (i.e., three compositional time-series), we conducted isometric log-ratio transform (ILR) (Egozcue et al., 2003) for pathogen spectrum to obtain the ILR components. The ILR components can be analyzed without the limitation of collinearity and restriction limited to 0 ~ 1. Setting the compositional time-series of other viruses, CVA16 and EV71 as p i , i = 1,2,3, then the ILR components in each month with other viruses as the reference are as follows: By default, if one of p i is equal to 0, all ILR components for this month will set to 0.
First, we described the number of cases and pathogen spectrum using the time-series plots. Second, we examined the association between case numbers and pathogen spectrum by regressing case numbers on the ILR components using poisson regression. The relative risk (RR) of ILR component was extracted.
Third, the driven factors (e.g., vaccination rate, temperature, humidity) were described using the timeseries plots. Fourth, we used the univariate and multivariate linear regression, regressing the ILR components on the driven factors to examine the association between pathogen spectrum and driven factors. To control for the unmeasured time-varying confounding, we used the natural cubic splines of calendar time with 4 df (common multiple number for 2 epidemic peaks and 4 seasons) per year to remove the long-term trends and seasonality (Xiao et al., 2017). Finally, the tted and forecasted pathogen spectrum were conducted to validate the internal and external accuracy of the climatepathogen model. In the present study, all data management and statistical analyses were conducted using R version 4.0.2 (R Foundation for Statistical Computing, Vienna, Austria; https://www.rproject.org/). Table 1 shows the characteristics of HFMD cases by serotyped and unserotyped in Guangzhou, China from 2016 to 2018. A total of 185,838 HFMD cases were reported to the Guangzhou CDC surveillance system from January 2016 to December 2018. Most HFMD cases were diagnosed in children at one year of age (29.82%), among males (60.36%), and in mild (100.00%). A total of 5,067 (2.73%) cases completed a serotype-speci c antigen test. Among the cases serotyped, children at one year of age (31.40%), among males (61.26%), in mild (99.90%) had similar distributions as unserotyped cases.  Figure 1 presents the seasonal trends in case numbers and pathogen spectrum for HFMD from 2016 to 2018 in Guangzhou, China. At the beginning of each calendar year (January), the case numbers continuously decreased to a low level, while the proportions of EV71 remained the predominant pathogen. During the middle of each calendar year (approximatedly April to October), the case numbers continuously reached to a high level, while the proportion of other viruses always stayed the predominant pathogen. The trends for CVA16 was comparable to that of EV71, with the proportions increasing annually. The results of regressing case numbers on pathogen spectrum shows that case numbers were negatively associated with both the ILRs of CAV16 (RR = 0.725, P < 0.001) and EV71 (RR = 0.421, P < 0.001). Figure 2 presents the time trends for vaccination rates and climate factors. We observed the crude relationship between these factors and case numbers, as well as pathogen spectrum for HFMD. Since the launch of EV71 vaccine in 2016, the vaccination rate among children under 5 years of age increased slightly to approximately 50% at the end of 2018. Temperature had an expected seasonal trend that was analogous to the trends for HFMD case numbers. A positive relationship was observed between temperature and the proportions of other viruses, while a negative relationship was observed for CVA16 and EV71. Similar ndings were also observed for humidity. Table 2 shows the association between HFMD pathogen spectrum and vaccination rate, as well as climate factors derived from univariate linear regressions. The vaccination rate (linear regression coe cient (β) = -0.235, P = 0.029), temperature (β = -0.062, P < 0.001), and humidity (β = -0.041, P = 0.040) were inversely related to the trends for EV71 proportions. Similarly, the trends for CVA16 proportions was inversely related to vaccination rate (β = -0.562, P < 0.001) and temperature (β = -0.056, P = 0.024), while the humidity (β = -0.009, P = 0.754) was not statistical signi cant.  *: statistical signi cance with P < 0.05; †: statistical signi cance with P < 0.10. Table 3 shows the association between HFMD pathogen spectrum and vaccination rate, as well as climate factors derived from multivariate linear regression. The vaccination rate was inversely related to the trends for CVA16 (β = -0.461, P = 0.001) and EV71 (β = -0.152, P = 0.098) proportions. The temperature was inversely related to the trends for CVA16 (β = -0.068, P = 0.031) and EV71 (β = -0.065, P = 0.004) proportions. However, the humidity was not found statistically signi cant association with CVA16 (β = 0.041, P = 0.246) and humidity (β = 0.011, P = 0.632) proportions. In Table 4, each unit increase in vaccination rate was associated with all changes for HFMD pathogen spectrum (13.84% increase for other viruses, 8.77% increase for CVA16, and 5.07% decrease for EV71). Each unit increase in temperature was associated with all changes for HFMD pathogen spectrum (2.52% increase for other viruses, 0.76% decrease for CVA16, and 1.76% decrease for EV71). Each unit increased in humidity was associated with all changes for HFMD pathogen spectrum (1.10% decrease for other viruses, 0.80% increase for CVA16, and 0.30% decrease for EV71).  Figure 3 presents tted trends in pathogen spectrum for HFMD derived from the multivariate linear regression. The seasonality for the tted trends was consistent with the observed trends. The trends for other viruses decreased slightly, while EV71 decreased signi cantly. By contrast, the trends for CVA16 increased annually. Figure 4 presents the predicted performance and forecasted trends in pathogen spectrum for HFMD based on the multivariate linear regression. During the stepwise forecasting (Jun 2018 to Dec 2018), the result for only one month in Nov 2018 was incorrectly predicted. In the following six months, CVA16 would become the predominant pathogen and then the secondary pathogen in Jun 2019. As expected, EV71 would maintain the lowest pathogen during the forecast period.

Discussion
Using a time-series design, this study examined the relationship between HFMD cases and pathogen spectrum, as well as the association between pathogen spectrum and climate factors. Our ndings revealed that HFMD cases were negatively associated with the ILR components for both ILR CAV16 and ILR EV71 . Additionally, change in HFMD pathogen spectrum was associated with climate factors. Each unit increase in temperature will result in 0.76% and 1.76% decrease in CVA16 and EV71, while each unit increase in humidity will cause 0.80% increase in CVA16 and 0.30% increase in EV71. To the best of our knowledge, this study is the rst to reveal the role of pathogen spectrum in an endemic area.
Pathogen spectrum is an important factor to consider during a HFMD epidemic. The ndings in this study showed a statistical association between HFMD cases and pathogen spectrum. The substantial cause of diseases is always due to serotypes sporadically emerging (Pons-Salort and Grassly, 2018). For instance, around the year 2015 CVA6 emerged as the main cause of HFMD worldwide (Bian et al., 2015). Having a better understanding of pathogen spectrum would support the appropriate choice of serotypes in multivalent vaccines for targeting a disease. For example, we observed that the proportion of EV71 decreased annually. Therefore, a multivalent vaccine besides the present EV71 vaccine is urgently needed.
Climate factors may affect the HFMD epidemic by in uencing the spectrum of pathogens. Previous studies have described the relationship between climate factors and HFMD in terms of viral activity and outdoor activities (Du et al., 2019). Findings in the present study are more biologically plausible due to showing that the trends for CVA16 proportions was inversely related to temperature and positively related to humidity. Therefore, we suggest a potential pathway for climate factors -pathogen spectrum -HFMD cases.
The present EV71 vaccine is effective in prevention of HFMD. We found that the trends for CVA16 and EV71 proportions were inversely related to the vaccination rate. This nding is consistent with our previous studies assessing the effectiveness of EV71 vaccine (Du et al., 2019). While an overall 1-year e cacy of monovalent EV71 vaccines of 95% was reported, a bivalent EV71/CVA16 vaccines is potentially more cost-effective [14]. For this reason, an expanded program for EV71 vaccination is immediately needed.
The models based on ILR transform are suitable for compositional data, such as pathogen spectrum. Applying ILR transform to compositional data have been validated in the scienti c areas of medical imaging (S. et al., 2014) and geosciences (Buccianti et al., 2006). The pathogen spectrum of a certain disease is a classical compositional data. In addition, the ILR components can be used as both dependent and independent variables providing great exibility for its applications.
While this study offers several strengths, including the practical value of ILR transform and the potential pathway of climate factors -pathogen spectrum -HFMD cases, there are several important limitations worth mentioning. First, this study is essentially ecological in nature and thus unable to rule out potential ecological fallacies. However, the associations we observed could prompt other studies of causal design. Second, cases were randomly sampled for laboratory con rmation/viral serotyping but was not performed for all cases. However, this sampling information was representative of all HFMD cases due to the similar distributions between serotyped and unserotyped cases. Our ndings could be applied to settings with similar population contact patterns and climate conditions where the potential pathway of climate factors -pathogen spectrum -HFMD cases might be similar.

Conclusions
Our study contributes to the limited knowledge on quantifying the associations for HFMD casespathogen spectrum and pathogen spectrum -climate factors. The results of this study suggest a potential pathway for climate factors -pathogen spectrum -HFMD cases. This study can serve as a reference for studies on the associations between infectious diseases and environmental factors. Moreover, our ndings are useful for targeted prevention and control of HFMD.

Declarations
Availability of data and material The cases data and vaccination data are available from Guangzhou, China CDC (http://www.gzcdc.org.cn/), which were used under license and not publicly available. Authors' contributions ZD and YH participated in the design, performed data analysis and interpretation, and drafted the manuscript. WL, JX, ZY, and JL participated in the design and results interpretation, and helped to nalize the manuscript. ZZ and YH participated in the design and supervised the study, participated in interpretation of results, and helped to nalize the manuscript. All authors have read and approved the contents of the nal version.