Supplementary Figure S1 shows the geographical location of the study area—Heilongjiang, Jilin, and Liaoning provinces. The three provinces are located in the northeastern of China and have medium levels of economic development and population size.

We obtained HFRS case surveillance data from the National Public Health Data Center of China ¹ for the study area covering the period from 2005 to 2018. All patients were diagnosed according to the HFRS management criteria issued by the Ministry of Health of the People’s Republic of China. We obtained the corresponding daily weather data, including air temperature and dew point temperature, from the China Meteorological Data Sharing Service (data.cma.cn). Pollutant information, including CO, NO₂, and O₃, was originally sourced from the National Oceanic and Atmospheric Administration (NOAA).

We used the R language implementation of the moving epidemic method (MEM) (package “mem”), which is available online for free. The method is based on a complex mathematical algorithm that can be summarized in three steps. The first step is the division of the pre-epidemic, epidemic, and post-epidemic periods. In the second step, the pre- and post-epidemic values of the historical seasons are used to calculate the baseline and epidemic thresholds. In the third step, the maximum values of n surveillance indicators during the epidemic period are selected separately to calculate different epidemic intensity thresholds. The unilateral 50%CI upper limit of the geometric mean of the n maximum surveillance indicators during the epidemic period is defined as the medium intensity threshold, the unilateral 90%CI upper limit as the high-intensity threshold, and the unilateral 95%CI upper limit as the very high-intensity threshold.

Distributed lag non-linear models (DLNM) have been widely used to assess the exposure–lag–response relationship between environmental factors and human diseases such as congenital heart disease, hand, foot, and mouth disease, and chronic sinusitis (8, 10–12). The model can be written as follows:

To analyze the lag and extreme effects of climate factors, air temperature, dew point temperature, wind direction, and wind speed were considered and applied to the cross-basis functions of a DLNM. Here, Y_t is the number of the HFRS cases in monthly t; α₁ is the intercept of the entire model; NS is a natural cubic spline that acts as a smooth function of the model; M represents the estimated climate or pollutants variable related to HFRS; and X_t represents other climate and pollutant variables involved in the pathogenesis of HFRS, for which non-linear confounding effects are adjusted. When constructing the meteorological factor model, ∑ ( X t ) does not exist, whereas in the pollution model, meteorological factors are used as confounding factors to construct ∑ ( X t ) . The NS was applied to adjust for the monthly confounding effects in the model. Month is a binary variable used to control the effect of time, and β represents regression coefficients. The optimal degrees of freedom (df) for the spline function were estimated using the Akaike information criterion for quasi-Poisson (Q-AIC) and minimum partial regression coefficient (PACF_min) criteria. The NS with 4 df was used for the climate factors, except for wind direction, which used 5 df during the period of low epidemic intensity. For both the high epidemic intensity period and the overall model, the NS with 4 df was applied to the climate and pollutant factors. The lag space was set to 3 df. The NS with 2–3 df/year was applied to the time variable in both pollutant and climate models. The climate model was constructed using the glm () function, while the pollution model was constructed using the gam () function.

Subsequently, a generalized additive model (GAM) was used to explore the interaction between the pollutants and the prevalence of HFRS. The model formula can be written as follows:

α₂ is the intercept; X₁ represents the AQI, whereas X₂ and X₃ denote the other two pollutants. s () indicates a penalized spline function. s (X₁, X₂) represents the spline function for the interaction between the parameters X₁ and X₂. X₁, X₂, and X₃ represent 6-month lagged variables. ∑ ( X t ) represents the factors of climate.

A Gaussian process (GP) can be regarded as an extended function of a multivariate Gaussian distribution, which can be applied to a wide range of variables. In a Gaussian process (GP), it is assumed that any finite set of data follows a multivariate Gaussian distribution. Prior beliefs concerning the relationships between variables are incorporated into these (an infinite number of) multivariate Gaussian distributions to create a model that represents the observational variance. The combination of multiple Gaussian distributions in a GP can effectively model non-linear relationships and is more versatile than traditional parametric models, which depend on fitting a global model. This is because multivariate Gaussians can represent local covariance patterns between individual sites (13).

Support vector machines (SVMs) are a non-probabilistic binary linear regression method. Given a set of training data labeled as belonging to one of two classes, the algorithm maps the data into a space and defines a hyperplane that maximizes the margin between the two classes to separate them. This plane is called the “maximal marginal hyperplane.” An algorithm uses a kernel approach to acquire non-linear mapping to the feature space if linear integration is impossible. Thus, the hyperplane of the feature space stands for the non-linear boundary of the determination in the input space (14). All model metrics are compared using traditional machine learning metrics such as RMSE, R², and MAE (15–17). A total of 75% of the dataset is used as the training set, while the remaining 25% is used as the test set. All analyses in our study were performed using R software (version 4.1.3).

A total of 59,431 HFRS cases were reported in the three eastern provinces of China from 2005 to 2018, showing a decreasing trend each year (Table 1). This was followed by the main epidemic area in Heilongjiang province, with a total of 28,074 cases until 2018. The incidence of influenza was primarily observed in the individuals aged 15–39 and 40–59 years, accounting for 86.42% of all cases.

Based on Table 1, however, there was a short-term rise in the cases from 2012 to 2014. We also performed a calculation of the prevalence threshold and determined from Supplementary Table S1 that the optimal parameter δ was 7.0 after the calculation of the popular threshold model. As shown in Table 2, the years with a high prevalence threshold were 2005–2007 and 2012–2014, while the years with a low prevalence threshold were 2008–2011 and 2015–2018. Based on the threshold model prediction shown in Table 2 and Figure 1, it was concluded that the epidemic months were primarily concentrated in November.

The summary statistics for all HFRS cases and environmental variables in northeastern China are shown in Supplementary Table S2. The Spearman’s rank correlation coefficient analysis showed that HFRS was significantly correlated with air temperature (r = −0.18, p < 0.05), dew point temperature (r = −0.23, p < 0.01), wind direction (r = 0.22, p < 0.01), and wind speed (r = 0.29, p < 0.01) (Supplementary Table S3). As shown in Supplementary Figure S2, these climate factors were associated with high relative risk at the lags above moderate levels, except for air temperature.

From the dose–response relationship shown in Supplementary Figure S3, air temperature showed mostly a U-shaped relationship with the risk of HFRS, both in general and across the different regions and age groups, while the other factors mostly showed an arch bridge-shaped relationship. In Liaoning province, air temperature, dew point temperature, and wind speed all showed a parabolic decreasing trend in their relationship with HFRS risk. As shown in Supplementary Table S5, the climate lag effect was weak during the low prevalence threshold period, with sensitivity mainly concentrated in the high prevalence areas of Heilongjiang province and the 0–14 years age group, where the main lag factor was low wind direction. As shown in Supplementary Table S6, the meteorological lag effect was higher during the high prevalence threshold period, with sensitivity mainly concentrated in the 0–14 years and 60 years and above age groups, where the main lag factors were cold air and hot dew point. When comparing the climatic lags during the low and high prevalence threshold periods (Supplementary Tables S5, S6), we found that low wind direction and windy conditions showed a dangerous lag effect on HFRS onset (OR > 0), while high wind direction and windless conditions showed a protective lag effect (OR < 0). In addition, air temperature showed protective effects at both low and high levels, while cold air showed a dangerous effect in the 0–14 years age group during the high prevalence threshold period (OR (95% CI): 3.2e+17(8.4e+08, 1.2e+26)). Cold dew point had a little lag effect, while hot dew point showed a protective effect during the low prevalence threshold period. However, this effect was reversed during the high prevalence period.

The Spearman’s rank correlation coefficient analysis showed that HFRS was significantly correlated with the AQI (r = 0.40, p < 0.05), PM_2.5 (r = 0.37, p < 0.05), and PM₁₀ (r = 0.40, p < 0.01) (Supplementary Table S4). As shown in Supplementary Figure S4, these factors were associated with high relative risk at the lags above high levels, except for PM₁₀. From the dose–response relationship shown in Figure 2, PM_2.5 mostly showed an arch bridge-shaped relationship, while the AQI and PM₁₀ mostly showed a U-shaped relationship with the risk of HFRS, both in general and across the different regions and age groups. In Jilin province and the 0–14 years age group, the AQI exhibited a parabolic decreasing trend, while PM_2.5 showed a parabolic increasing trend. As shown in Figure 3, in terms of the total pollution lags, the effects of the low-level pollutants were mainly concentrated in the long-term lag conditions (3–6 months), while the effects of the high-level pollutants were mainly concentrated in the short-term lag conditions (1–2 months). In terms of the lagging trend, PM_2.5 differed from the other pollution factors. As shown in Table 3, except for high-level PM₁₀, the lag effect of the other pollution factors was more pronounced, and the sensitivity was mainly concentrated in Liaoning province and the age group of 40–59 years. Among these, we found that low levels of the AQI and PM₁₀ and high levels of PM_2.5 showed a dangerous lag effect on the onset of HFRS (OR > 0), while extremely high levels of PM_2.5 (P95) showed a protective effect. In addition, high levels of the AQI and PM₁₀ and low levels of PM_2.5 showed a protective lag effect (OR < 0). However, at extremely high levels of the AQI (P95), a dangerous effect was observed.

From Supplementary Figure S5, we can see that the AQI interacted with PM_2.5 and PM₁₀ in relation to HFRS incidence. PM₁₀ was weakly positively correlated with the risk of HFRS, while PM_2.5 showed the opposite relationship. From the interaction effect shown in Figure 4, we found that low AQI combined with high levels of PM_2.5 and PM₁₀ had the greatest impact on HFRS onset. The results from the test in Supplementary Table S7 indicate that the model involving the interaction between PM_2.5 and the AQI performed better (R² = 44.1%). From Supplementary Table S8 and Table 4, the model fit was best in Liaoning province among the different regions (R² > 70%) and in the 15–39 age group. In addition, the GPR model showed the same fit as that of the SVM model. In the GPR model, the prediction results were good, except for the polydot kernel function. In the SVM model, good prediction results were observed with the radial and sigmoid kernel functions. Based on the SVM-radial model for exploring the importance of the variables related to HFRS, the priority order was the pollutant factors (in the order of the AQI, PM₁₀, and PM_2.5), followed by the climatic factors (in the order of windspeed, dew point temperature, and air temperature).