An Epidemiological Study to Investigate Links between Atmospheric Pollution from Farming and SARS-CoV-2 Mortality

Exposure to atmospheric particulate matter and nitrogen dioxide has been linked to SARS-CoV-2 infection and death. We hypothesized that long-term exposure to farming-related air pollutants might predispose to an increased risk of COVID-19-related death. To test this hypothesis, we performed an ecological study of five Italian Regions (Piedmont, Lombardy, Veneto, Emilia-Romagna and Sicily), linking all-cause mortality by province (administrative entities within regions) to data on atmospheric concentrations of particulate matter (PM2.5 and PM10) and ammonia (NH3), which are mainly produced by agricultural activities. The study outcome was change in all-cause mortality during March–April 2020 compared with March–April 2015–2019 (period). We estimated all-cause mortality rate ratios (MRRs) by multivariate negative binomial regression models adjusting for air temperature, humidity, international import-export, gross domestic product and population density. We documented a 6.9% excess in MRR (proxy for COVID-19 mortality) for each tonne/km2 increase in NH3 emissions, explained by the interaction of the period variable with NH3 exposure, considering all pollutants together. Despite the limitations of the ecological design of the study, following the precautionary principle, we recommend the implementation of public health measures to limit environmental NH3 exposure, particularly while the COVID-19 pandemic continues. Future studies are needed to investigate any causal link between COVID-19 and farming-related pollution.

Several studies have investigated effects of environmental and meteorological factors on the dissemination and severity of viral respiratory infections [3][4][5][6]. A 2018 paper provided evidence that air pollution from coal burning had exacerbated the "Spanish flu" pandemic of 1918 [7]. Although most COVID-19 patients develop mild or no symptoms, more severe and life-threatening symptoms such as pneumonia (often leading to acute respiratory distress syndrome), vascular inflammation and thrombosis, myocarditis, and cardiac arrhythmia, also occur, typically associated with excessive inflammation and cytokine storm [8].
Long-term and short-term exposure to atmospheric particulate matter ≤10 µm (PM 10 ) and ≤2.5 µm (PM 2.5 ) have been linked to respiratory and cardiovascular events [9]. A recent US study found that an increase of just 1 µg/m 3 in long-term exposure to PM 2.5 was associated with an 8% increase in the COVID-19 death rate [10]. Studies in Italy have identified associations between COVID-19 and COVID-19-related death and exposure to PM 10 and PM 2.5 [11][12][13][14], while tropospheric nitrogen dioxide (NO 2 ) in northern Italy was associated with levels of SARS-CoV-2 infection [15]. The COVID-19 pandemic in Italy started in the Po Valley of northern Italy-one of the most polluted areas in the world [16]where intensive livestock rearing and heavy use of fertilizers make major contributions to atmospheric pollution [17]. In December 2019 in the Po Valley Regions of Lombardy, Piedmont, Emilia-Romagna and Veneto, there were, respectively, 1,543,639; 824,801; 627,627 and 824,112 cattle, and 3,984,633; 1,121,723; 1,377,527 and 717,557 pigs. In comparison, the island Region of Sicily had 387,619 cattle and 45,152 pigs [18]. As regards nitrogen fertilizer use, the 2019 figures (tonne/km 2 /year) were 12.92 in Lombardy, 8.27 in Piedmont, 18.64 in Veneto and 19.50 in Emilia-Romagna, compared to 5.20 in Sicily [19].
Moreover, many international research groups have conducted descriptive analyses evaluating the impact of the COVID-19 pandemic on air quality during lockdown, mostly in urban areas. They reported a temporary improvement in air quality and a decrease in atmospheric emissions from the transportation sector [20,21].
Agricultural activity produces many pollutants, mainly NH 3 and PM 10 and PM 2.5 , also, whereas NH 3 can be used as a proxy for all the pollutants produced by intensive agricultural activities [22][23][24]. In Italy, agriculture is the main source of NH 3 emissions, with an estimated 362.18 kilotonnes/year, accounting for 94.3% of the total [22]. In Lombardy in 2017 (latest year for which figures are available) the proportion of NH 3 emissions from agriculture varied by province (68-99%) [25].
Analysis of all-cause mortality data indicates that COVID-19-related deaths are underestimated [2,26]. A study that analysed 22 countries reported that in Italy COVID-19-related deaths were underestimated by 30% with similar underestimates in several other countries [27]. For this reason, all-cause mortality may represent the epidemiologic indicator of choice for assessing the mortality impact of COVID-19 [2,28,29] from which excess deaths can be estimated as proxy of COVID-19 mortality [28]. Moreover, factors such as mobility, economic activities, social interactions, level of globalization of an investigated area, and total international import and export (used as synthetic parameter), have been considered important drivers of virus spread [30].
We hypothesized that long-term exposure to farming-related air pollutants might predispose to an increased risk of COVID-19-related death. The aim of this study was, 3 of 12 therefore, to investigate if exposure to farming-related atmospheric pollution may worsen the effect of SARS-CoV-2 on mortality. To this end, we performed an ecological study to find out whether the results justified additional analytic studies to identify modifiable risk factors to treat for limiting the effect of the COVID-19 on mortality. We included in our model the effect of total international import-export and of gross domestic product as a proxy of economic development by province.

Study Design
We conducted an ecological study to investigate whether interaction between SARS-CoV-2 infection and exposure to farming-related air pollutants worsened the effect of SARS-CoV-2 on mortality at the level of the provinces (administrative entities within regions) in the Italian Regions of Lombardy, Emilia-Romagna, Piedmont, Veneto and Sicily (see Figure 1). We investigated the period of the first COVID-19 wave (March-April 2020) in comparison with the same two-month period in the years 2015 to 2019. Mortality differences between these two periods were considered a proxy for mortality due to COVID-19. We hypothesized that long-term exposure to farming-related air pollutants might predispose to an increased risk of COVID-19-related death. The aim of this study was, therefore, to investigate if exposure to farming-related atmospheric pollution may worsen the effect of SARS-CoV-2 on mortality. To this end, we performed an ecological study to find out whether the results justified additional analytic studies to identify modifiable risk factors to treat for limiting the effect of the COVID-19 on mortality. We included in our model the effect of total international import-export and of gross domestic product as a proxy of economic development by province.

Study Design
We conducted an ecological study to investigate whether interaction between SARS-CoV-2 infection and exposure to farming-related air pollutants worsened the effect of SARS-CoV-2 on mortality at the level of the provinces (administrative entities within regions) in the Italian Regions of Lombardy, Emilia-Romagna, Piedmont, Veneto and Sicily (see Figure 1). We investigated the period of the first COVID-19 wave (March-April 2020) in comparison with the same two-month period in the years 2015 to 2019. Mortality differences between these two periods were considered a proxy for mortality due to COVID-19.

Pollutant Exposure
We assessed atmospheric concentrations of PM 2.5 , PM 10 , NH 3 and NO 2 . The first three pollutants were considered because linked to farming, while NO 2 was included because it could be a confounder or effect modifier. Data on atmospheric PM 2.5 , PM 10 and NO 2 levels were accessed from the national monitoring system coordinated by the Italian National Institute for Environmental Protection and Research (ISPRA) [31]. We used data from the monitoring stations managed by Regional Environmental Protection Agencies (ARPAs) provided at a province level. To estimate long-term exposure, we averaged daily concentrations from 2016 to 2019. We excluded from the analysis two provinces for which PM 10 and NO 2 annual concentrations were not available for the period under study. We did not use 2020 data, as they could have biased estimates because of changes in pollutant levels caused by lockdown [20,21,32,33]. As NH 3 concentrations provided by from ARPAs are monitored by too few stations to be representative of the areas being studied, we estimated concentrations of atmospheric NH 3 at the province level as well [25,[34][35][36][37][38]. More in depth, ARPAs estimate NH 3 emissions for each province and for each known source (e.g., pigs) by multiplying the average number pigs per year present in the province by an estimate of the average production by each animal (as kg NH 3 /animal/year). A similar method is applied to each type of livestock present and also other NH 3 sources (e.g., fertilizer use). Total NH 3 emissions were obtained by summing estimated emissions from all known sources. We next derived provincial exposure indexes from each ARPA estimate by dividing the total estimated quantities of emissions (tonne/year) by the area (km 2 ) of the province, thereby taking into account variations in the province area. We used data for the latest available year as proxy for long-term exposure: 2017 for Lombardy and Emilia-Romagna Backes et al. [39] showed that the atmospheric concentration patterns of NH 3 are in line with NH 3 emission patterns, indicating the provincial exposure indexes are acceptable indicators of atmospheric NH 3 . To further investigate exposure indexes as indicators of atmospheric pollution levels for NH 3 , we downloaded available data from the Italian National Statistics Institute (ISTAT) [18,19] on numbers of cattle and pigs and fertilizer used per km 2 by region. For each one of the twenty Italian regions, we downloaded data about NH 3 emissions from ISPRA [22], and we built regional exposure indexes in the same way we did for provinces. We then assessed correlations (Pearson's r) between NH 3 emissions and cattle numbers, pig numbers and fertilizer use per km 2 by region. We could only compare data at the regional level because of the unavailability of data on numbers of pigs and cattle and nitrogen fertilizer use at the provincial level.

Meteorological Variables
We hypothesized that temperature influenced COVID-19 rates over the short-term and so might have an effect on mortality. We therefore used temperature and humidity, as measured by ARPA monitoring stations in the main town of each province, as proxies for the average temperature and average relative humidity in each province, in March and April in 2020, and in March-April 2019 also.

Additional Covariates
Social deprivation is also linked to SARS-CoV-2 infection and death [40]. However, social deprivation indexes at the provincial level were not available, so we assumed that the gross domestic product per capita was a rough proxy for social deprivation. International trade has been hypothesized to be linked to COVID-19 spread [30], and therefore, we used these variables provided by ISTAT [18,19] to create an indicator as the ratio between the sum of the international import and export activities in Euros by provinces, divided by the size of the provincial population. We also introduced in the model the gross domestic product per capita (GDPc) by province (ISTAT), as a proxy indicator of population wealth and of socio-economic status.
High body mass index (BMI) and diabetes are recognized as major risk factors for SARS-CoV2 death [41,42] and could bias estimates of links between farming-related pollutants and SARS-CoV-2 mortality. Provincial level data were unavailable, so we accessed regional level data from ISTAT [43] and assessed correlations (Pearson's r) with NH 3 to provide indications as to the possible effects of BMI and diabetes on our findings. Variation in lockdown periods across provinces could also bias the results. However, in Italy during the first COVID-19 wave, days of lockdown were closely similar throughout the country. Similarly, regulations to limit contagion (mask-wearing, social distancing) applied to the whole country during the two-month study period.

Outcome
All-cause daily mortality data provided by ISTAT [44] were the primary study outcome. As mortality data are available by municipality (administrative entities within provinces), we calculated total mortality by province, for March-April 2020, as the sum of the daily death counts in all the municipalities in each province. We also estimated average mortality in March-April of 2015-2019, as the average over the five years of the sum of daily death counts in all the municipalities of each province.

Statistical Methods
We first sought correlations (Pearson's r) between levels of the atmospheric pollutants being studied (PM 2.5 , PM 10 , NH 3 and NO 2 ). We next analyzed associations between changes in total mortality from March-April 2015-2019 to 2020 and atmospheric pollutant levels, humidity, temperature, and population density (study covariates). To this end, we used negative binomial regression models, that included a population size offset, and which estimated mortality rate ratios (MRRs) in relation to study covariates. MRR is the ratio of the mortality rate for a specific value of a covariate relative to reference. Thus, for the covariate "period" MRR was the mortality in March-April 2020 relative to mortality in March-April in 2015-2019 and is a good proxy for the increase in mortality due to COVID-19. We were particularly interested in MRRs for the interaction between period and the other covariates (mainly pollutant levels) as these are an estimate of the mortality associated with the pollutant over and above that due to COVID-19 alone.
We first ran basic models that included period, one covariate and an interaction term between period and covariate. We next ran complete models that included all covariates and all interaction terms between period and covariates.
Information on PM 10 levels was available for all provinces, while PM 2.5 data were unavailable for four provinces. We therefore ran two sets of models, one that included all 43 provinces but excluded PM 2.5 and another that included PM 2.5 and the 39 provinces for which PM 2.5 was available. The R statistical package [45], version 4.0.2, was used to perform the analyses. As regards emissions, mean across-province exposure index estimates for NH 3 was 2.41 tonne/km 2 /year. The NH 3 exposure index ranged from 0.12 (Verbano-Cusio-Ossola, Piedmont) to 10.3 tonne/km 2 /year (Cremona, Lombardy) (data not shown). Table 1 reports Pearson correlation coefficients for PM 10 , PM 2.5 , NO 2 , NH 3, importexport and gross domestic product. Table 2 presents modelling results for all 43 provinces included in the study: the MRRs represent percentage increases in overall death rate in relation to increments of 1 µg/m 3 in PM 10 , and NO 2 , 1 tonne/km 2 increments in NH 3 , a 1 • C increment in atmospheric temperature, a 1% increment in relative humidity, a one unit increase in population density, a 1000 Euros increment in international import-export per capita and a 1000 Euros increment in GDPc (Gross Domestic Product per capita).   * MMRs are percentage increases in overall death rate associated with 1 µg/m 3 increase in atmospheric PM10, 1 µg/m 3 increase in atmospheric NO 2 , 1 tonne/year/km 2 increases in emissions of NH 3 , 1 • C increase in atmospheric temperature, 1% increase in relative humidity, one unit increase in population density, a 1000 Euros increment in international import-export per capita and a 1000 Euros increment in Gross Domestic Product per capita. § Basic models include only period, a single covariate and its interaction with period. ‡ Complete models include all covariates together with their interactions with period.

As compared to mortality documented in
The basic regression models identified a significantly increased in MMRs for interactions between period and PM 10 (MRR: 1.030; 95%CI: 1.009-1.052), NO 2 (MRR: 1.017; 95%CI: 1.001-1.032) and NH 3 (MRR: 1.093; 95%CI: 1.039-1.149), indicating that high levels of these pollutants were significantly linked to an increase in total mortality over and above that due to period (proxy for presence of . None of the individual covariates was significantly associated with MRR. In the complete model, MRRs for the interaction between NH 3 and period (MRR: 1.069; 95%CI: 1.006-1.136) remained significant, indicating that an increase of 1 tonne/km 2 /year in NH 3 emissions was significantly associated with a 6.9% increase in all-cause death over and above that associated with period (proxy for COVID -19). MRRs for interaction of period with PM 10 and NO 2 were reduced, compared to the basic model, and became non-significant. Table 3 shows results for the 39 provinces for which PM 2.5 data were available. The MRRs were similar to those obtained for all 43 provinces. As regards the complete model, only the MRR for interaction between period and NH 3 was significant (1.072; 95%CI: 1.001-1.151) corresponding to a 7.2% increase (over and above that for period alone) in mortality for each tonne/km 2 /year increase in NH 3 emissions. * MMRs: percentage increases in overall death rate associated with 1 µg/m 3 increase in atmospheric particulate matter (PM 2.5 or PM 10 ), 1 µg/m 3 increase in atmospheric NO 2 , 1 tonne/km 2 /year increases NH 3 , 1 • C increase in atmospheric temperature, 1% increase in relative humidity, one unit increase in population density, a 1000 Euros increment in international import-export per capita and a 1000 Euros increment in Gross Domestic Product per capita. § Basic models include only period, a single covariate and its interaction with period.ˆComplete models include all covariates and their interactions with period.
Correlation analysis found the following correlations: 0.92 between provincial exposure index for NH 3 and number of pigs per km 2 ; 0.98 between provincial exposure index for NH 3 and number of cattle per km 2 ; 0.72 between the provincial exposure index for NH 3 and fertilizer use per km 2 . These findings suggest that emission levels are a reasonable surrogate for atmospheric levels.
In addition, correlations were: −0.87 between BMI and the regional exposure index for NH 3 and −0.77 between diabetes and the regional exposure index for NH 3 . Given the strength of the association of COVID-related mortality increase (6.9%) due to NH 3, it is unlikely that BMI and diabetes would have had biased this association.

Discussion
Using basic models, which excluded most covariates, we found that higher NH 3 emissions-whose major source is agriculture-were associated with all-cause mortality increases in 2020 compared to previous years, over and above those due to period (proxy for COVID-19-related mortality). Higher PM 10 , PM 2.5 and NO 2 levels were also associated with all-cause mortality increases over and above those due to period.
Using complete models, which included all other studied pollutants as covariates, higher NH 3 levels remained associated with a significant all-cause mortality increase: a 1-tonne/km 2 increase in NH 3 was associated with a 6.9% increase in mortality over and above that due to period. The fact that interactions of period with NO 2 , PM 10 and PM 2.5 were not associated with excess mortality in the complete models is likely due to the fact that the levels of these pollutants correlated highly with each other.
A major study assumption was that the provincial exposure index for NH 3 based on pollutant emissions is a good proxy for the levels of that pollutant in the atmosphere. This assumption is supported by the study of Backes et al. [39], which showed that atmospheric concentration patterns of NH 3 are in line with NH 3 emission patterns. To further investigate provincial exposure indexes as indicators of atmospheric pollution levels, we assessed correlations (Pearson's r) between NH 3 emissions and cattle numbers, pig numbers and nitrogen fertilizer use per km 2 by region (provincial data unavailable): we found high or good correlations (0.98, 0.92 and 0.72, respectively), indicating that estimated emissions are consistent with the emission sources present and suggesting that these estimated emissions (standardized by area size) are a reasonable proxy for atmospheric NH 3 concentrations.
It is known that BMI and diabetes are major risk factors for COVID-19 infection and death. However, provincial data for BMI and diabetes were not available and could not be included in our models as covariates. We indirectly assessed whether the association between farming-related pollutants and SARS-CoV-2 mortality were biased by BMI and diabetes by estimating correlations between BMI and NH 3 and between diabetes and NH 3 at the regional level. As noted, given the strength of the association of the COVID-related mortality increase (6.9%) due to NH 3 , it is unlikely that BMI, diabetes or income would have exerted a major bias on this association. Any influence of population density and educational level on associations was taken into account by entering these as covariates in the models.
To our knowledge, this is the first study to investigate relationships between farmingrelated atmospheric pollution and COVID-19-related mortality. However, a 2020 study conducted in Italy found that the number of SARS-CoV-2 infections per unit of population depended on the type of rural landscape [46]. The study classified rural landscapes into (A) urban and periurban with high-intensity agriculture; (B) high-intensity agriculture; (C) medium-intensity agriculture (hilly areas); (D) low-intensity agriculture (high hills and mountains). They found that areas of less intensive agriculture (C and D) had fewer COVID-19 cases (per 10,000 inhabitants) than areas of intensive agriculture (A and B). Notwithstanding the markedly different methodologies, both our study and the mentioned one [46] identified intensive farming as a risk factor for COVID-19. Another study conducted in Italy [47] linked the short-term exposure to some atmospheric pollutants and hypothesized also a possible role for NH 3 associated with COVID-19 spread. Despite the differences in outcome-spread versus mortality and short-term versus long-term one-our study and the latter one seem to point in the same direction.
A study conducted in Italy by Bontempi et al. [30] identified international importexport as the major driver for COVID-19 spread. In our study, international import-export in the complete analysis showed a possible (although not statistically significant) association, but less than NH 3 and comparable to NO 2 . One explanation could be that international import-export may influence the place where SARS-CoV-2 starts in a country (in Italy, for example, the Codogno municipality), but once SARS-CoV-2 has started, then it rapidly spreads everywhere, independent of international import-export. A difference between our study and the one by Bontempi et al. was the different metrics used to take import-export into account. When developing our study, we used a variable created by dividing the total import-export by provincial population size instead of the rough total amount of import-export.
Similarly, the 2020 Lancet Countdown report highlighted a link between SARS-CoV-2 infection, environmental degradation and climate change and identified agriculture as one of the main sources of pollutant emissions [48].
Other studies have reported associations between exposure to livestock farming and respiratory diseases. A study in the Netherlands found that people living in high-livestockdensity areas had a significantly higher prevalence of pneumonia than those living in low-density areas [49]. A German study [50] found that people exposed to higher levels of NH 3 had poorer pulmonary health and were more likely to be sensitized to ubiquitous allergens than persons with lower NH 3 exposure. A study performed in an agricultural region of Washington State (US) found that industrial-scale animal feeding operations were associated with higher daily outdoor NH 3 levels and that forced expiratory volume in one second was lower with each interquartile increase in the previous day's NH 3 exposure; however, no association with asthma symptoms was observed [51]. The European Academy of Allergy and Clinical Immunology position paper [52] highlighted that increased risks of developing respiratory diseases and allergies were associated with animal farming and emphasized the need to protect workers from such risks.
The data provided by the above-cited studies [49][50][51][52] suggest that a link between excess COVID-19-related mortality and farming-related atmospheric pollutants is biologically plausible. The specific association we found between COVID-19-related mortality and NH 3 might be explained by the hypothesis that atmospheric NH 3 leads to the formation of an alkaline aerosol which triggers a conformation change in the SARS-CoV-2 spike that facilitates fusion of the viral envelope with the plasma membrane of target cells [53]. It is noteworthy that clusters of SARS-CoV-2 infection have been reported in association with slaughterhouses [54], which are known to have elevated levels of NH 3 [55].
We used change in all-cause total mortality in 2020 as a proxy for SARS-CoV-2 mortality, and to our knowledge, our study is the first to do this. A strength of this approach is that it avoids loss of deaths due to underreporting of the SARS-CoV-2 mortality [26]. However, the pandemic may have stretched the health resources, so that more people than usual may have died from non-COVID-19 causes, somewhat inflating the all-cause mortality.
The main limitations of the study are that, for the atmospheric pollutants that were measured, we used levels measured at the monitoring station in the chief town of the province as a proxy for exposure of the entire provincial population. In addition, since no measured levels were available for NH 3 , we used estimated emissions at the provincial level obtained from the regional environmental protection agencies as proxies for exposure of the provincial populations.

Conclusions
Although this study was ecological in design and could not identify causal links between farming-related atmospheric pollutants and COVID-19 mortality, if we follow the precautionary principle, our findings support the implementation of public health measures to limit environmental NH 3 exposure from the agricultural sector, particularly while the COVID-19 pandemic continues. Of interest, covering tanks containing slurry can speedily reduce NH 3 emissions from agriculture [56,57]. However, medium-term structural changes are required to reduce farming-related pollution. Consumption of meat, particularly red meat, is recognized as a threat to human health and well-being in part because of the effect of livestock rearing on the environment [48]. As the desire to consume meat encourages the use of intensive farming, resulting in increased emission of farming-related pollutants, we need to move towards a more sustainable agriculture with reduced livestock raising [46,48].