key: cord-0823489-a19y66fe
authors: Landier, Jordi; Paireau, Juliette; Rebaudet, Stanislas; Legendre, Eva; Lehot, Laurent; Fontanet, Arnaud; Cauchemez, Simon; Gaudart, Jean
title: Cold and dry winter conditions are associated with greater SARS-CoV-2 transmission at regional level in western countries during the first epidemic wave
date: 2021-06-17
journal: Sci Rep
DOI: 10.1038/s41598-021-91798-9
sha: 72a05d5f81dfebea6580e75f0049c2af4fb9bb22
doc_id: 823489
cord_uid: a19y66fe

Higher transmissibility of SARS-CoV-2 in cold and dry weather conditions has been hypothesized since the onset of the COVID-19 pandemic but the level of epidemiological evidence remains low. During the first wave of the pandemic, Spain, Italy, France, Portugal, Canada and USA presented an early spread, a heavy COVID-19 burden, and low initial public health response until lockdowns. In a context when testing was limited, we calculated the basic reproduction number (R(0)) in 63 regions from the growth in regional death counts. After adjusting for population density, early spread of the epidemic, and age structure, temperature and humidity were negatively associated with SARS-CoV-2 transmissibility. A reduction of mean absolute humidity by 1 g/m(3) was associated with a 0.15-unit increase of R(0). Below 10 °C, a temperature reduction of 1 °C was associated with a 0.16-unit increase of R(0). Our results confirm a dependency of SARS-CoV-2 transmissibility to weather conditions in the absence of control measures during the first wave. The transition from summer to winter, corresponding to drop in temperature associated with an overall decrease in absolute humidity, likely contributed to the intensification of the second wave in north-west hemisphere countries. Non-pharmaceutical interventions must be adjusted to account for increased transmissibility in winter conditions.

Region selection. The six countries (USA, Canada, Spain, Italy, France and Portugal) included 128 regions or states. After applying exclusion criteria, the analysis was performed over 63 regions with sufficient exponential growth and complete weather data ( Fig. 1 , Figure S1 ): 24 states in the USA, two states in Canada, 11 regions in France, 12 autonomous communities in Spain, 13 regions in Italy and one region in Portugal (hereafter 'regions').

Data description. The maximal daily death count was 39 deaths/day in median (interquartile range (IQR) = 18-83, max = 779) and occurred 25 days (median; IQR = 19-43) after the start of the lockdown (Fig. 2) . Of note, larger reductions in human mobility patterns after lockdown as assessed from Google mobility led to shorter delay (Spearman correlation coefficient = 0.63, p < 0.0001, Figure S3 ). The delay between the start of lockdown and the peak in daily death count was reduced to 19 days (median; IQR = [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] for South European countries with nationwide lockdowns and strong mobility reductions (Italy, France, Spain and Portugal, > 60% in average, Figure S3 ). R 0 was estimated from death counts over an exponential growth period that had a median duration of 11 days (IQR = 9-14, range = 5-19) . The R 0 estimation period started 5 days (median across regions; IQR = 0-8.5) and ended 16 days after the date of lockdown (IQR = 11-21, range . Figure S4 presents detailed calculation periods for all regions.

The median R 0 value was 2.58 (IQR = 2.08-2.66). R 0 estimates were lowest (< 1.5) in two regions in France and one in Spain, the USA and Italy (respectively in Centre-Val de Loire, Nouvelle Aquitaine, La Rioja, Alabama and Abruzzo). R 0 was highest (> 4.0) in New York (USA), Lombardia and Piemonte (Italy), Castilla-La Mancha (Spain) and Ontario (Canada) (Fig. 3A,B) . R 0 values exhibited significant spatial autocorrelation (Moran's I = 0.20, p = 0.0087).

Demographic and climatic variables were heterogeneous between countries (Figs. 3, 4). Regional population density (excluding areas with < 5 inhabitants/km 2 ) ranged between 82 and 1056 inhabitants/km 2 (median = 271, Fig. 3C ,D) and was overall higher in Italy (Fig. 4B) . Weather covariates were introduced as mean, minimum or maximum value over the estimated transmission period. Across the 63 regions, mean AH during the transmission period was 4.98 g/m 3 www.nature.com/scientificreports/ median (IQR = 7.1-11.5, range = − 2.0- 19.9 ). Distance to the first region with 10 cumulative deaths was 406 km in median (IQR = 228-794) and much larger in the USA compared to European countries (Fig. 4F ).

Factors associated with R 0 . Each variable was included in a univariate model assuming linear (Fig. 5 , Table S2 ) or non-linear (spline smoothing, Table S3 ) relationships with R 0 . We identified significant relationships between R 0 value and temperature, AH or DP, average cumulative precipitation per day, but not with mean wind speed over the transmission period. www.nature.com/scientificreports/ www.nature.com/scientificreports/ Multivariable models were constructed according to a directed acyclic graph ( Figure S2 ). After adjustment for distance to the first affected region in the country, percentage of population over 80 and population density, there was a consistent relationship between increasing temperature, AH or DP, and decreasing R 0 . No relationship was found with average cumulative precipitation per day.

Mean temperature led to the model with the largest deviance explained compared to models including minimum or maximum temperature, or any AH or DP measurements (Table 1 ). In this model, a tenfold (+ 1 log10 unit) increase in population density was associated with a + 0.67 R 0 unit increase (95% Confidence Interval (95%CI) = 0.05-1.28). The proportion of population aged 80 years and older was not significantly associated with R 0 . The relationship between mean temperature and R 0 was not linear. A strong, nearly linear drop of approximately 1.0 R 0 unit was observed between 2.5 and 10 °C, and a plateau beyond 10 °C (Figs. 6, 7) . The residuals from this model did not exhibit significant spatial autocorrelation (Moran's I = 0.054, p = 0.215).

Mean AH and mean DP values exhibited similar profiles with increasing values associated with decreasing R 0 values (Table 1 , Fig. 6 ). It seemed that the relationship between R 0 and AH or DP did not reach a plateau. Assuming a linear relationship, a 1 g/m 3 higher AH translated in a 0.15 unit lower regional R 0 , and a 1 °C higher DP translated in a 0.08 unit lower R 0 (Table S4 ). The residuals from these model did not exhibit significant spatial autocorrelation (Moran's I = 0.096, p = 0.105 for AH, Moran's I = 0.0942, p = 0.11 for DP).

In order to evaluate the importance of the 3-week delay between meteorological variables measurement (during estimated infection period) and R 0 calculation (period when corresponding deaths occurred) we studied the relationship between R 0 and meteorological variables using different delays, from 2 weeks (lag = − 1) to 8 weeks (lag = 5), with the 3-week delay used in the main analysis corresponding to the reference value (lag = 0). Using lagged weather summary values as linear univariate predictors, our statistical model found similar relationships between R 0 and temperature variables, but 0-, 1-and 5-week lags led to the strongest effect for mean AH and mean DP ( Figure S5 ). Using lagged weather summary values as non-linear predictors in the multivariate model, the shapes of the relationships between temperature, AH and DP remained similar, with a nearly linear drop reaching a plateau at values corresponding to milder winter weather/climate ( Figure S6 ). Overall, correlations were strong between weather summary observations at the different lags (Figure S7 ). Correlations of lagged temperature and humidity observations was highest between − 1, 0, 1 and 5-week lagged observations compared to other lags.

When setting the upper limit of the R 0 calculation window to 18 days after date of lockdown, 10 additional regions were excluded because they failed to meet the 10 daily deaths limit over the narrower R 0 calculation window and the analysis was conducted on 53 regions. After adjusting for population density, population over 80 and distance to the first region affected, the non-linear relationship between temperature and R 0 remained unchanged, reaching a plateau around 10 °C (spline p value = 0.0475). The shape of the relationship between AH, respectively DP, and R 0 remained similar but did not reach statistical significance (p value = 0.19). Later AH or DP summary values, corresponding to 1 week later than the estimated transmission period (i.e. 2-week lag from www.nature.com/scientificreports/ Figure 6 . Non-linear effects for weather parameters obtained in the multivariable model for R 0 analysis (see Table 1 www.nature.com/scientificreports/ R 0 calculation period), appeared to restore the relationship with R 0 (p = 0.0661, respectively p = 0.0667), although not reaching statistical significance threshold. Fitting continent-specific splines also resulted in low statistical power: only 37 regions remained in South Europe and 26 in North America. After adjusting for population density, population over 80 and distance to the first region affected, increasing humidity variables (AH or DP) were associated with a decreasing R 0 (continentspecific spline p values were 0.0915 (North America) and 0.0988 (South Europe) for AH, respectively 0.0506 and 0.1078 for DP, not reaching statistical significance). The relationship between R 0 and mean temperature was markedly different between North America (temperature increase associated with R 0 decrease, p = 0.005) and South Europe (p = 0.17). Assuming continent-specific effects of the distance to the first region affected on R 0 indicated that distance played a significant role in South Europe, but not in North America, without modifying the effect of temperature, AH or DP on R 0 .

In this study, we analysed SARS-CoV-2 propagation parameter R 0 during the first wave of the pandemic in 63 regions of six north western hemisphere countries. We showed that R 0 values were influenced by population density (+ 0.6 for a tenfold increase in density), by proximity with the first epidemic focus of the country or coast for USA (− 0.3 for a tenfold increase in distance to the first region to record 10 COVID-19 deaths), and by weather or climate conditions. For regions with mean temperatures below 10 °C during the transmission period, a linear association was observed with R 0 values: a 1 °C increase in temperature between regions was associated with a 0.16-unit decrease in R 0 . A 1 g/m 3 increase in mean AH was associated with a 0.15-unit decrease in R 0 ( Figure 7) . Similar results were obtained with DP, with a 1 °C increase associated with a 0.08-unit decrease in R 0 (Table S4 ). After adjusting for major confounders and spatial autocorrelation, our results indicate that weather conditions brought a significant contribution to drive the magnitude of the first wave, even if it was limited by an initial heterogeneous spread of the virus, which protected regions located furthest away from the first foci.

This study relied on a regional scale analysis and accounted for different dynamics within the same country. The overall epidemic wave at country level was actually the sum of diverse dynamics, as is obvious for large countries but also true for Spain, Italy or France, where regions were heterogeneously affected, as confirmed by serological studies 3, 20 . Likewise, weather or climate heterogeneity between regions of a given country was large. By analyzing distinct spatial units with heterogeneous population density and weather, we were able to assess the effect of parameters that may be otherwise confounded by country-level parameters such as response strategy, timing of the analysis period compared to the progression of the epidemic, but also age structure of the population.

The six countries were selected due to their homogeneous location and the low efficacy of their early counter epidemic responses. We are therefore as close as possible of conditions allowing R 0 estimation. This fourparameter hierarchical generalized additive model provides an explanation of up to 42% of R 0 variability across 63 affected regions of the northern hemisphere, thus providing useful insights on the drivers of the first wave, and allowing to estimate the contribution of population density and weather conditions for the next waves, when the effect of local introduction is no longer relevant.

We analyzed R 0 , a parameter derived from the exponential growth rate, which is more robust to variations in reporting standards between countries or regions and requires less stringent assumptions on reporting standards comparability. It is also less biased compared to indicators based on cumulative counts or cumulative incidence rates. Indeed, the cumulative death (or case) count in a region depends on: the duration since the start of the epidemic in the region, the basic reproduction number (subsequently the effective reproduction number) and the number of seeding events: regions were not affected simultaneously, and their numbers of seeding events were likely different (but nearly impossible to quantify). The definition of a time period over which to sum counts to fully account for the epidemic requires many assumptions. In order to try and quantify the effect of climate using case/death counts, it would be necessary to account for factors influencing R 0 as we did here, but also for factors influencing the number of seeding events and the duration of the epidemic (such as the date of introduction, the intensity of lockdown…). Our conservative assumptions excluded regions without a clearly defined exponential growth phase from this analysis, using an arbitrary 10 deaths threshold (for starting at 10 cumulative deaths and for exponential growth phase at 10 daily deaths). Given what is known of the case-fatality ratio, this threshold corresponded to >1000 infections. This may have excluded regions with a lower growth, but identifying a relationship between R 0 and temperature, AH or DP on the regions with the highest R 0 remains relevant irrespectively of the relationship in low growth regions.

We acknowledge several limits. The first is the necessity to rely on death counts to estimate R 0 during this wave due to the lack of information on actual infections from limited testing and the lack of consistency between countries for hospitalization counts. This assumption is however similar to usual assumptions for R 0 estimation based on diagnosed cases, since true number of infections remains unknown and detection occurs at variable delays after infection.

A second limit is the use of a single summary weather observation over the assumed transmission period, which may fail to express the dynamical aspects of weather. The choice of fixed, 1-week increments to study the effect variations in the weather summary window may be too coarse to capture effects for the narrowest exponential growth periods. The sensitivity analysis showed a slight increase in effect when considering weather summary values calculated one week earlier than the estimated transmission period, but the 18-day sensitivity analysis showed a stronger negative relationship between R 0 and AH summary values occuring 1-week later than the estimated transmission period. This study could also only evaluate the effect of a limited range of weather conditions on SARS-CoV-2 transmission, since the number of observations with a mean temperature below 2.5 °C or above 15 °C was low. This prevents from analyzing the effect of higher temperatures such as during www.nature.com/scientificreports/ autumn. This restriction was however necessary to achieve a minimal homogeneity of the studied regions in terms of their exposure and response to the pandemic. Finally, this analysis includes only 63 regions and may lack statistical power. The necessity to ensure that the epidemic growth of deaths was sufficient led to the exclusion of regions that may be affected, but not enough for daily deaths counts to reach 10 deaths/day. Finally, this study showed an association between SARS-CoV-2 R 0 and temperature/AH, but due to the strong correlation between AH and temperature in the seasonal conditions analysed here, it is difficult to determine which parameter is more important and they could not be analysed in combination 21 . The continent-specific analysis suggests that the relationship between AH and R 0 was more stable than that of temperature, which was strong in the US but less so in South Europe. We could not conclude whether the relationship results from a direct role of weather on individual susceptibility to viral invasion (e.g. dry nasal mucosa from indoors heating and outdoors cold) or on viral persistence/survival, or from an indirect role of climate on human behaviors (e.g. regions with cold winter favor more indoors living conditions and lead to bigger infection opportunities). Interestingly, relative humidity and temperature levels associated with experimental decrease in virus survival corresponded to the range of AH values associated with R 0 decrease in this study ( Figure S9) 8 .

UV irradiance was identified as an important parameter to inactivate viral particles and could thus also bring an important contribution to limiting viral transmission. UV data was not included in reported data from this weather-station network. However, strong correlations exist between temperature, latitude and UV irradiance, which would be difficult to disentangle in an aggregated analysis. If warmer, more humid winter conditions directly influence transmission by limiting virus survival, it is possible that part of this effect originates from higher UV irradiance.

Our study shows an important dependency of SARS-CoV-2 transmission to weather/climate, with each 1°C reduction in mean regional temperature below 10°C leading to a 0.16-unit increase in R 0 , or each 1 g/m 3 reduction in AH leading to a 0.15-unit increase in R 0 . Northern hemisphere countries experienced a new wave of SARS-CoV-2 infections during autumnal transition from summer to winter (colder temperatures and lower absolute humidity levels), while still actively maintaining control strategies 22 . When planning to adjust the level of restrictions on social activities, public health strategies need to account for the increased transmissibility of SARS-CoV-2 when/where cold and dry winter conditions are prevalent.

Study design. In order to compare the drivers of the epidemic dynamics between regions accurately, we calculated the R 0 of the virus for each region affected by the first wave of COVID-19 epidemic in six countries, using the dynamics of the daily death counts. The study focused on the four countries hit earliest in Europe and North America (Italy, Spain, France and USA). Canada and Portugal provided geographical contiguity. These six countries shared similar location in the western part of the northern hemisphere, approximately between 25 and 50° of latitude ( Figure S1 ) defining winter conditions at the start of the first wave of the epidemic. Furthermore, they also all underwent significant SARS-CoV-2 transmission in a context of low public health response.

This analysis was conducted at the first administrative subdivision of the countries, corresponding to states in USA and Canada, autonomous communities in Spain, and regions in Italy (regioni), Portugal (região) and France (regions), all referred to as "region" in this study. www.nature.com/scientificreports/ Confirmed case counts were insufficiently reliable due to overall lack of tests and different testing strategies between countries, between regions of a given country and between periods for a given region. Hospitalization counts were not available consistently at regional level across countries. Overall, COVID-19 deaths were preferred as they were less likely to have different definitions within the same region or country and to undergo significant changes in their definition over the study period. R 0 is an indicator of the speed of progression of the outbreak, and was therefore less likely to be biased compared to indicators based on cumulative counts or cumulative incidence rates.

Estimating R 0 values based on deaths counts relies on the minimal assumption that the infection-fatality rate was constant over the study period (~ 1 month) for a given region, which defined a proportional relationship between infections and deaths. As a result, R 0 remains comparable across regions reporting death counts consistently over time, even if the percentage of deaths reported across regions is variable (for example, if only deaths occurring in the hospital are reported in one region, compared to another where all deaths are accounted for). As a result, comparability in COVID-19 related deaths reporting standards was not necessary: the only prerequisite was that death counts are reported using the same method in a given region over the period analysed.

Data. Deaths from COVID-19. Regional level data on COVID-19 deaths were retrieved from data shared by national health ministries and/or public health agencies of Spain, Italy, Portugal, France, USA, and Canada (Table S1 ). Deaths were reported as daily new death counts or as a cumulative number.

Population and geographical data. Population structure by one-or five-year age groups, sex and region was retrieved from open access data shared by national institutes or administrations responsible for national statistics or demography (Table S1 ). The percentage of the region population aged ≥ 70 or ≥ 80 years was calculated in order to adjust for differences in way of life (e.g. rural regions have older population than metropolitan areas in Europe).

Region shapefiles by country were obtained from national geographical authorities, open source datasets or as provided in the coronavirus open data packages proposed by several national health agencies (Table S1 ). The region surface was either obtained from the geographical layer (land surface area), or calculated from the polygon extent. We estimated the percentage of each region surface with ≥ 5 inhabitants per km 2 based on WorldPop 2020 raster dataset 23 . Population density was estimated as the population in the region divided by the surface after excluding areas < 5 inhabitants/km 2 in order to limit the underestimation of population density when high heterogeneity existed between urban centers and sparsely populated territories within the same region (deserts, mountains, polar regions).

Weather/climate. Weather station reports since 1 January 2020 were obtained from US National Oceanic and Atmospheric Administration (NOAA) through the R-package {worldmet}. We extracted hourly temperature, relative humidity (RH), dew point temperature (DP), precipitation and windspeed observations for each station. Hourly absolute humidity (AH, in g/m 3 ) was calculated from RH and temperature based on the Clausius-Clapeyron formula 24 . Daily minimum, maximum and mean values were calculated for temperature, AH, DP, as well as cumulative sum of precipitation and mean wind speed. Days with ≥ 5 missing hourly record (no observation of any parameter for > 20% of the day) were excluded to avoid biasing the period estimation with daily reports excluding night time temperatures.

Humidity is a key parameter for respiratory viruses, and was analysed using temperature-independent physical parameters (directly observed DP or calculated AH). Relative humidity (RH) was excluded from the analysis, since it is dependent on both temperature and AH through the Clausius Clapeyron formula, it would not characterize weather conditions specifically. Detailed properties of DP, AH and RH have been discussed extensively by Babin 21 .

Each region was attributed stations based on geographic location. All available observations in weather stations of the region contributed to the regional daily average. Weather stations can be located in mountains or sparsely populated locations where they record weather conditions that differ strongly from actual populated areas of the same region. To avoid this bias, observations were assigned population-based weights: we estimated the population located within 10 km of each weather station using WorldPop 2020 data; and for each day and each region, we calculated the total population within 10 km of any station reporting data for that day. Daily station observations were weighted in proportion of the population around each station relative to the total population. Stations located within 10-km of each other were included in a single buffer and each station was assigned equal weight within the buffer.

For each weather parameter, the regional summary value was calculated over the transmission period, during which infections were assumed to have occurred. We assumed that the transmission period had the same duration as the R 0 calculation period, and occurred 3 weeks earlier according to the delay between infection and death ( Figure S8 ).

Using this approach, the weather parameters averaged over the assumed transmission period and over each region included: minimum, mean and maximum average values for temperature, AH, DP, average cumulative precipitation per day and average wind speed.

Date of lockdown definition. Lockdown date was defined using Google mobility data as the date when a decrease > 25% in workplace localization was reported and sustained over 3 days in the region 25 . All references to a «date of lockdown» hereafter refer to this definition. This definition matched national lockdowns in European countries. This simplification was necessary due to the heterogeneity in social distancing measures taken www.nature.com/scientificreports/ at regional (state) level in the USA and Canada. The objective was to exclude periods where transmission would start to slow down due to these measures.

Distance to first region with 10 cumulative deaths. In order to account for the proximal spread of the epidemic, we considered that the first region affected in a country could represent a significant source of imported cases for other regions, and that the effect would decrease with the distance, reflecting spatial autocorrelation due to proximity in the spread of the epidemic. For each of the four European countries considered, we identified the first region with 10 cumulative deaths. In Portugal, two regions reached ≥ 10 deaths on the same day, and the region with the highest count (14 vs. 12) was selected. In the USA and Canada, the distance to the first region above 10 cumulative deaths was calculated separately for East Coast and West Coast, using the limit between Central and Mountain time zones. This was necessary due to the early start of the epidemic in the state of Washington, which reached 10 deaths by 2 March 2020, 16 days before the next state (New York on 18 March). Within each country, the euclidean distance in kilometers between the main city in each region and the main city in the first region above 10 cumulative deaths was calculated.

This distance was a coarse but necessary simplification. Including more foci or accounting for actual passenger flux (rather than Euclidean distance) would rely on major assumptions requiring integration of complex data to describe connectivity between regions and/or importation rates. In addition, 10 cumulative deaths corresponded to a large number of clinical cases ensuring that by the time a region reached this level in the country, local case circulation was likely to exceed international case importation. R 0 estimation. Daily death counts were smoothed using a 5-day moving average filter to account for irregularities in data transmission and publication, which did not present a weekend effect. A 5-day span was a compromise between smoothing accidents or irregular reporting delays and flattening the slope which might lead to underestimations.

For each region, the exponential growth period was estimated using a log(deaths) = f(time) representation and r (the exponential growth rate) was extracted as the coefficient of a Poisson regression.

The exponential growth period was assumed to start when a linear increase was identified on the log(deaths) = f(time) representation for the region. A lower limit was set at the date when 10 cumulative deaths were reached. The end of the exponential growth period was defined as the date when incidence growth slowed down (indicated by a change in the slope of the log(deaths) = f(time) representation) in order to ensure that R0 was calculated over a period not impacted by the lockdown. If no inflection was identified, the end of the exponential growth period was set at maximum 28 day after lockdown, as explained below (see also Figure S8 ).

The lower boundary of 10 cumulative deaths was set to avoid early stochasticity and limitations in available recording of the first COVID-19 deaths at regional level. The upper boundary of 28 days after lockdown was defined to avoid the influence of lockdown measures on the growth rate of death counts, since we aimed to estimate SARS-CoV-2 R 0 prior to implementation of major interventions. At individual level, the median delay between infection and death was 18 days, with a large interquartile range (IQR) of 9 to 24 days 27-33 . At regional level, we assumed that reported deaths corresponded to transmission events which had occurred in median three weeks earlier, and defined a 28-day boundary to cover the upper boundary of the IQR. R 0 was calculated for each region using the generation time method assuming a gamma distribution with parameters 7 and 5.2 for SARS-CoV-2 generation time 34, 35 . In order to improve the adjustment of the regression, the start and end dates of the calculation period were allowed to shift by up to 2 days (+/− 1 day or + 1/+ 2 days for start date if the calculation period began at the date of 10 cumulative deaths) using the built-in function sensitivity.analysis() of the package {R0} 36 .

Regions not connected to mainland country (islands, Alaska…) were not included to avoid effects related to isolation.

Regions reporting less than 10 cumulative deaths from SARS-CoV-2 before lockdown + 28 days did not have a sufficient number of deaths to exclude early stochasticity and were not included. Regions where daily death counts remained low (< 5 deaths/day) were also not included, since a low number of deaths could be linked to either the absence of an epidemic start (only sporadic cases) or an R 0 close to 1. Finally, regions with a smoothed daily death count < 10 deaths/day during the linear growth phase were excluded from the analysis. This avoided including regions where imported cases had so far failed to trigger an epidemic growth as "low risk" regions: as a result, only regions with a clearly established exponential growth phase were considered for this analysis (given current knowledge of the case-fatality ratio, 10 deaths correspond to at least 1000 cases).

One Italian region was excluded because there was no weather station data available except for 2 stations > 2000 m altitude.

Analysis of the relationship between R 0 and weather parameters. A directed acyclic graph (DAG) was constructed using Dagitty v3.0 web-based application (http:// www. dagit ty. net/ dags. html) in order to visualise the relationships between R 0 and explanatory covariates ( Figure S2 ). Dependence and independence assumptions were verified using Spearman correlation coefficient. AH or DP are also strongly correlated to temperature 21 . The different weather covariates were observed during the estimated transmission period, which corresponded to the R 0 calculation period lagged by 3 weeks to account for delay between infection and death. Weather covariates were included separately in the models. AH, DP and temperature were included as mean www.nature.com/scientificreports/ values over the period, but also minimum and maximum in order to account for the effect of extreme days. A generalized additive mixed model (gamm) regression was used to evaluate the effects of climate, population and other determinants on the value of R 0 , using a Gaussian distribution and the identity link function (package {mgcv}). A country-level random effect was included to account for within-country correlations. Canada presented with only two regions and was grouped with the USA for random effect, while the single region in Portugal was grouped with Spain. Significance threshold was set at 0.05. Models were compared using the percentage of deviance explained and Akaike's information criterion (AIC). Univariate analyses were conducted assuming linear and non-linear effects, using B-splines to model non-linear effect of covariates. Linear approximation was found acceptable for log10(population density) since it minimized AIC while causing a minor loss in percentage deviance explained over splines. Non-linear effects were retained for weather variables and distance to the first focus (minimal AIC and large increase in percentage deviance explained for splines compared to linear). We verified presence/absence of spatial autocorrelation in R 0 values and in the final model residuals by calculating Moran's I. Sensitivity analyses. First, we assessed the importance of the 3-week delay for weather variables (corresponding to the delay between transmission period and R 0 calculation period based on death count exponential growth). For this, we tested a variety of lags from − 1 to 5 weeks from that estimated transmission period (i.e. 2-8 weeks from the R 0 calculation period). The longer lags were likely irrelevant for actual transmission but aimed at identifying climate trends rather than weather influence. We included lagged weather variables as linear explanatory variables in the univariable hierarchical model or as non-linear explanatory variables in the multivariable hierarchical generalised additive model. Second, we assessed the effect of the 28-day upper boundary to the exponential growth period by using a narrower window ending 18 days after date of lockdown. In this analysis, the exponential growth period ended at the date when incidence growth slowed down (indicated by a change in the slope of the log(deaths) = f(time) representation) or 18 days after lockdown if no inflexion was identified. We recalculated R 0 for regions where the linear growth period retained for the main analysis extended beyond the 18-day limit, and followed the same plan as the main analysis.

Third, we assessed possible continent specific effects by fitting continent-specific splines for weather variables or for the distance to the first affected region in the final multivariate models.

All data used in this study is available from public data sources listed in Table S1 . The only exception corresponded to early death counts at regional level in France which were obtained directly from the French Public Health Agency but have been publicly released since. The data table of regional-level values (R 0 estimates, weather summary, population density etc.) used for the hierarchical generalized additive model analysis of the relationship between regional level R 0 and weather covariates is provided in a supplementary csv file.

No custom code was used for this study beyond usual application of standard functions of softwares or R packages. Specific code sections wil be made available upon reasonable demand. 

COVID-19 in New Zealand and the impact of the national response: A descriptive epidemiological study

Response to COVID-19 in South Korea and implications for lifting stringent interventions

Prevalence of SARS-CoV-2 in Spain (ENE-COVID): A nationwide, population-based seroepidemiological study

Population-based seroprevalence surveys of anti-SARS-CoV-2 antibody: An up-to-date review

The estimation of the basic reproduction number for infectious diseases

Seasonality of respiratory viral infections

Effect of environmental conditions on SARS-CoV-2 Stability in human nasal mucus and sputum

Mechanistic theory predicts the effects of temperature and humidity on inactivation of SARS-CoV-2 and other enveloped viruses

Airborne SARS-CoV-2 is rapidly inactivated by simulated sunlight

Transmissibility and transmission of respiratory viruses

The effect of climate on the spread of the COVID-19 pandemic: A review of findings, and statistical and modelling techniques

Comment on A. annua and A. afra infusions vs. Artesunate-amodiaquine (ASAQ) in treating Plasmodium falciparum malaria in a large scale, double blind, randomized clinical trial

Effects of temperature and humidity on the spread of COVID-19: A systematic review

Impact of climate and ambient air pollution on the epidemic growth during COVID-19 outbreak in Japan

Association of social distancing, population density, and temperature with the instantaneous reproduction number of SARS-CoV-2 in counties across the United States

Seasonality and uncertainty in global COVID-19 growth rates

Susceptible supply limits the role of climate in the early SARS-CoV-2 pandemic

Immune life history, vaccination, and the dynamics of SARS-CoV-2 over the next 5 years

A framework for research linking weather, climate and COVID-19

4.5% of population in metropolitan France developped antibodies against SARS-CoV-2: First results from the national survey EpiCov

Use of weather variables in SARS-CoV-2 transmission studies

A cross-sectional and prospective cohort study of the role of schools in the SARS-CoV-2 second wave in Italy

Impact of meteorological conditions and air pollution on COVID-19 pandemic transmission in Italy

R: A language and environment for statistical computing. R Foundation for Statistical Computing

Pathophysiology, transmission, diagnosis, and treatment of coronavirus disease 2019 (COVID-19): A Review

Epidemiology and transmission of COVID-19 in 391 cases and 1286 of their close contacts in Shenzhen, China: A retrospective cohort study

Presenting characteristics, comorbidities, and outcomes among 5700 patients hospitalized with COVID-19 in the New York City Area

Risk factors associated with mortality among patients with COVID-19 in intensive care units in Lombardy, Italy Supplemental content

es/ QueHa cemos/ Servi cios/ Vigil ancia Salud Publi caREN AVE/ Enfer medad esTra nsmis ibles/ Docum ents/ INFOR MES/ Infor mes COVID-19/Informe n o 32. Situación de COVID-19 en España a 21 de mayo de 2020

Characteristics of SARS-CoV-2 patients dying in Italy Report based on available data on December 16th

Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: A retrospective cohort study

Estimating the burden of SARS-CoV-2 in France

How generation intervals shape the relationship between growth rates and reproductive numbers

The R0 package: A toolbox to estimate reproduction numbers for epidemic outbreaks

Factors associated with the spatial heterogeneity of the first wave of COVID-19 in France: a nationwide geo-epidemiological study

The authors would like to thank Laurax Simgiane Ferbliegas from Hab' lab Marseille and L. Jae for help and support in retrieving multi-country regional data.

J.L., J.G., S.C. and A.F. designed the study. J.L. conducted the analysis with methodological contributions of J.P., E.L., and L.L. J.L. wrote the first version of the manuscript with contributions of S.R. and J.G. All authors interpreted results and contributed to the manuscript. All authors revised and approved the manuscript.

No specific funding was obtained to conduct this study.

The authors declare no competing interests.

The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s41598-021-91798-9.Correspondence and requests for materials should be addressed to J.L.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.