Introduction

The COVID-19 pandemic has imposed significant health and economic burdens all over the world1. A better understanding of the factors affecting the COVID-19 epidemic is critical to the design of tailored public health and social measures (PHSMs), e.g., travel restrictions, school closures, cancellation of public events and gatherings, etc. and much attention has been given to the impact of meteorological factors on the COVID-19 transmissibility and severity.

Over the last couple of decades, essential factors related to the transmission of viral respiratory diseases have been investigated such as the highly predictable seasonal pattern of influenza epidemics2. These epidemiological studies are supported by laboratory evidence that low temperature and/or humidity improve the stability of influenza virus3, impair the human innate immune system4 and contribute to the aerosol evaporation5,6.

Since the COVID-19 pandemic hit, many research groups worldwide have aimed to reveal the relationships between temperature and COVID-19 transmission. Previous studies examined the hypothesis that high temperature, humidity, wind speed, and ultraviolet (UV) radiation might have a reduction in transmission7,8,9,10,11,12,13,14,15. For example, in recent studies conducted worldwide, Jie et al. found that temperatures below 21 °C, relative humidity, and wind speed were negatively correlated with the number of daily cases for one year in 188 countries10. Simiao et al. investigated that 1° increase in absolute latitude is associated with a 4.3% increase in COVID-19 cases per million inhabitants which is consistent with the hypothesis that high temperature and UV radiation can contribute to the reduction in transmission11. Some of these investigated the possibility that the transmissibility is associated with temperature, where the transmissibility is often translated into the number of positive cases7,10,11,14; however, these studies did not fully account for the transmission dynamics influenced by PHSMs of various intensities. In this context, Yiqun et al. investigated the association of increased effective reproduction number and lower temperature (within the 20–40 °C range), lower humidity, and lower UV radiation15. Yet there is still divergence in the literature as another study indicated higher temperatures are not significantly associated with a reduction in total cases or effective reproduction number of COVID-19 in Canada16.

A few of the earlier studies have explored the association between temperature and mortality7,17,18 as an indicator of the clinical severity. Previous studies showed a 1% increase in temperature was associated with a 1.19% decrease in daily new COVID-19 deaths in 166 countries17 and about 6% lower mortality in the subsequent 30 days from the first death in the OECD countries and US states19. Non-linear relationships of temperature and mortality have also been investigated in the time-series analysis as such the daily cumulative relative death risk decreased by 12.3% for every 1.0 °C increase in temperature20. However, day-to-day fluctuations in the number of deaths are also vulnerable to the epidemic dynamics. Transmission dynamics of infectious diseases should also be considered when performing the regression models because observation of each case with a contagious disease is not independent, which characteristics is referred to as dependent happening and explicitly distinguishable from other non-communicable diseases; otherwise, such inferences get largely biased21,22.

The present study explored the association between temperature and both the transmissibility and the severity of COVID-19 from early 2020 to early 2021. We used the effective reproduction number \((R_{t} )\), defined as the mean number of secondary cases generated by a single primary case, to quantify the transmissibility of the ongoing epidemic in Tokyo. To explore the association between temperature and severity, we used case fatality risk (CFR), an epidemiological measurement of severity. Crude CFR calculated from the ratio of the cumulative number of deceased cases to the cumulative number of confirmed cases can underestimate the actual CFR when cases are increasing and overestimate it when they are decreasing due to the time that passes from the onset illness to death. Such issues are also known as right censoring. We, therefore, estimated the time-delay adjusted CFR for every illness onset date, which accounted for the delay.

Results

The daily mean temperature from 15th February 2020 to 28th February 2021 in Tokyo is shown in Supplementary Fig. S1. Supplementary Fig. S2 (a) and (b) show the epidemic curves for confirmed cases by two age groups (under 70 s and over 70 s) and deaths, respectively. The epidemic curve and estimated median value of \(R_{t}\) with 90% credible intervals (CrI) from 15th February 2020 to 28th February 2021 are shown in Fig. 1. Analysing the impact of temperature on \(R_{t}\), the overall cumulative exposure–response relationship of temperature on \(R_{t}\) was non-linear, with lower temperature leading to higher RR (Fig. 2A). The RR corresponding to temperature at the first percentile (3.3 °C) was 1.3 (95% confidence interval (CI): 1.1–1.7). Figure 2B shows the three-dimensional plot of RR with temperature and lags up to 7 days. We found that the greatest risk of cold effects occurs in the day of exposure, increasing in 3–7 days of exposure.

Figure 1
figure 1

Transmission dynamics from 15th February 2020 to 28th February 2021 in Tokyo, Japan. Blue line represents median, blue shading represents 95% credible intervals of the estimated effective reproduction number from 15th February 2020 to 28th February 2021. Green bars show the observed number of COVID-19 cases with confirmed dates in Tokyo.

Full size image
Figure 2
figure 2

Overall and three-dimensional plots of relative risks with the reference at 15.3 °C. (A) The three-dimensional plot of the association between daily mean temperature (°C) and the effective reproduction number over the lags of 7 days. The reference value of temperature was median temperature (15.3 °C). (B) The estimated overall effects of mean temperature (°C) over 7 days on \(R_{t}\). Blue line shows the mean relative risks, and 95% confidence intervals are shown in the gray shadings.

Full size image

Figure 3 shows the temporal variation of time-delay adjusted CFR and unbiased CFR from 25th May 2020 to 28th February 2021. As of 28th February 2021, the time-delay adjusted daily CFR and the unbiased CFR were 8.21% (95% CI: 4.50–12.9) and 2.42% (95% CrI: 2.41–2.43) (Supplementary Fig. 4), respectively. Figure 3 illustrates the temporal deviations from the baseline value of CFR, i.e., the unbiased CFR. To examine the potential for the temperature to contribute to changes in CFR, we estimated the overall effect of temperature with the reference of 18.6 °C (Fig. 4A). The harmful effect was seen to increase as temperature increased from the reference, and moderately cold temperatures were associated with high RRs of CFR. The three-dimensional plot of RR with temperature and lags for CFR is displayed in Fig. 4B, cold temperatures have an obvious impact on the day of exposure (lag day 0, the illness onset day) and we found a week delayed effect on both high and cold temperatures. The outcomes of the possible confounders are shown in Supplementary Table S4.

Figure 3
figure 3

Temporal variation of time-delay adjusted case fatality risks (CFR) with unbiased CFR from 25th May 2020 to 28th February 2021 in Tokyo, Japan. The mean values of time-delay adjusted daily case fatality risks (CFR) from 25th May 2020 to 28th February 2021 are shown with a purple line. The shade region represents the 95% confidence intervals. The blue dot line shows the unbiased case fatality risk as 2.42% (95% credible interval: 2.41–2.43). If the time-delay adjusted daily CFR gets higher or lower, it is caused by random noises or other variables which have causal relationships. The unbiased case fatality risk plays a key role as a reference of the daily CFR.

Full size image
Figure 4
figure 4

Overall and three-dimensional plots of relative risks with the reference at 18.6 °C. (A) The three-dimensional plot of the association between daily mean temperature (°C) and time-delay adjusted case fatality risks (CFR) over the lags of 14 days. The reference value of temperature was median temperature (18.6 °C). (B) The estimated overall effects of mean temperature (°C) over 14 days on CFR. Blue line shows the mean relative risks, and 95% confidence intervals are shown in the gray shadings.

Full size image

Similar results were obtained in sensitivity analysis under different lags and adjustment of several meteorological variables (Supplementary Figs. S6–S11) for the transmissibility and severity analysis. We assessed the impact on the overall effects and delayed effect and successfully checked the robustness of the primary analysis.

Discussion

The present study was the first to comprehensively quantify the association between temperature and the epidemiological dynamics of COVID-19 in Tokyo using the effective reproduction number and time-delay adjusted daily CFR considering the lagged effect of temperature. Though the epidemiology of COVID-19 is differentiated by the substantial transmissibility and severity which are measured by reproduction number and CFR, there is no study to explore the contribution of temperature using both of the two rigorous epidemiological measurements appropriately, specifically the association between temperature and CFR, to our best knowledge.

\(R_{t}\) rose explicitly at low temperatures; for example, RR of 3.3 °C (the 1st percentile of temperatures, defined as extremely cold temperature) was estimated as 1.3 (95% CI: 1.1–1.5) in the 0–1 days lag from an infected date (Supplementary Table S1), with median of all the temperatures (15.3 °C) as reference temperature. This indicates that the cold effects appear in short lags and the overall effect is more plausible in low temperature (Fig. 2A). Our results are consistent with most of the previous works in that COVID-19 incidence decreases as temperature increased and non-linear associations of temperature and COVID-19 transmissibility were observed15,23,24. In winter, human behaviors such as less ventilations in the rooms and close contact might contribute to high transmissibility25. In a recent study, relative impacts of meteorological factors compared to PHSMs were reported to be small in the early stage of the local epidemics as 2.4% and 2.0% of the variation in effective reproduction number are attributable to temperature and absolute humidity while 13.8% are explained by government response in 409 cities across 26 countries23. Another study in the USA investigated attributable fractions of temperature, specific humidity, and UV as 3.7%, 9.4%, and 4.4%, respectively and observed the differences of each factor depending on seasons15. Impacts of meteorological factors on both transmission and severity might depend on climate zones and seasons and further investigations in different climate conditions are helpful to understand the contributions of meteorological factors in accordance with varying phases of pandemic as an introduction of vaccination and distribution of variants.

The exposure–response relationships of population mobilities and meteorological factors with \(R_{t}\) (Supplementary Table S2) are consistent with previous work12,13,26. For example, the residual and workplace mobility changes were not slightly related to the fluctuation of \(R_{t}\) while the recreation mobility change was significant. This relationship in Tokyo were reported in the previous study26. Analysis including the other possible meteorological factors, i.e., solar radiation and wind speed, have a slight difference compared to the main analysis (Fig. 2) as the shape of the exposure–response outcome shows (Supplementary Fig. S9).

As for the severity, we found that low temperatures had a strong association with high CFR in the short lag periods. For example, RR of 2.3 °C (the 1st percentile of temperatures, defined as extremely cold temperature) and 5.8 °C (the 10th percentile of temperatures, defined as extremely cold temperature) were 2.0 (95% CI: 1.2–3.5) and 2.8 (95% CI: 1.7–4.7), respectively, in 0–2 days lag periods from illness onset dates (Supplementary Table S3), with a median of all the temperatures (18.6 °C) as reference temperature. While cold effects appeared in the short lag period and showed a slight gradual decline, extremely high temperatures were associated with higher CFR from few days after the illness onset, and those effects were stably maintained for two weeks (Supplementary Table S3).

Plausible mechanisms explaining the association between temperature and high CFR of COVID-19 remain undetermined. Even though many studies have postulated seasonal variations and the impact of temperature in the transmissibility of infectious diseases, little is known about the association of temperature and severity of the contagious disease. Since the common infectious respiratory diseases such as the influenza virus, circulate in the cold season, the impact of high temperature on severity is yet to be explored. One study indicated that immune response against the common cold virus can be impaired under the environment with cold temperature27. As a previous study on the impact of heat effect showed, extreme high temperature dampens physiological responses when the body temperature exceeds its normal range28. These phenomena might have contributed to the result observed in the current study. For high CFR of moderate temperature, one of the possible explanations is the difference in human movement which can affect the exposure level with ambient temperature; however, further investigation is needed.

There are several limitations to be noted. First, the effect of air particulates on increase in \(R_{t}\) and CFR were not considered in our study due to the data availability although some previous studies explored the relationship between air particulates and COVID-19 transmissibility and severity in other areas of the world29,30. Second, the assumption for epidemiological time-delays (e.g., generation time and incubation period, etc.) was imposed not to be contracted by meteorological factors or interventions. Third, we did not consider the uncertainty of variables in the regression models, which were derived from the mathematical models. Forth, age-specific CFR was not estimated in the present study due to scarce data. However, we attempted to make the best use of the data with back-projected incidence on illness onset dates and controlled the estimated incidence as a confounder in the regression model. Fifth, since we used Tokyo data and other geographical locations were not analysed, the results may vary under different climates. Thus, future studies in multiple locations with our proposed approach are needed. Lastly, the change of ascertainment rate and other confounders (e.g., comorbidity) were not considered in the regression model for CFR. However, we used the empirical data for severity from 25th May 2020 when the state of emergency was lifted, i.e., the end of the first wave of the epidemic in Tokyo to avoid higher ascertainment bias and overestimation of CFR.

Despite such limitations, we believe that the present study has provided comprehensive and valuable insights on the association between temperature and the characteristics of COVID-19. Higher transmissibility is likely to be seen at low temperatures, while higher severity is likely to present at high and moderately low temperatures. Our findings have important implications to public health responses, as exposure to cold and hot temperatures under a surge of COVID-19 may have the major impact on the dynamics and reduction of the burden. We also successfully provided a framework to explore the impact of meteorological factors on the transmissibility and virulence of directly transmitted diseases. Our proposed approach will be applicable for future studies on the relationships between meteorological factors and directly transmitted diseases.

Method

Epidemiological and meteorological data

The data used in the present study were from lab-confirmed, illness onset, and death cases in Tokyo. Meteorological data of temperature, relative humidity, ultraviolet radiation, and wind speed were also analysed.

We used data from 16th February 2020 to analyse the transmissibility with the regression model, as it was the earliest date of the limited publicly available dataset. To analyse severity, data in the regression model were used from 25th May 2020, when the first state of emergency was lifted in Tokyo because CFR may be underestimated given the under-ascertainment rate, and downward ascertained trend early in the epidemic. To avoid the influence of the different infectivity and severity between the previous strain and other evolved strains, e.g., variant Alpha of SARS-CoV-2 (B.1.1.7), we cut off the period in both analyses after March 202131,32.

The daily number of confirmed cases, illness onset cases, and deaths with COVID-19 in Tokyo were collected from 16th January 2020 to 19th March 2021. Confirmed data with age (decades) were also collected from 16th February 2020 to 7th April 2021. To address the measurement of overwhelmed medical situations, we obtained the daily number of cases of emergency transportation whose destination had not been determined within 20 min from the start of the Emergency Medical Services team’s request, or who had been refused by at least five medical institutions. To deal with the impact of human mobility, we resorted to Google’s COVID-19 Community Mobility Reports33, which provides three data-streams on movement in Tokyo: “residual”, “retail and recreation”, and “workplace”. All measures quantify the percentage of deviation from a baseline which indicates the median value for the day of the week during the 5 weeks from 3rd January 2020 to 6th February 2020.

Daily weather data (mean temperature (°C), relative humidity (%), solar radiation as an ultraviolet (MJ/m2), and mean wind speed (m/s)) were obtained from the Japan Meteorological Agency.

Effective reproduction number R t

The daily \(R_{t}\) estimates were derived from the daily number of confirmed cases and implemented in the “EpiNow2” package in R v4.0.2 which method accounted for the week effect and the smoothed renewal process with an appropriate Gaussian process with a squared exponential kernel34. The distribution of generation time was adopted from the earlier work35.

Nonlinear and delayed effect of temperature on R t

Non-linear and delayed effects of temperature on the transmissibility of COVID-19 were identified by using the generalized additive Gaussian model with the distributed lag non-linear model36,37.

$$log\left( {E\left( {R_{t} } \right)} \right) = \alpha + cb.temp + s\left( {time,7} \right) + intervention_{t} + mobilities_{t}$$
(1)

where \(cb.temp\) represents the nonlinear and delayed exposure-lag-response relationship between the daily \(R_{t}\) and temperature as a form of cross-basis spline function. We used a natural cubic spline with four equally spaced internal knots in the log scale in the cross-basis function38, accounting for up to 7 days of lag for the temperature to examine the lag effect from infection to secondary infection, which is referred to as generation time35. Four degrees of freedom (df) of lag were chosen by Akaike Information Criteria (AIC) to find the best-fit df for predicting missing observations, i.e., unobserved temperatures. \(s\left( . \right)\) is a natural cubic spline function. The median value of temperature for calculating relative risk (RR) was 15.3 °C. We controlled calendar dates for seasonality or long-term trend (\(time\)) as a confounder (Supplementary Fig. S5). Seven df per 380 days to \(time\) were chosen. In addition, \(R_{t}\) would be also influenced by the suppression or mitigation strategies, and other social behavioral changes due to increase in individual awareness of infection26. Therefore, we used mobility data, specifically classified into recreation, work, and residual place based on Google mobility data, assuming the three types of places as major possible sites of infection as the variables in the model involved some non-pharmaceutical interventions. To compensate above-mentioned issues other than human mobility, we reflected three categorical variables (\(intervention_{t}\)) as 0/1/2. \(intervention_{t}\) was imputed as 0 when there were interventions with low intensity on the day \(t\), 1 was denoted when the shortened business hours were requested by the Tokyo Metropolitan Government, and 2 was denoted when the state of emergency was declared. Here we did not include a variable for week effect because the framework to estimate \(R_{t}\) has implicitly accounted for the week effect34. Distributed lag non-linear model was implemented via the “dlnm” package in R v4.0.2.

Time-delay adjusted case fatality risk (CFR)

Subsequently, the association between temperature and the severity of COVID-19 was explored using CFR as a proxy of severity, and the unbiased CFR and daily CFR were estimated39. Unbiased CFR is time consistent value while daily CFR is fluctuated on every illness onset date and both accounted for the delay from illness onset to death. We assumed \(f_{s} = F_{s} - F_{s - 1}\) for \(s > 0\) where \(F_{s}\) is cumulative density function of the time-delay. The empirical time-delay distribution was fitted to lognormal, Weibull, gamma, and exponential distributions and best fit gamma distribution with mean 16.6 days and standard deviation 118.4 days by the lowest value of AIC (Supplementary Fig. S3). Here let \(\delta_{t}\), \(d_{t}\), and \(j_{t}\) be the number of illness onset dates of deaths, deceased dates of deaths, and daily new cases on day \(t\), respectively. To adjust for the time delay, we developed a framework to estimate daily CFR on an illness onset date. Then the time-delay adjusted daily CFR \(\pi_{{t_{i} }}\) on a time point \(t_{i}\) with observation \((i = 1,2, \ldots , 299)\), i.e., from 25th May 2020 to 28th February 2021, was modeled as

$$\pi_{{t_{i} }} \sim Beta\left( {shape1 = \delta_{{t_{i} }} + 1,shape2 = j_{{t_{i} }} - \delta_{{t_{i} }} + 1} \right)$$
(2)
$$d_{t} \sim Poisson\left( {d^{\prime}_{t} } \right)$$
(3)
$$d^{\prime}_{t} = \mathop \sum \limits_{s = 1}^{t - 1} \delta_{s} f_{t - s}$$
(4)

The daily CFR was modelled to be generated by beta posterior distribution (Eq. (2)). We convoluted \(f_{t}\) with \(\delta_{t}\) to obtain the expected number of illness onset dates of deceased cases \(d^{\prime}_{t}\) and \(d_{t}\) was assumed to follow a Poisson distribution (Eqs. (3) and 4). To deal with the latent variable caused by the convolution, the non-parametric back-projection based on Expectation–Maximization-Smoothing algorithm40,41 was conducted by using the “surveillance” package in R v4.0.2.

In addition, unbiased CFR was estimated as the baseline of the daily CFR estimates. \(\pi\) denoted the parameter representing the unbiased CFR on the latest day \(t\), the likelihood of the estimate \(\pi\) was given as

$${\text{L}}\left( {\pi ;j_{t} ,\theta } \right) = \mathop \prod \limits_{{t_{i} }} \left( {\begin{array}{*{20}c} {\mathop \sum \limits_{t = 1}^{{t_{i} }} j_{t} } \\ {D_{{t_{i} }} } \\ \end{array} } \right)\left( {\pi \frac{{\mathop \sum \nolimits_{t = 2}^{{t_{i} }} \mathop \sum \nolimits_{s = 1}^{t - 1} j_{t - s} f_{s} }}{{\mathop \sum \nolimits_{t = 1}^{{t_{i} }} j_{t} }}} \right)^{{D_{{t_{i} }} }} \left( {1 - \pi \frac{{\mathop \sum \nolimits_{t = 2}^{{t_{i} }} \mathop \sum \nolimits_{s = 1}^{t - 1} j_{t - s} f_{s} }}{{\mathop \sum \nolimits_{t = 1}^{{t_{i} }} j_{t} }}} \right)^{{{\mathop \sum \limits_{t = 1}{t_{i} }} j_{t} - D_{{t_{i} }} }}$$
(5)

where \(t_{i}\) and \(D_{{t_{i} }}\) represent and the cumulative number of deaths until the reported day \(t_{i}\), respectively39,42. The parameter was estimated by using Markov chain Monte Carlo (MCMC) method in a Bayesian framework with the flat prior \(\left( {Uniform\left( {0,1} \right)} \right)\). We employed Hamiltonian Monte Carlo algorithm with No-U-Turn-Sampler and obtained five chains of 600 thinned samples from 30,000 MCMC iterations where the first 1000 samples of the chains were discarded as burn-in. The MCMC simulations were performed using the “rstan” package in R v4.0.2.

Nonlinear and delayed effect of temperature on time-delay adjusted CFR

We fitted a gamma regression combined with DLNM to estimate the association between temperature and the time-delay adjusted daily CFR \(\pi_{{t_{i} }}\) with illness onset dates taking into account the delays in effect of temperature.

$$\pi_{{t_{i} }} \sim Gamma\left( {\mu_{{t_{i} }} } \right)$$
(6)
$$log\left( {\mu_{{t_{i} }} } \right) = \beta + cb.temp + hospital_{t} + age_{t} + DOW + holiday + s\left( {time,5} \right)$$
(7)

where \(cb.temp\) represents cross-basis spline function of temperature by a natural cubic spline with four equally spaced internal knots in the log scale in each cross-basis function, accounting for up to 14 days of lag to temperature to examine the period between infection to illness onset, i.e., incubation period which has previously been explored elsewhere43. We considered the 99% upper bound of the incubation period. We also adjusted for the days of the week (\(DOW\)), holidays \((holiday)\), and calendar days for seasonality and long-term trend (\(time\)). The smooth function of date (\(time\)), to allow for changes due to seasonal effects and demographic shift or other slow change not captured in the covariates, comprised a natural cubic spline of date with five degrees of freedom per 299 days (Supplementary Fig. S6). \(\beta\) is the intercept. The median value of temperature for calculating RR was 18.6 °C. Daily age distribution of infected cases with an illness onset day is also critical for CFR as a confounder, i.e., age and age-specific infection fatality risk has an exponential relationship44. Because only age distribution with reported dates was publicly available, we back-projected the illness onset date of cases who were over 70 years and in all age groups from the reported dates of cases to calculate the proportion of the daily number of cases over 70 years out of the daily number of cases in all age groups. The time delay between illness onset to reporting is fit as Weibull distribution and the parameters were adopted from the previous study41. In addition, we used the time-series data describing the pressure on medical institutions as \(hospital_{t}\) because whether the healthcare system is overloaded or not is a critical factor for CFR.

We conducted sensitivity analysis corresponding to the length of lag and possible meteorological confounders to assess the robustness of the models. As for the lag, the maximum lag day of temperature was set to 5 and 6 to examine the sensitivity of the effect in DLNM for the analysis of transmissibility. For the severity, the maximum lag day of temperature was set to 10 and 12. Regarding meteorological factors as confounders, relative humidity, wind speed, and ultraviolet were included for the analysis of transmissibility, while we considered only relative humidity for the analysis of severity.