Spatiotemporal patterns and risk mapping of provincial hand, foot, and mouth disease in mainland China, 2014–2017

Background Hand, foot, and mouth disease (HFMD) has remained a serious public health threat since its first outbreak in China. Analyzing the province-level spatiotemporal distribution of HFMD and mapping the relative risk in mainland China will help determine high-risk provinces and periods of infection outbreaks for use in formulating new priority areas for prevention and control of this disease. Furthermore, our study examined the effect of air pollution on HFMD nationwide, which few studies have done thus far. Methods Data were collected on the number of provincial monthly HFMD infections, air pollution, meteorological variables, and socioeconomic variables from 2014 to 2017 in mainland China. We used spatial autocorrelation to determine the aggregate distribution of HFMD incidence. Spatiotemporal patterns of HFMD were analyzed, risk maps were developed using the Bayesian spatiotemporal model, and the impact of potential influencing factors on HFMD was assessed. Results In our study, from 2014 to 2017, the HFMD annual incidence rate in all provinces of mainland China ranged from 138.80 to 203.15 per 100,000 people, with an average annual incidence rate of 165.86. The temporal risk of HFMD for 31 Chinese provinces exhibited cyclical and seasonal characteristics. The southern and eastern provinces had the highest spatial relative risk (RR > 3) from 2014 to 2017. The HFMD incidence risk in provinces (Hunan, Hubei, and Chongqing) located in central China increased over time. Among the meteorological variables, except for the mean two-minute wind speed (RR 0.6878; 95% CI 0.5841, 0.8042), all other variables were risk factors for HFMD. High GDP per capita (RR 0.9922; 95% CI 0.9841, 0.9999) was a protective factor against HFMD. The higher the birth rate was (RR 1.0657; 95% CI 1.0185, 1.1150), the higher the risk of HFMD. Health workers per 1,000 people (RR 1.2010; 95% CI 1.0443, 1.3771) was positively correlated with HFMD. Conclusions From 2014 to 2017, the central provinces (Hunan, Hubei, and Chongqing) gradually became high-risk regions for HFMD. The spatiotemporal pattern of HFMD risk may be partially attributed to meteorological and socioeconomic factors. The prevalence of HFMD in the central provinces requires attention, as prevention control efforts should be strengthened there.

Background: Hand, foot, and mouth disease (HFMD) has remained a serious public health threat since its first outbreak in China.Analyzing the province-level spatiotemporal distribution of HFMD and mapping the relative risk in mainland China will help determine high-risk provinces and periods of infection outbreaks for use in formulating new priority areas for prevention and control of this disease.Furthermore, our study examined the e ect of air pollution on HFMD nationwide, which few studies have done thus far.
Methods: Data were collected on the number of provincial monthly HFMD infections, air pollution, meteorological variables, and socioeconomic variables from to in mainland China.We used spatial autocorrelation to determine the aggregate distribution of HFMD incidence.Spatiotemporal patterns of HFMD were analyzed, risk maps were developed using the Bayesian spatiotemporal model, and the impact of potential influencing factors on HFMD was assessed.
Results: In our study, from to , the HFMD annual incidence rate in all provinces of mainland China ranged from . to .per , people, with an average annual incidence rate of . .The temporal risk of HFMD for Chinese provinces exhibited cyclical and seasonal characteristics.The southern and eastern provinces had the highest spatial relative risk (RR > ) from to .The HFMD incidence risk in provinces (Hunan, Hubei, and Chongqing) located in central China increased over time.Among the meteorological variables, except for the mean two-minute wind speed (RR .; % CI ., .
), all other variables were risk factors for HFMD.High GDP per capita (RR .
) was a protective factor against HFMD.The higher the birth rate was (RR .
), the higher the risk of HFMD.Health workers per , people (RR .; % CI .

Introduction
Hand, foot, and mouth disease (HFMD) is a widespread infectious illness primarily attributable to enterovirus infection (1).The majority of cases indicate that HFMD is a self-limiting disease.However, a small percentage of patients will have neurological complications and potentially death (2).For more than 20 years, HFMD has occurred frequently in many countries, such as Australia, Malaysia, and China (3).On average, 96 people per 100,000 are infected with HFMD each year during the period from 2008 to 2021 (4).Furthermore, studies have revealed that the average cost per case of HFMD ranges from $990 to $3,000, indicating that HFMD poses a serious economic burden as well (5).As of today, the specific antiviral therapeutic regimen for HFMD remains unavailable.Thus, public health measures related to prevention and control are still important for the management of HFMD.
There have been numerous prior studies on the spatiotemporal pattern of HFMD (6)(7)(8)(9)(10)(11). Nevertheless, the majority of the related researches have been carried out in just one province or local region with a small spatial size, such as Shaanxi Province (12), Xinjian (6), Henan Province (13), or the Ili River Valley Region (14).Analysis of the spatial distribution and temporal trend of HFMD incidence nationwide helps to determine high-risk provinces and periods, thus identifying provinces where prevention and control efforts are most needed and maximizing and rationalizing the utilization of limited public health resources.Wang et al. used data from mainland China from 2008 to 2012 to detect spatiotemporal clusters (15).In 2014, the Ministry of Health of China provided new guidance for HFMD prevention and control.Six years have passed since the last relevant policy was issued in 2008.In the first half of 2016, the first inactivated vaccine for EV-A71 was marketed in mainland China (16).The prevalence and spatiotemporal characteristics of national HFMD may be impacted by these interventions to some extent (17)(18)(19).However, studies on spatiotemporal patterns that contain fully recorded incidence data on nationwide HFMD cases from 2012 to the present are limited.There was only one study regarding spatiotemporal cluster detection but no risk estimation (20).Therefore, it is necessary to estimate the smoothed map of HFMD risk in recent years to obtain the spatiotemporal risk of HFMD across the country and to observe whether changes have occurred in high-risk provinces.
Identifying risk factors for HFMD is valuable in controlling outbreaks (21).The correlations between HFMD incidence and meteorological factors (22)(23)(24)(25), air pollution (26), and socioeconomic factors (27,28) have been supported by many previous studies.However, most studies on underlying HFMD explanatory variables have concentrated on a specific region or city, and few studies have focused on all of China.In addition, few studies have been conducted to assess the effect of air pollution on the risk of HFMD infection from a national perspective.The impacts of air pollution and meteorological and socioeconomic factors on HFMD have not been explored nationwide.Using Bayesian spatiotemporal models, spatiotemporal variations can be identified simultaneously, the impacts of possible explanatory variables can be assessed (29), and uncertainties due to residual confounding factors can be controlled.
The purpose of the present study was to map the spatiotemporal risk based on a Bayesian spatiotemporal model using province-level HFMD data from 2014 to 2017 in mainland China.The associations between air pollution, meteorological conditions, socioeconomic factors, and HFMD incidence were also analyzed.

Methods . Study area
With complex and diverse topographic conditions and vast land area, China is located in the eastern part of Asia.It has several different climate zones, from temperate to subtropical.The topography from the west and east is from high to low.In mainland China, there are 31 provincial administrative regions (Figure 1).

. Data sources
We used monthly province-level HFMD case data from January 2014 to December 2017 from China's Public Health Science Data Center (30) and National Health Commission (31).
The raw monthly meteorological data were first collected by us from the China Meteorological Data Service Center (32) for each province's meteorological monitoring stations.First, the location information of the prefecture-level cities in each province was incorporated, and the monthly average data of all the meteorological monitoring stations were used to obtain the monthly average meteorological data of 334 prefecture-level cities in mainland China via the ordinary kriging interpolation method.Subsequently, the annual population of prefecturelevel cities within each province was used as a weight to derive the weighted monthly average meteorological data for the corresponding provinces.For the four municipalities, the data from those monitoring sites were averaged to calculate the monthly mean concentrations of meteorological factors.The meteorological variables included the monthly mean temperature ( • C), monthly mean rainfall (10 mm), monthly mean relative humidity (%), monthly mean air pressure (hPa), monthly mean hours of sunshine (h), and monthly mean two-minute wind speed (m/s).
Daily air quality indices (AQI) for each prefecture-level city for the same period were collected from the website (33), and monthly statistics for each prefecture-level city were obtained after averaging.Resident population data for each prefecture-level city for each year were incorporated as weights, and a weighted average was used to obtain monthly average AQI statistics for each province (municipality).
Provincial socioeconomic indicators, which were selected from the aspects of the economy, demographic factors, population mobility, and health resources, were obtained from the National Bureau of Statistics (34) and included per capita gross domestic product (GDP pc , $ 140.60), health workers per thousand people, hospital beds per 10,000 people, the number of village clinics (10,000), total passenger traffic on public transport (100 million), and birth rate (‰).

. Statistical analysis
In this study, spatial autocorrelation and Moran's I were analyzed.Spatial autocorrelation analysis is composed of both global and local indicators.Global spatial autocorrelation analysis mainly studies the spatial distribution of attribute values of regional spatial objects.Whether the distribution of attribute values between two adjacent spatial areas is correlated or random can be obtained through local indicators of spatial association (LISA).
Moran I is computed in Equation 1: where y i or y j represents the number of HFMD cases in province i or j.The average number of cases in 31 provinces is y, and ω ij denotes the spatial proximity matrix of two provinces, LISA is computed in Equation 2: I i = 0 indicates that the distribution of attribute values between adjacent spatial units is random.I i > 0 indicates that the distribution of attribute values of the spatial object and the neighboring spatial objects has a positive correlation ("high-high" or "low-low"), and I i < 0 indicates that the distribution has a negative correlation ("high-low" or "low-high").
Bivariate spatial autocorrelation, usually represented by the bivariate Moran I statistic, characterizes the correlation between one variable and the spatial lag of another.We used the bivariate local indicators of spatial autocorrelation (BiLISA) to capture the relationship between the value of one variable in one spatial unit and the average of another variable in neighboring areas.We drew cluster maps of BiLISA with P < 0.05 between meteorological variables, AQI, socioeconomic variables and HFMD incidence.
The BiLISA is computed in Equation 3: where x i represents the explanatory variables in province i.We constructed two Bayesian spatiotemporal models, investigated the fluctuations in HFMD incidence in 31 provinces in mainland China in spatial and temporal dimensions and analyzed the influencing factors.The monthly number of HFMD cases follows a negative binomial distribution, denoted as: where E it , θ it and y it represent i (i = 1, 2, . . .31) provinces t (t = 1, 2, . . .48) months of the population, incidence, and the number of cases, respectively.
which is modeled as follows: Equation 5 is the model used to estimate the spatial and temporal effects.The intercept β 0 is equal to the logarithm of the average HFMD incidence in all inland provinces; u i is the spatially structured random effect, which reflects the spatial autocorrelation among spatial area units; v i is the spatial heterogeneity; γ t is the time structure random effect, taking month as the analysis scale unit, reflecting the overall trend of change in the time effect over the study period; and ϕ t reflects the temporal heterogeneity.δ it indicates the interaction effect between space and time.
Equation 6 is the model with potential influencing factors added to Equation 4, where β k represents the model coefficient and p represents the number of variables.
β 0 and β k adopted uninformative prior distributions.The prior distribution of u i is intrinsic conditional autoregression (ICAR) (29), and v i ∼ N(0, σ 2 v ).The reciprocal of σ 2 u and σ 2 v as unknown variance parameters are changed into precision parameters τ u and τ v , respectively, and both of their prior distributions are gamma (1, 0.0005).γ t follows a second-order autoregressive prior distribution (29).The unstructured time effect follows a normal distribution with a zero mean, ϕ t ∼ N(0, σ 2 ϕ ).We assumed the spatiotemporal interaction term to follow a normal distribution with a zero mean; for example, δ it ∼ N(0, τ δ ).
The preliminary Spearman correlation analysis indicated that rainfall was relatively strongly correlated with the mean temperature (0.78, P < 0.001) and mean relative humidity (0.75, P < 0.001).Therefore, the initial variable of mean rainfall was removed.Next, the multi-collinearity among other variables was estimated by calculating the variance inflation factor (VIF) indicator.Normally, explanatory variables (VIF > 10) were not included in the Bayesian model due to their severe collinearity.We chose the deviance information criterion (DIC) metric to evaluate

. Epidemiological characteristics of HFMD
Table 1 demonstrates the descriptive statistics of the yearly HFMD cases and their incidence in this study.The annual incidence rates varied between 138.80 and 203.15 per 100,000 people.In mainland China, 2014 represented the highest annual incidence of HFMD during the study period, while the lowest number of cases occurred in 2017.
The spatial distribution of the average annual incidence in each province is shown in Figure 2. The mean incidence of HFMD differed across provinces and varied widely.Guangxi Province had the highest incidence, with an average of 563.01 per 100,000 individuals, followed by Hainan Province (491.97/100,000),Guangdong Province (361.27/100,000), and Fujian Province (242.13/100,000).The incidence of HFMD in the eastern (Jiangsu, Anhui, Fujian, Shanghai, and Shandong Provinces) and southern coastlands (Guangdong, Guangxi, and Hunan Provinces) was higher than that in northwestern (Xinjiang, Qinghai, and Gansu Provinces) and northeastern (Heilongjiang, Jilin, Liaoning, etc.) regions of China.
The monthly incidence of HFMD in all provinces generally exhibited seasonal characteristics (Figure 3).As shown in Figures 3A, B, the peak months of HFMD incidence varied in different provinces each year.Provinces with two epidemic peaks each year were located mainly in southern (Guangdong, Guangxi, and Hainan provinces), central (Hunan, Chongqing, and Hubei provinces), and eastern (Jiangsu, Anhui, Fujian, Shanghai, and Shandong provinces) China, with a high peak in May-June and a low peak in September-October.However, as shown in Figure 3C, northern China (e.g., Heilongjiang, Jilin, and Liaoning) had an epidemic peak only in June-July. .

Spatial autocorrelation analysis
The global spatial autocorrelation analysis indicated a positive spatial autocorrelation of the provincial HFMD incidence from 2014 to 2017 (Table 2); the Moran I value varied between 0.24 and 0.30, and all the values were statistically significant (P < 0.01).
Figure 4 shows the results of the LISA for each province across mainland China from 2014 to 2017.During the study period, Guangxi, Hainan, and Hunan provinces, located in central and southern China, were the main high-high (HH) regions, and only Hainan Province showed a strong aggregated distribution of HFMD outbreaks in 2015.However, Hunan Province was just the opposite and was the HH province except in 2015.Guizhou and Jiangxi, the two provinces with low-high (LH) clusters, were adjacent to Hunan.Only Inner Mongolia was a low-low (LL) cluster in 2016.

. Spatial autocorrelation between variables
The results of the bivariate local spatial autocorrelation (BiLISA) analysis of the explanatory variables are displayed in Figure 5. Provinces with statistically significant BiLISA values are labeled in different colors, according to the findings, while all other provinces are labeled in gray.Except for the average 2-minute wind speed, we found that the LL and high-low (HL) provinces were situated mostly in Inner Mongolia and northeast and northwest China based on cluster maps of the BiLISA of the relationships between meteorological variables, air pollution, socioeconomic variables, and the incidence of HFMD.This suggests that regardless of whether the meteorological, air pollution and socioeconomic variables of the provinces (e.g., Xinjiang, Inner Mongolia and Jilin) located in these regions were at higher or lower levels, the incidence of HFMD in their neighboring provinces was at lower levels.For the monthly average meteorological variables (except for sunshine hours), the HH provinces were located mainly in central and southern China, where the temperature, rainfall, relative humidity, air pressure, and 2-minute winds in Hunan, Guangdong, and Jiangxi and the incidence of HFMD in neighboring provinces were at high levels.For air pollution and socioeconomic variables (except for health workers per 1,000 people), the HH regions varied by exploratory variables but were still located in south-central China.Especially for Hunan Province, the meteorological variables, air pollution and socioeconomic variables in this region and the incidence of HFMD in neighboring provinces were at high levels. .

Spatiotemporal pattern
The spatial relative risk (RR) values for HFMD in each of the 31 provinces during the study period are shown in Figure 6.The risk maps showed vast differences in spatial RR across provinces.The spatial RR (>3.0) was consistently highest in the three southern provinces (Guangdong, Guangxi, and Hainan) shaded in deep red in Figure 6.The spatial RR of Hunan Province, located in

Discussion
Obtaining HFMD risk across the Chinese mainland is of great value for the overall allocation and optimization of public health resources.However, studies of national-scale spatiotemporal analyses of HFMD data in China are limited.In this study, we focused on characterizing the global spatiotemporal pattern of HFMD risk and producing a provincial risk map for the Chinese mainland from 2014 to 2017.Our study is the first to explore the association between AQI and HFMD on a national scale.During the study period, southern Chinese provinces had a stable high risk of HFMD, while the risk in eastern provinces fluctuated.From 2014 to 2017, the central provinces (Hunan, Hubei, and Chongqing) gradually became high-risk areas.
Meteorological and socioeconomic variables partially explain this spatiotemporal pattern.
The temporal relative risk of HFMD from 2014 to 2017 indicated the cyclical and seasonal characteristics of the disease in mainland China.HFMD had an epidemic cycle of 2 years.There was evidence that the decline in population immunity between adjacent outbreaks and the increase in the number of newborns led to this cyclical pattern (35).The epidemic of HFMD in mainland China was characterized by two seasonal peaks at the beginning of summer and autumn, with the former being strong and the latter weak, which was consistent with the findings previous research (2).The seasonal prevalence in southern and northern mainland China, however, was inconsistent.We found that only one epidemic peak in northern China occurred in June-July.The onset of HFMD was affected by latitude to some extent according to previous studies (36); for example, tropical and temperate regions located at low latitudes had a higher incidence of HFMD.The duration of the epidemic periodicity of HFMD varies with latitude, with southern regions located at low latitudes showing the strongest semiannual periodicity, but northern regions at high latitudes showing a 1-year periodicity (1).Furthermore, meteorological factors and socioeconomic conditions in northern and southern China showed great discrepancies due to different geographical topography and development, and regional studies also showed different epidemiological patterns of HFMD (37).This explains the discrepancy in the seasonal characteristics of HFMD between northern and southern China.
The spatial RR of provinces located in the eastern coastlands remained at a high level (ranging from 1.5 to 10) during the study period in the spatial dimension.However, from 2014 to 2017, high-risk provinces in China were gradually concentrated in the southern and central parts of the country.In our study, three coastal provinces in southern China, namely Guangdong, Guangxi, and Hainan, consistently exhibited the highest risk (>3.0) of HFMD.The predominant climate types in the eastern and southern coastal provinces are subtropical and tropical monsoon climates.These provinces, which have adequate rainfall, heat, and humidity, are more likely to experience HFMD outbreaks (38)(39)(40).The high population density and mobility of people in coastal provinces, with many opportunities for interpersonal communication and promote the transmission of the virus, resulting in a higher risk of HFMD spread (37,41,42).Moreover, the spatial RR map revealed that the central provinces of China have gradually become high-risk regions for HFMD.Driven by the country's macroeconomic policies, since the time that the central China rise strategy was proposed, plenty of supportive policies have been given to urbanization, infrastructure development and transportation network improvement in the central regions, which has led to high levels of migration (43).Furthermore, the central provinces have abundant heat, sunshine duration and high humidity because of short, mild winters and long summers.A suitable natural environment and socioeconomic conditions contribute to the spread and reproduction of the virus in the central region (44,45).These findings indicated that the epidemic of HFMD in the central provinces is becoming more serious and that prevention and control efforts should be strengthened.
This study suggested that temperature, air pressure, relative humidity, and sunlight duration all increase the risk of HFMD, which was consistent with previous studies (46, 47).Hot and humid environments contributed to the breeding of enteroviruses (48,49).Higher wind speed and HFMD risk had a negative correlation.Higher wind speeds may accelerate evaporation, making it difficult for the virus to multiply and survive in dry food and environments (50,51).The infection rate of certain airborne diseases decreased as ventilation rates in the environment increased according to the  results of a systematic review (52).Hence, high ventilation resulting from high wind speed reduces the opportunities for HFMD.Air pollution (AQI) was not related to the risk of contracting HFMD.The risk of HFMD decreased with better economic conditions (GDP pc , RR: 0.9889), which was consistent with previous studies (53).The higher the per capita GDP of the region is, the larger the budget for disease prevention and control.The infected individuals can be treated in time, and the less likely it is that other uninfected children will contract HFMD.Birth rate (RR = 1.2010) was a risk factor for HFMD incidence.In particular, health workers per 1,000 people (RR 1.2010; 95% CI 1.0443, 1.3771) increased the risk of HFMD.Areas with more healthcare workers tend to have more healthcare resources and better healthcare systems.Medical facilities can report disease immediately to higher health authorities, resulting in a tendency toward higher risk of morbidity (54).In this study, bivariate local spatial autocorrelation was utilized to investigate the relationships between the explanatory variables of each province and the HFMD incidence in neighboring provinces.The analysis showed a polarization trend, with LL and HH provinces located in northwestern and southern China, respectively.The burden of HFMD incidence in southern China (Guangxi, Hainan) as well as in provinces along the eastern coast of China (Fujian, Zhejiang) and the monthly meteorological variables (including temperature, rainfall, relative humidity, barometric pressure, and 2-minute wind speed) in south-central China were high levels, which can be attributed to the fact that these areas were all in the same climatic zone of China and had similar climatic conditions.The socioeconomic variables of Guangdong, Jiangxi, and Hunan and the incidence of HFMD in surrounding provinces were at high levels, which may partially explain the emergence and persistence of provinces at high risk of HFMD in and around southcentral China.Moreover, Hunan Province, a HH region, provides further evidence of the gradual transformation of the central region into a high-risk area for HFMD.
Several limitations existed in our study.First, the meteorological data we used may not be consistent with the climatic conditions of HFMD patients at the time of infection, as these data were obtained from meteorological monitoring stations in each province and could not be specific to individuals.Secondly, some mild cases of HFMD were underreported because infected individuals did not always go to the hospital, resulting in surveillance systems not recording the actual number of HFMD patients.

Conclusion
Within the present study, we analyzed spatiotemporal patterns using national provincial HFMD data from 2014 to 2017, estimated spatial smoothing maps of risk, and explored potential influencing factors.The risk of HFMD was characterized by cyclical and seasonal patterns in mainland China.The central provinces (Hunan, Hubei, and Chongqing) gradually became high-risk areas for HFMD during 2014-2017.Prevention and control efforts should be strengthened in the central region of China.

FIGURE
FIGUREProvincial-level administrative regions in mainland China.
i and j. n is the number of provinces.The larger the absolute Moran I value is, which takes values from −1 to 1, the stronger the correlation.Positive values (P < 0.05) indicate positive correlations, and negative values (P < 0.05) indicate negative correlations.

FIGURE
FIGURESpatial distribution of the average incidence of HFMD by province in mainland China, -.

FIGURE
FIGURE Monthly incidence of HFMD in mainland China from to .(A) Monthly incidence of HFMD in provinces located on the southern and eastern coasts of mainland China.(B) Monthly incidence of HFMD in provinces located in central and south-western mainland China.(C) Monthly incidence of HFMD in provinces located in western and northern mainland China.

FIGURE
FIGURESpatial association cluster map of HFMD incidence in mainland China, -.

FIGURE
FIGURESpatial distribution of the BiLISA between provincial HFMD incidence and influencing factors.

FIGURE
FIGUREThe spatial RR map of HFMD incidence in mainland China from to .

FIGURE
FIGURETemporal risk trend of HFMD incidence from to in mainland China.
TABLE Annual cases and incidence of HFMD in mainland China from to .s goodness-of-fit and select the best Bayesian model.For the same data, the smaller the model DIC value is, the better.A P value < 0.05 indicated that the conclusion was statistically significant.The spatial autocorrelation analyses were conducted using the "spdep" package (R 4.1.2) and GeoDa software version 1.16.0.16.We used ESRI ArcGIS software version 10.8 (Environmental Systems Research Institute, USA) to map the spatial distribution of the data.The Bayesian spatiotemporal models were conducted in R 4.1.2with the "INLA" package.
TABLE Global spatial autocorrelation analysis and significance test results.
TABLE Posterior distribution of intercept terms and influencing factors.