Next Article in Journal / Special Issue
Climate Change and West Nile Virus in a Highly Endemic Region of North America
Previous Article in Journal
Educational Differences in Smoking among Adolescents in Germany: What is the Role of Parental and Adolescent Education Levels and Intergenerational Educational Mobility?
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Modeling Monthly Variation of Culex tarsalis (Diptera: Culicidae) Abundance and West Nile Virus Infection Rate in the Canadian Prairies

1
Large Animal Clinical Sciences, Western College of Veterinary Medicine, University of Saskatchewan, 52 Campus Drive, Saskatoon, SK S7N 5B4, Canada
2
Department of Veterinary Microbiology, Western College of Veterinary Medicine, University of Saskatchewan, 52 Campus Drive, Saskatoon, SK S7N 5B4, Canada
3
Saskatchewan Ministry of Health, 3475 Albert Street, Regina, SK S4S 6X6, Canada
4
Environment Canada, Science & Technology Branch, 115 Perimeter Road, Saskatoon, SK S7N 0X4, Canada
*
Author to whom correspondence should be addressed.
Int. J. Environ. Res. Public Health 2013, 10(7), 3033-3051; https://doi.org/10.3390/ijerph10073033
Submission received: 7 June 2013 / Revised: 15 July 2013 / Accepted: 16 July 2013 / Published: 22 July 2013
(This article belongs to the Special Issue Epidemiology of West Nile Virus)

Abstract

:
The Canadian prairie provinces of Alberta, Saskatchewan, and Manitoba have generally reported the highest human incidence of West Nile virus (WNV) in Canada. In this study, environmental and biotic factors were used to predict numbers of Culex tarsalis Coquillett, which is the primary mosquito vector of WNV in this region, and prevalence of WNV infection in Cx. tarsalis in the Canadian prairies. The results showed that higher mean temperature and elevated time lagged mean temperature were associated with increased numbers of Cx. tarsalis and higher WNV infection rates. However, increasing precipitation was associated with higher abundance of Cx. tarsalis and lower WNV infection rate. In addition, this study found that increased temperature fluctuation and wetland land cover were associated with decreased infection rate in the Cx. tarsalis population. The resulting monthly models can be used to inform public health interventions by improving the predictions of population abundance of Cx. tarsalis and the transmission intensity of WNV in the Canadian prairies. Furthermore, these models can also be used to examine the potential effects of climate change on the vector population abundance and the distribution of WNV.

1. Introduction

Since the introduction of West Nile virus (WNV) into eastern North America in 1999 [1], WNV has become an endemic disease in the most of southern Canada, especially in the Canadian prairie provinces of Alberta, Saskatchewan, and Manitoba, which have had the highest human infection rate in Canada. Of the 2,315 cases in Canada in 2007, more than 98% occurred in these three prairie provinces [2]. Saskatchewan alone accounted for over half the human cases reported in Canada. Research focusing on the drivers of WNV occurrence in the prairie ecosystem is needed to inform WNV control and public health interventions.
West Nile Virus is primarily transmitted and amplified among local avian fauna and ornithophilic mosquito vectors, with occasional spillover into mammalian populations through mosquito blood feeding [3,4,5,6]. Therefore, the WNV infection rate in mosquito vectors is commonly utilized as an indicator of pathogen transmission intensity [7,8] and has been demonstrated to be a better indicator for WNV activity than surveillance of dead or infected birds [9]. Surveillance of mosquito infection rate also provides sufficient lead time for intervention and management of mosquito borne diseases [7,9].
Although WNV has been isolated from at least 59 mosquito species in North America, only a small portion including mosquitoes belonging to the genera Culex (Diptera: Culicidae) have been shown to be competent vectors [5]. Culex tarsalis Coquillett is considered to be the main vector of WNV in the Canadian prairies [3,10]. It is one of the most efficient WNV vectors evaluated in the laboratory studies [4] and the predominant species in the Canadian prairies during the summer WNV season [3]. Several biological features of Cx. tarsalis also facilitate the transmission of WNV in the enzootic cycles. Culex tarsalis can vertically transmit WNV to its offspring [11], it takes several blood meals, and it produces multiple generations per season [3]. Furthermore, it is known to feed on both avian and mammalian hosts and plays the role of the “bridge vector” which transmits WNV out of its enzootic cycle to humans and other mammalian species [12].
Environmental variables such as temperature, precipitation and habitat type influence both Cx. tarsalis population abundance and its WNV infection rate [13,14,15,16]. In addition, environmental factors may influence seasonal and spatial overlap among key hosts in the sylvatic amplification cycles. For instance, drought-induced congregation of mosquitoes and birds on shrinking wetland habitats may enhance the transmission of WNV [17]. Therefore, understanding how environmental and biotic factors influence abundance of mosquito population and WNV infection rate in mosquitoes is important in predicting the risks of WNV.
The Canadian prairies are at the northern limit of WNV distribution in the Western hemisphere. Climate and habitat suitability for Cx. tarsalis determine the distribution of this mosquito and WNV in the prairie provinces [18]. Future alterations in climate or habitat might expand the current spatial and temporal distribution of Cx. tarsalis and WNV [19,20]. A predictive model using environmental factors as the primary explanatory variables could be applied to evaluate the potential effects of climate change.
Many studies have been conducted to clarify the effects of environmental and biotic factors on the risk of WNV and predict the distribution of WNV in North America since its incursion. However these models usually cannot be applied to predict risks in other regions, particularly when those regions have different ecological dynamics or primary vector species [21]. Previous work had evaluated the effects of climate factors on the distribution of WNV by constructing a weekly model to predict the weekly variation of WNV infection rate in the Canadian prairies [22]. The objectives of this present study were to clarify how environmental and biotic factors affect the abundance and WNV infection rate of Cx. tarsalis on a monthly scale and compare this to the weekly model. In addition, due to the limitation of time scale of the climate change dataset, it was necessary to construct a monthly model fitted with the monthly dataset, the finest time scale for evaluating the effects of climate change [23].

2. Materials and Methods

2.1. Mosquito Data

WNV infection in pooled female Cx. tarsalis was determined using reverse transcription polymerase chain reaction (RT-PCR) [24,25]. Data on mosquito trap locations, abundance as the number of Cx. tarsalis per trap night, and WNV infection from across the prairie provinces were obtained from May to September for 2005 to 2008 from the Public Health Agency of Canada (PHAC). The original data were supplied to PHAC from Alberta Environment, Manitoba Public Health and Healthy Living, and Saskatchewan Ministry of Health.
Mosquito sampling sites were distributed within the different health regions across the southern half of the prairie provinces (Figure 1). The new standard miniature light traps with photocell controlled CO2 release (Model 1012-CO2; John W. Hock Company, Gainesville, FL, USA) were used for mosquito sampling. The mosquito collection period generally began in late May and lasted until the end of August (in Manitoba, the collection period ended the first week of September). During each week, the trap was operated for one night at each of the mosquito collection sites in Alberta and Manitoba, but for one to four nights per week at Saskatchewan collection sites. Monthly mean of Cx. tarsalis abundance, individuals per sampling night, was estimated for each collection site.
WNV infection rate (per 1,000 Cx. tarsalis) was computed using PooledInfRate (version 3.0), A Microsoft® Excel plug-in [26] by Maximum Likelihood (ML-IR) and minimum infection rate (MIR) methods [26,27]. WNV infection rates in the mosquito population were usually low, especially in the early transmission season. In this period, estimations of arbovirus infection rate in mosquitoes with values of zero were commonly recorded. Therefore, to achieve reasonable detection probability, large mosquito sample size for screening was required [28,29,30]. Gu and Novak [30] suggested that for a medium detection probability of 0.5, 693 mosquitoes are required.
Figure 1. Distribution of mosquito sampling sites in the Canadian prairies provinces of Alberta, Saskatchewan, and Manitoba for the period from 2005 to 2008. (a) Location of prairie provinces (grey spot) in Canada. (b) The distribution of sampling sites across the Canadian prairies ecozone.
Figure 1. Distribution of mosquito sampling sites in the Canadian prairies provinces of Alberta, Saskatchewan, and Manitoba for the period from 2005 to 2008. (a) Location of prairie provinces (grey spot) in Canada. (b) The distribution of sampling sites across the Canadian prairies ecozone.
Ijerph 10 03033 g001
In addition, we also found extremely high infection rates in some records where a positive result was obtained from pooled test which was comprised of only a few mosquitoes. Excluding observations from our dataset with low Cx. tarsalis sample sizes in the late season of WNV (usually in the early September) might also remove the records with high infection rate; however, female Cx. tarsalis are usually inseminated and preparing for hibernation in this period in which they do not take a blood meal [3] and the risk of WNV transmission is considered to be low. Based on these findings, we excluded observations with samples less than 100 female Cx. tarsalis per site per week from the analysis to prevent potential outliers and incorrect estimation of WNV infection rate resulting from small sample size [30]. The monthly mean WNV infection rate in Cx. tarsalis was calculated for each collection site.

2.2. Land Cover

The land cover dataset was derived from the Advanced Very High Resolution Radiometer (AVHRR) sensor operating on board the United States National Oceanic and Atmospheric Administration satellites. AVHRR Land Cover Digital Data was downloaded from Natural Resources Canada. Although the satellite image was taken in 1995, the land use in this area was considered to have remained relatively stable until the start of the study period [31].
Data for the prairie provinces were extracted and converted to a single 1 km2 GIS raster layer. Eleven different types of land cover categories were included in the original dataset. According to habitat utilization by Cx. tarsalis [18], land cover was simplified into forest including deciduous, transitional coniferous and mixed forests, water, barren land, agricultural land including cropland and rangeland, and urbanized area. To analyze the association between land cover, Cx. tarsalis population, and infection rate, we used 20 km radius buffer zones around collection sites, defined according to the flight distance of Cx. tarsalis [10,32] to determine the percentage of primary composition of land cover type within the buffer zones.

2.3. Weather Data

Daily mean temperature, daily maximum temperature, daily minimum temperature, and daily total precipitation were downloaded from the National Climate Data and Information Archive, Environment Canada. The daily weather datasets were used to create various predicting variables (Table 1). Daily maximum and minimum temperature were used to calculate the accumulative degree days [33] by the single sine method [34] for current month, and two and three months accumulative degree days (Table 1). We also created the monthly mean degree days predictors by using monthly mean maximum and minimum temperature [35]. For estimating the monthly mean degree days, the low temperature threshold for WNV amplification in Cx. tarsalis was set as 14.3 °C [15]. Finally, to determinate possible nonlinear effects of weather variables on Cx. tarsalis abundance and WNV infection rate, this study also created second order polynomial variables using the centered value of temperature and precipitation. Weather predictors of each climate station were interpolated by the inverse distance weighted method to create prairie-wide climate raster layers in ArcGIS. There was a total of 473 climate stations located in the Canadian prairies which were evenly distributed across the study area.
Table 1. Descriptions of variables used in both Cx. tarsalis abundance and infection rate models and relationships between explanatory variables based on the Pearson correlation and principal component analysis. Variables with the same arabic number indicated that the Pearson correlations are larger than 0.8 or have factor loading larger than 0.6 in each component of principal component analysis.
Table 1. Descriptions of variables used in both Cx. tarsalis abundance and infection rate models and relationships between explanatory variables based on the Pearson correlation and principal component analysis. Variables with the same arabic number indicated that the Pearson correlations are larger than 0.8 or have factor loading larger than 0.6 in each component of principal component analysis.
VariablesCorrelation
(>0.8)
PCA component
(Factor loading >0.6)
Variables description
LMMGLMMLMMGLMM
Monthly mean temperature1111Monthly mean temperature of the month of mosquito data collection
1 month lagged mean temperature22221 month lagged mean monthly temperature
2 month lagged mean temperature22222 months lagged mean monthly temperature
3 months mean temperature2222Including mosquito collection month, and previous 1 and 2 months
Winter mean temperature 433From December to February
Monthly mean degree days3333Monthly mean degree day of the mosquito data collection month
2 months total of monthly mean degree days3333Created by summing the monthly mean degree days of the month of mosquito data collection and previous month
3 months total of monthly mean degree days-3-3Created by summing the monthly mean degree days of the month of mosquito data collection, previous one and two months. Not applied in the LMM
Mean temperature fluctuations33, 433Monthly mean maximum temperature minus monthly mean minimum temperature
1 months accumulative degree days1111The accumulative degree days of data collection month
2 months accumulative degree days2222The accumulative degree days of data collection month and previous months
3 months accumulative degree days-2-2The accumulative degree days of data collection month and previous one and two months. Not used in the LMM
1 month lagged mean precipitation-4-41 month lagged monthly mean daily total precipitation.
Monthly total precipitation Monthly total precipitation
1 month lagged total precipitation 4441 month lagged monthly total precipitation
2 month lagged total precipitation 222 month lagged monthly total precipitation
Total precipitation of previous year 44Annual total precipitation of previous year
3 months total precipitation 44The total precipitation of mosquito collection month, and previous one and two months
LMM: Linear mixed model for predicting Cx. tarsalis abundance; GLMM: Generalized linear mixed model for predicting WNV infection rate in Cx. tarsalis; “-”: variable is not used in the model construction.

2.4. Data Analysis

Counts of Cx. tarsalis per trap site per night were transformed by ln(y + 1) to normalize the data distribution prior to analysis. Models to predict Cx. tarsalis abundance using environmental factors were constructed using linear mixed model (PROC MIXED, SAS ver. 9.2, SAS Institute, Cary, NC, USA).
A generalized linear mixed model with a log link function (PROC GLIMMIX, SAS ver. 9.2, SAS Institute, Cary, NC, USA) was then used to develop a model for Cx. tarsalis infection rate prediction. A negative binomial distribution was chosen in the Cx. tarsalis infection rate model according to a preliminary analysis where the formula “Pearson Chi-Square divided by degrees of freedom (DF)” value was close to one. This formula was used to evaluate overdispersion.
The Pearson correlation test was used to test for multi-collinearity between explanatory variables. If the correlation between any pair of variables was larger than 0.8, the more significant variable (lower −2 log likelihood) was chosen for further model construction [36]. We also conducted principal component analysis for determining the relationship between explanatory variables and compared the results with those of the Pearson correlation. Components selected were based on eigenvalues (>than 1) and the cumulative percentage of variance accounted for by the components (>than 80%). Factor loading larger than 0.6 of each variables was considered to be highly correlated with the component [37]. Health regions where mosquito collection traps were located were used as a random effect for both Cx. tarsalis abundance and infection rate models. The variable parameters were estimated by the restricted maximum likelihood and Gauss-Hermite quadrature method for Cx. tarsalis abundance and infection rate models, respectively [38].

2.5. Model Selection and Validation

Explanatory variables selected in the models of Cx. tarsalis abundance and WNV infection rate were based on the Wald test with p threshold value set at 0.05 [38,39] The backward stepwise method was adopted for variable selection. The AICc (AIC value corrected for finite sample sizes) values and AICc weights were used to assess the models fit and select the best fitted model [38,39,40]. Four-fifths of the records were randomly selected from the dataset for model construction and the remaining one-fifth of the records were used for model validation. Standardized residuals and standardized Pearson residuals were estimated for the validation dataset of Cx. tarsalis abundance and infection rate models, respectively. Values of Root Mean Square Error (RMSE) were calculated for training and validation datasets to validate and compare the predictability of the model on each dataset. Moran’s I test was used to test the spatial autocorrelation of residuals for the final models.
The final model of WNV infection rate in Cx. tarsalis was applied to create maps of the predicted WNV infection rate in the Canadian prairies using ArcGIS 10 (Environmental System Research Institute, Redlands, CA, USA). To compare the predicted Cx. tarsalis infection rate and actual human WNV incidence in the prairie provinces, maps of human WNV incidence (cases per 100,000 individuals) were created for the entire transmission season (May-September) for each health region for 2005 to 2008. Human WNV case data was retrieved from prairie provincial governmental websites [41,42,43].

3. Results

3.1. Descriptive Statistics

Out of 309 mosquito collection sites, 96 were in Alberta, 38 in Saskatchewan and 175 in Manitoba (Figure 1). Mosquito collections were conducted from late May to early September and the highest observed mean abundance of Cx. tarsalis was in July (for 2005 to 2007) or August (for 2008) (Figure 2). The highest observed mean Cx. tarsalis infection rates were in August for 2005 to 2007; there was no obvious trend observed in 2008 (Figure 2).
Figure 2. Temporal trends of (a) Monthly mean temperature (unit 1°C). (b) Monthly total precipitation (unit 1 mm). (c) Monthly mean abundance of Cx. tarsalis, ln(y+1) transformed, compared to monthly mean WNV infection rate in female Cx. tarsalis in the Canadian prairies. Error bar indicates the standard deviation of mean Cx. tarsalis infection rate.
Figure 2. Temporal trends of (a) Monthly mean temperature (unit 1°C). (b) Monthly total precipitation (unit 1 mm). (c) Monthly mean abundance of Cx. tarsalis, ln(y+1) transformed, compared to monthly mean WNV infection rate in female Cx. tarsalis in the Canadian prairies. Error bar indicates the standard deviation of mean Cx. tarsalis infection rate.
Ijerph 10 03033 g002
July was the warmest month with temperature around 20 °C. The mean temperature in July 2007 was higher than other years, ranging from one to three Celsius degree, and 8.6% higher than mean temperature in July of study period. The range of monthly temperature fluctuation was from 6.3 to 20 °C in the Canadian prairie ecozone. In 2007, total precipitation in May was highest compared to other years with a mean temperature closer to 2006 but higher than 2005 and 2008 in the Canadian prairies (Figure 2). Agriculture land cover predominated the buffered zones of mosquito collection sites; 94.3% of sampling sites had agriculture land as the primary land cover type. Mean percentage of agriculture land was 86.1% in each buffer zone. The other two primary land cover types identified were forest (2.21%) and water (3.52%) (Figure 3).
Figure 3. Distribution of land cover types in the Canadian prairies.
Figure 3. Distribution of land cover types in the Canadian prairies.
Ijerph 10 03033 g003

3.2. Constructed Models

Based on the results of Pearson correlation and principal component analysis, explanatory variables could be grouped into four highly correlated components (Table 1). The cumulative variances explained by these four components of principal component analysis were 84.6% and 86.1% for Cx. tarsalis abundance and infection rate models, respectively.
In the final Cx. tarsalis abundance model, increased mean monthly temperature, 1 month lagged mean temperature, total precipitation, and temporal lags of precipitation from 1 to 2 months were significantly associated with increased Cx. tarsalis abundance, while an inverse association between forest land cover and Cx. tarsalis abundance was observed (Table 2).
In the final WNV infection rate model, we found that increasing Cx. tarsalis abundance, and 1 month lagged temperature were associated with increased WNV infection rate, while one month lagged mean precipitation, 3 months total precipitation, and water land cover were inverse associated with infection rate (Table 3).
Table 2. Estimated coefficients of explanatory variables in the constructed models of Cx. tarsalis abundance. Single variable indicates the explanatory variables which are assessed individually. Final model represents the final fitted model with the lowest AICc value. Full model is model fitted with all created explanatory variables.
Table 2. Estimated coefficients of explanatory variables in the constructed models of Cx. tarsalis abundance. Single variable indicates the explanatory variables which are assessed individually. Final model represents the final fitted model with the lowest AICc value. Full model is model fitted with all created explanatory variables.
VariablesSingle variableFinal modelFull model
Coef.95% CICoef.95% CICoef.95% CI
Intercept−3.48 *−4.05 to −2.91−3.93 *−4.6 to −3.25
Weather
Monthly mean temperature0.25 *0.22 to 0.260.22 *0.2 to 0.250.22 *0.19 to 0.25
1 month lagged temperature0.08 *0.07 to 0.10.07 *0.05 to 0.090.06 *0.04 to 0.09
Winter mean temperature−0.04 *−0.06 to −0.01−0.03 *−0.06 to −0.01
Monthly mean degree days0.28 *0.23 to 0.33 0.032−0.04 to 0.1
Monthly total precipitation−0.003 *−0.004 to −0.0010.0033 *0.002 to 0.0050.0032 *0.002 to 0.005
1 month lagged precipitation0.006 *0.005 to 0.0070.0042 *0.003 to 0.0050.0037 *0.002 to 0.004
2 month lagged precipitation0.005 *0.004 to 0.0060.0033 *0.002 to 0.0040.003 *0.002 to 0.005
Land cover 1
Forest−0.48 *−0.91 to −0.04−0.54 *−0.9 to −0.17−0.59 *−0.95 to −0.22
Water−0.11 *−0.47 to −0.260.03−0.28 to 0.34
Coef.: estimated variable coefficient; * P < 0.05; 1 agriculture land was used as a reference group.
Table 3. Estimated coefficients of explanatory variables in the constructed models of WNV infection rate. Single variable indicates the explanatory variables which are assessed individually. Final model represents the final fitted model with the lowest AICc value. Full model is model fitted with all created explanatory variables.
Table 3. Estimated coefficients of explanatory variables in the constructed models of WNV infection rate. Single variable indicates the explanatory variables which are assessed individually. Final model represents the final fitted model with the lowest AICc value. Full model is model fitted with all created explanatory variables.
VariablesSingle variableFinal modelFull model
Coef.95% CICoef.95% CICoef.95% CI
Intercept−2.26 *−4.47 to −0.05−1.64−5.64 to 2.37
Cx. tarsalis abundance0.16 *0.03 to 0.300.55 *0.31 to 0.790.58 *0.28 to 0.87
Weather
Monthly mean temperature−0.14 *−0.2 to −0.08−0.04−0.18 to 0.10
1 month lagged temperature0.25 *0.21 to 0.290.32 *0.22 to 0.410.32 *0.21 to 0.42
Winter mean temperature0.23 *0.06 to 0.400.01 *−0.13 to 0.15
3 months total of monthly mean degree days0.20 *0.16 to 0.24−0.10 *−0.2 to −0.01−1.10−0.21 to 0.002
Monthly total precipitation−0.015 *−0.02 to −0.01−0.01−0.02 to 0.003
1 month lagged mean precipitation−0.48 *−0.56 to −0.39−0.27 *−0.36 to −0.18−0.43 *−0.62 to −0.24
2 month lagged total precipitation0.003−0.001 to 0.01−0.01−0.02 to 0.003
3 months total precipitation−0.085−0.11 to 0.06−0.05 *−0.08 to −0.020.013−0.06 to 0.08
Land cover 1
Forest−1.3 *−1.84 to −0.76−0.43−1.27 to 0.41
Water−1.31 *−2.8 to −0.182−1.52 *−2.56 to −0.47−1.61 *−2.85 to −0.38
Coef.: estimated variable coefficient; *: P < 0.05; 1: agriculture land was used as a reference group.
In addition, we found the inverse association between the 3-months total of monthly mean degree days and WNV infection rate when the lagged mean temperature was controlled. Time lagged mean temperature was the distorter variable which distorted the coefficient of the 3-months total of monthly mean degree days from a positive association when this variable was tested alone to negative association.
Second-order polynomial variables of temperature and precipitation were not significantly associated with Cx. tarsalis abundance and infection rate. The RMSE of training and validation dataset was 0.97 and 1.0, respectively. Both RMSE values were close to one and close to each other which indicated the accuracy and precision of model predictability. There was no significant spatial autocorrelation detected for residuals of Cx. tarsalis infection rate by Moran’s I test (p = 0.65).
We used the Cx. tarsalis infection rate of 20 per 1,000 as a criterion to represent the high risk area in our predictive maps based on the mean mosquito infection rate in August 2007, a major epidemic period of WNV in the Canadian prairies. In our predictive maps, the southern part of the Canadian prairies was generally at higher risk, especially in southeast Alberta, southwest to southeast Saskatchewan and southwest Manitoba (Figure 4).
Figure 4. Maps of predicted WNV infection rate in female Cx. tarsalis per 1,000 in July and August 2005–2008 in the Canadian prairies and log transformed human incidence (cases per 100,000 individuals) for each health region in the entire WNV transmission season.
Figure 4. Maps of predicted WNV infection rate in female Cx. tarsalis per 1,000 in July and August 2005–2008 in the Canadian prairies and log transformed human incidence (cases per 100,000 individuals) for each health region in the entire WNV transmission season.
Ijerph 10 03033 g004

4. Discussion

This study analyzed the effects of environmental drivers on Cx. tarsalis abundance and infection rate of WNV from 2005 to 2008. Most of the geographic regions at highest risk based on predicted WNV infection rate in Cx. tarsalis were consistent with the distribution of human WNV cases between 2005 and 2008. The slight differences between the distribution of mosquito infection rate and human incidence were expected, and can potentially be explained by social and economic factors, population density, risk perception, mosquito control programs, and the intensity of human case detection in different health regions. However, mosquito infection rate as an indicator for early season forecasting or predicting of WNV human incidence has been previously demonstrated [9]. Mosquito infection rate was also influenced by variations in climate factors. Of particular interest here, understanding the relationship between climate factors and mosquito infection on monthly scale could be adopted to predict WNV activity under different climate change scenarios [23].
In this study, more explanatory variables were evaluated and compared with our previously published weekly model [22]. The variables associated significantly with WNV infection rate in this study were similar to the weekly model; however, on the monthly time scale, we found an inverse relationship between mean degree days and WNV infection rate when time lagged temperature was controlled. In addition, we demonstrated the effects of land cover composition in a 20 km radius buffer zone on the abundance of Cx. tarsalis and WNV infection rate. We also found a better predictability of monthly model than weekly model based on the values of Root Mean Square Error (data not shown).
In both Cx. tarsalis abundance and WNV infection rate models on the monthly time scale, increasing temperature and lagged temperature significantly increased Cx. tarsalis abundance and infection rate. The mean temperature in May and June, 2007 was higher than 2005 and 2008, with the highest mean temperature in July during study period. High environmental temperature with the antecedent highest precipitation in May and Cx. tarsalis abundance in June might have contributed to the outbreak of WNV in 2007. Increasing environmental temperature shortens the maturation time required for Cx. tarsalis and the extrinsic incubation period of virus. Furthermore, it also accelerates the gonotrophic cycle and affects mosquito survival. In combination, these relationships influence virus transmission by increasing the contact rate between mosquito and host [15,33].
This study found the distorted effect of time lagged mean temperature on the 3 months total of monthly mean degree days for the model of WNV infection rate. Increasing 3 months total of monthly mean degree days can decrease the WNV infection rate in Cx. tarsalis when the variable of time lagged mean temperature was controlled. In addition, the Pearson correlation and principal component analysis revealed that mean temperature fluctuations, and 3 months total of monthly mean degree days were highly correlated. These variables were all estimated based on the maximum and minimum temperature and indicated the degree of temperature fluctuation. These findings indicated that increasing temperature fluctuation could decrease the Cx. tarsalis infection rate in the environment with similar time lagged mean temperature. The effect of high temperature fluctuation has been proposed to limit the midgut infection of flaviviruses in mosquitoes by preventing the virus from entering the midgut epithelial cells or limiting initial replication of virus in these midgut cells [44]. Adverse effect of high temperature fluctuations have also been observed in the transmission of dengue virus by Aedes aegypti [45] and western equine encephalomyelitis by Cx. tarsalis [46,47].
The probable contact rate between an infected competent vector and a susceptible host is essential for maintaining the enzootic cycle of an arbovirus in a geographic area, although herd immunity of host population might dampen the transmission of WNV [48]. A minimum threshold of vector and susceptible host interaction is needed to allow for virus transmission [49]. Therefore, population abundance of suitable competent vectors is usually an indicator of the prevalence of pathogen occurrence [49,50,51]. A positive association has been demonstrated for other arboviruses such as Japanese encephalitis, western equine encephalitis and St. Louis encephalitis [49,52,53]. Furthermore, abundance of Cx. tarsalis, Cx. p. quinquefasciatus, Cx. pipiens and Cx. restuans has also been used as indicators to predict human WNV risk in different regions [54,55,56].
In contrast to temperature, the influence of precipitation on the enzootic cycle of WNV is paradoxical. Increased precipitation is generally believed to create standing water suitable for mosquito breeding and increase mosquito population and the risk of vector borne disease. This reason likely accounts for why precipitation and time lagged precipitation were positively associated with mosquito abundance in our study. In 2007, total precipitation in May was highest compared to other years with the mean temperature closer to 2006 but higher than 2005 and 2008 in the Canadian prairies. The high precipitation might have contributed to the explosive occurrence of Cx. tarsalis in June. Other researchers have demonstrated that preceding droughts can increase the incidence of WNV in human population in the western United States [57] and the Canadian prairies [58]. Potential explanations for the effects of precipitation on WNV incidence include drought induced decreases in mosquito competitors or predators and increase the abundance of Cx. tarsalis, increased congregation of susceptible avian hosts and mosquito vector on dwindling wetland habitat, and changes in the composition of the avian host community [17,59,60]. In our study, precipitation was positively associated with Cx. tarsalis abundance, similar to studies in California and the northern Great Plain habitat in South Dakota [14,16] but an inverse association between precipitation and WNV infection rate in Cx. tarsalis was observed. These findings suggest that other factors, such as the aggregation of competent hosts and Cx. tarsalis [59], or alterations to composition of the avian community might better explain the effects of precipitation on WNV infection rate in Cx. tarsalis in the Canadian prairies [61].
Nonlinear relationships between climate factors and arthropod vectors are commonly observed. For instance, Reisen et al. [14] found that spring Cx. tarsalis population was positively correlated with temperature in winter and spring, whereas summer abundance was correlated negatively with spring temperature and not correlated with summer temperature in California. A unimodal relationship between precipitation and WNV incidence was demonstrated in the northern Great Plains and the optimal total precipitation from May to July was approximately 200 mm [62]. In our study, all the second-order polynomial variables were not associated with Cx. tarsalis abundance or infection rate. Total precipitation in May to July of 2005 to 2008 ranged from 79 to 332 mm and observed high risk areas had the lowest precipitation values. The warmest summer mean temperature was around 20 °C in July (Figure 2), which was lower than the mean temperature of 30 °C for the same month in California [56]. Low summer environmental temperatures across the Canadian prairies could be the main reason for our findings. In this temperature range, the development rate of Cx. tarsalis and the replication rate of WNV in the Cx. tarsalis are linearly related with environmental temperature [13,15].
Ezenwa et al. [61] found a negative association between WNV infection rate among Culex mosquitoes and wetland coverage. Wetland area has been positively associated with bird species diversity and thus it is possible that this represents an example of the dilution effect, in which increased bird diversity led to overall decreases in the mosquito infection rate [61,63]. In our models, Cx. tarsalis abundance was not significantly different between wetland and agriculture as the primary composition of land cover types; however, a significant negative association was found between wetland and Cx. tarsalis infection rate. This result reflects the effect of possible underlying factors, such as the differences of bird species composition between different land cover types and also indicates that vector abundance alone was not sufficient to predict the intensity of WNV occurrence.

5. Conclusions

Our study clarifies the relationship between environmental factors and the abundance of Cx. tarsalis and infection rate of WNV in the Canadian prairies. The observed association between environmental temperature and WNV infection rate could provide sufficient time to predict WNV occurrence and initiate disease control and public health interventions. In addition, warmer temperature and increased variability in precipitation are already being observed and are projected to accelerate in the Canadian prairies [64]. Predictive monthly models for vector-borne diseases are critical tools for public health and wildlife management in a future of rapid climate change. These models could be adopted to assess the effects of changing climate conditions on WNV in the Canadian prairies [23].

Acknowledgements

We thank the Pilot Infectious Disease Impact and Response System (PIDIRS)/program of Public Health Agency of Canada (PHAC) for their funding support, and Environment Canada (climate) and Public Health Division of Manitoba Health (mosquito) for providing useful data.

Conflict of Interest

The authors declare no conflict of interests.

References

  1. Nash, D.; Mostashari, F.; Fine, A.; Miller, J.; O’Leary, D.; Murray, K.; Huang, A.; Rosenberg, A.; Greenberg, A.; Sherman, M.; et al. The outbreak of West Nile virus infection in the New York city area in 1999. N. Engl. J. Med. 2001, 344, 1807–1814. [Google Scholar] [CrossRef]
  2. Public Health Agency of Canada. West Nile MONITOR-2007 Human Surveillance. Available online: http://www.phac-aspc.gc.ca/WNV-VNW/mon-hmnsurv-2007.eng.php (accessed on 23 February 2013).
  3. Curry, P. Saskatchewan mosquitoes and West Nile virus. Blue Jay 2004, 62, 104–111. [Google Scholar]
  4. Turell, M.J.; Dohm, D.J.; Sardelis, M.R.; O’Guinn, M.L.; Andreadis, T.G.; Blow, J.A. An update on the potential of North American mosquitoes (Diptera: Culicidae) to transmit West Nile virus. J. Med. Entomol. 2005, 42, 57–62. [Google Scholar] [CrossRef]
  5. Hayes, E.B.; Komar, N.; Nasci, R.S.; Montgomery, S.P.; O’Leary, D.R.; Campbell, G.L. Epidemiology and transmission dynamics of West Nile virus disease. Emerg. Infect. Dis. 2005, 11, 1167–1173. [Google Scholar] [CrossRef]
  6. Epstein, P.R. West Nile virus and the climate. J. Urban Health 2001, 78, 367–371. [Google Scholar] [CrossRef]
  7. Eldridge, B. Strategies for surveillance, prevention, and control of arbovirus diseases in western North America. Am. J. Trop. Med. Hyg. 1987, 37, 77–86. [Google Scholar]
  8. Gu, W.; Lampman, R.; Krasavin, N.; Berry, R.; Novak, R. Spatio-temporal analyses of West Nile virus transmission in Culex mosquitoes in northern Illinois, USA, 2004. Vector Borne Zoonotic Dis. 2006, 6, 91–98. [Google Scholar] [CrossRef]
  9. Brownstein, J.S.; Holford, T.R.; Fish, D. Enhancing West Nile virus surveillance, United States. Emerg. Infect. Dis. 2004, 10, 1129–1133. [Google Scholar] [CrossRef]
  10. Yiannakoulias, N.W.; Schopflocher, D.P.; Svenson, L.W. Modelling geographic variations in West Nile virus. Can. J. Public Health. 2006, 97, 374–378. [Google Scholar]
  11. Goddard, L.B.; Roth, A.E.; Reisen, W.K.; Scott, T.W. Vertical transmission of West Nile virus by three California Culex (Diptera: Culicidae) species. J. Med. Entomol. 2003, 40, 743–746. [Google Scholar] [CrossRef]
  12. Tempelis, C.H.; Reeves, W.C.; Bellamy, R.E.; Lofy, M.F. A three-year study of the feeding habits of Culex tarsalis in Kern county, California. Am. J. Trop. Med. Hyg. 1965, 14, 170–177. [Google Scholar]
  13. Reisen, W.K. Effect of temperature on Culex tarsalis (Diptera: Culicidae) from the Coachella and San Joaquin valleys of California. J. Med. Entomol. 1995, 32, 636–645. [Google Scholar]
  14. Reisen, W.K.; Cayan, D.; Tyree, M.; Barker, C.M.; Eldridge, B.; Dettinger, M. Impact of climate variation on mosquito abundance in California. J. Vector Ecol. 2008, 33, 89–98. [Google Scholar] [CrossRef]
  15. Reisen, W.K.; Fang, Y.; Martinez, V.M. Effects of temperature on the transmission of West Nile virus by Culex tarsalis (Diptera: Culicidae). J. Med. Entomol. 2006, 43, 309–317. [Google Scholar] [CrossRef]
  16. Chuang, T.W.; Hildreth, M.B.; Vanroekel, D.L.; Wimberly, M.C. Weather and land cover influences on mosquito populations in Sioux Falls, South Dakota. J. Med. Entomol. 2011, 48, 669–679. [Google Scholar] [CrossRef]
  17. Wang, G.; Minnis, R.; Belant, J.; Wax, C. Dry weather induces outbreaks of human West Nile virus infections. BMC Infect. Dis. 2010, 10, 38–44. [Google Scholar] [CrossRef]
  18. Jenkins, D.W. Bionomics of Culex tarsalis in relation to Western equine encephalomyelitis. Am. J. Trop. Med. Hyg. 1950, 30, 909–916. [Google Scholar]
  19. Hongoh, V.; Berrang-Ford, L.; Scott, M.; Lindsay, L. Expanding geographical distribution of the mosquito, Culex pipiens, in Canada under climate change. Appl. Geogr. 2012, 3, 53–62. [Google Scholar]
  20. Lafferty, K.D. The ecology of climate change and infectious diseases. Ecology 2009, 90, 888–900. [Google Scholar] [CrossRef]
  21. Brooker, S.; Hay, S.I.; Bundy, D.A.P. Tools from ecology: Useful for evaluating infection risk models? Trends Parasitol. 2002, 18, 70–74. [Google Scholar] [CrossRef]
  22. Chen, C.-C.; Epp, T.; Jenkins, E.; Waldner, C.; Curry, P.S.; Soos, C. Predicting weekly variation of Culex tarsalis (Diptera: Culicidae) West Nile virus infection in a newly endemic region, the Canadian prairies. J. Med. Entomol. 2012, 49, 1144–1153. [Google Scholar] [CrossRef]
  23. Chen, C.C.; Jenkins, E.; Epp, T.; Waldner, C.; Curry, P.S.; Soos, C. Climate change and West Nilevirus in a highly endemic region of North America. Int. J. Environ. Res. Public Health 2013, 3052–3071. [Google Scholar]
  24. Lanciotti, R.S.; Kerst, A.J.; Nasci, R.S.; Godsey, M.S.; Mitchell, C.J.; Savage, H.M.; Komar, N.; Panella, N.A.; Allen, B.C.; Volpe, K.E.; et al. Rapid detection of West Nile virus from human clinical specimens, field-collected mosquitoes, and avian samples by a TaqMan reverse transcriptase-PCR assay. J. Clin. Microbiol. 2000, 38, 4066–4071. [Google Scholar]
  25. Drebot, M.A.; Lindsay, R.; Barker, I.K.; Buck, P.A.; Fearon, M.; Hunter, F.; Sockett, P.; Artsob, H. West Nile virus surveillance and diagnostics: A Canadian perspective. Can. J. Infect. Dis. 2003, 14, 105–114. [Google Scholar]
  26. Biggerstaff, B. PooledInfRate, Version 3.0: A Microsoft Excel Add-In to Compute Prevalence Estimates from Pooled Samples; Centers for Disease Control and Prevention: Fort Collins, CO, USA, 2006. [Google Scholar]
  27. Chiang, C.L.; Reeves, W.C. Statistical estimation of virus infection rates in mosquito vector populations. Am. J. Epidemiol. 1962, 75, 377–391. [Google Scholar]
  28. Bernard, K.A.; Maffei, J.G.; Jones, S.A.; Kauffman, E.B.; Ebel, G.; Dupuis, A. West Nile virus infection in birds and mosquitoes, New York State, 2000. Emerg. Infect. Dis. 2001, 7, 679–685. [Google Scholar]
  29. Gu, W.; Lampman, R.; Novak, R.J. Assessment of arbovirus vector infection rates using variable size pooling. Med. Vet. Entomol. 2004, 18, 200–204. [Google Scholar] [CrossRef]
  30. Gu, W.; Novak, R.J. Short report: Detection probability of arbovirus infection in mosquito populations. Am. J. Trop. Med. Hyg. 2004, 71, 636–638. [Google Scholar]
  31. Huffman, T.; Ogston, R.; Fisette, T.; Daneshfar, B.; Gasser, P.; White, L.; Maloley, M.; Chenier, R. Canadian agricultural land-use and land management data for Kyoto reporting. Can. J. Soil Sci. 2006, 86, 431–439. [Google Scholar] [CrossRef]
  32. Reisen, W.K.; Lothrop, H.D. Population ecology and dispersal of Culex tarsalis (Diptera: Culicidae) in the Coachella Valley of California. J. Med. Entomol. 1995, 32, 490–502. [Google Scholar]
  33. Zou, L.; Miller, S.N.; Schmidtmann, E.T. A GIS tool to estimate West Nile virus risk based on a degree-day model. Environ. Monit. Assess. 2007, 129, 413–420. [Google Scholar] [CrossRef]
  34. Allen, J. A modified sine wave method for calculating degree days. Environ. Entomol. 1976, 5, 388–396. [Google Scholar]
  35. Schrag, A.; Konrad, S.; Miller, S.; Walker, B.; Forrest, S. Climate-change impacts on sagebrush habitat and West Nile virus transmission risk and conservation implications for greater sage-grouse. GeoJournal 2011, 76, 561–575. [Google Scholar] [CrossRef]
  36. Dohoo, I.; Martin, W.; Stryhn, H. Model-Building Strategies. In Veterinary Epidemiologic Research, 2nd ed.; VER Inc.: Charlottetown, PEI, Canada, 2009; pp. 365–390. [Google Scholar]
  37. Guadagnoli, E.; Velicer, W.F. Relation of sample size to the stability of component patterns. Psychol. Bull. 1988, 103, 265–275. [Google Scholar] [CrossRef]
  38. Bolker, B.; Brooks, M.; Clark, C.; Geange, S.; Poulsen, J.; Stevens, M.; White, J. Generalized linear mixed models: A practical guide for ecology and evolution. Trends Ecol. Evol. 2009, 24, 127–135. [Google Scholar] [CrossRef]
  39. Dohoo, I.; Martin, W.; Stryhn, H. Mixed Models for Continuous Data. In Veterinary Epidemiologic Research, 2nd ed.; VER Inc.: Charlottetown, PEI, Canada, 2009; pp. 553–578. [Google Scholar]
  40. Burnham, K.P.; Anderson, D.R. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach; Springer Verlag: New York, NY, USA, 2002; pp. 49–96. [Google Scholar]
  41. Government of Alberta-Health and Wellness, West Nile Virus. Available online: http://www.health.alberta.ca/health-info/WNv-evidence.html (accessed on 8 September 2012).
  42. Government of Manitoba-Manitoba Health, West Nile Virus. Available online: http://www.gov.mb.ca/health/wnv/stats.html (accessed on 8 September 2012).
  43. Government of Saskatchewan-Public Health, West Nile Virus. Available online: http://www.health.gov.sk.ca/wnv-surveillance-results-archive (accessed on 8 September 2012).
  44. Kramer, L.D.; Ebel, G.D. Dynamics of flavivirus infection in mosquitoes. Adv. Virus Res. 2003, 60, 187–232. [Google Scholar] [CrossRef]
  45. Lambrechts, L.; Paaijmans, K.P.; Fansiri, T.; Carrington, L.B.; Kramer, L.D.; Thomas, M.B.; Scott, T.W. Impact of daily temperature fluctuations on Dengue virus transmission by Aedes aegypti. Proc. Natl. Acad. Sci. USA 2011, 108, 7460–7465. [Google Scholar] [CrossRef]
  46. Kramer, L.; Hardy, J.; Presser, S. Effect of temperature of extrinsic incubation on the vector competence of Culex tarsalis for western equine encephalomyelitis virus. Am. J. Trop. Med. Hyg. 1983, 32, 1130–1139. [Google Scholar]
  47. Reisen, W.; Meyer, R.; Presser, S.; Hardy, J. Effect of temperature on the transmission of western equine encephalomyelitis and St. Louis encephalitis viruses by Culex tarsalis (Diptera: Culicidae). J. Med. Entomol. 1993, 30, 151–160. [Google Scholar]
  48. Kwan, J.L.; Kluh, S.; Reisen, W.K. Antecedent avian immunity limits tangential transmission of West Nile virus to humans. PLoS One 2012, 7, e34127. [Google Scholar] [CrossRef]
  49. Reeves, W.C. Ecology of mosquitoes in relation to arboviruses. Annu. Rev. Entomol. 1965, 10, 25–46. [Google Scholar] [CrossRef]
  50. Murray, M.D. Influences of vector biology on transmission of arboviruses and outbreaks of disease: The Culicoides brevitarsis model. Vet. Microbiol. 1995, 46, 91–99. [Google Scholar] [CrossRef]
  51. Saugstad, E.S.; Dalrymple, J.M.; Eldridge, B.F. Ecology of arboviruses in a Maryland freshwater swamp: Population dynamics and hibitat distribution of potential mosquito vectors. Am. J. Epidemiol. 1972, 96, 114–122. [Google Scholar]
  52. Arunachalam, N.; Murty, U.S.N.; Narahari, D.; Balasubramanian, A.; Samuel, P.P.; Thenmozhi, V.; Paramasivan, R.; Rajendran, R.; Tyagi, B.K. Longitudinal studies of japanese encephalitis virus infection in vector mosquitoes in Kurnool district, Andhra Pradesh, South India. J. Med. Entomol. 2009, 46, 633–639. [Google Scholar] [CrossRef]
  53. Wegbreit, J.; Reisen, W.K. Relationships among weather, mosquito abundance, and encephalitis virus activity in California: Kern County 1990–1998. J. Am. Mosq. Control Assoc. 2000, 16, 22–27. [Google Scholar]
  54. Kilpatrick, A.M.; Kramer, L.D.; Campbell, S.R.; Alleyne, E.O.; Dobson, A.P.; Daszak, P. West Nile virus risk assessment and the bridge vector paradigm. Emerg. Infect. Dis. 2005, 11, 425–429. [Google Scholar] [CrossRef]
  55. Kwan, J.L.; Kluh, S.; Madon, M.B.; Reisen, W.K. West Nile virus emergence and persistence in Los Angeles, California, 2003–2008. Am. J. Trop. Med. Hyg. 2010, 83, 400–412. [Google Scholar] [CrossRef]
  56. Reisen, W.K.; Carroll, B.D.; Takahashi, R.; Fang, Y.; Garcia, S.; Martinez, V.M.; Quiring, R. Repeated West Nile virus epidemic transmission in Kern County, California, 2004–2007. J. Med. Entomol. 2009, 46, 139–157. [Google Scholar] [CrossRef]
  57. Landesman, W.J.; Allan, B.F.; Langerhans, R.B.; Knight, T.M.; Chase, J.M. Inter-annual associations between precipitation and human incidence of West Nile virus in the United States. Vector Borne Zoonotic Dis. 2007, 7, 337–343. [Google Scholar] [CrossRef]
  58. Epp, T.Y.; Waldner, C.L.; Berke, O. Predicting geographical human risk of West Nile virus-Saskatchewan, 2003 and 2007. Can. J. Public Health. 2009, 100, 344–348. [Google Scholar]
  59. Shaman, J.; Day, J.F.; Stieglitz, M. Drought-induced amplification and epidemic transmission of West Nile virus in southern Florida. J. Med. Entomol. 2005, 42, 134–141. [Google Scholar] [CrossRef]
  60. Chase, J.M.; Knight, T.M. Drought-induced mosquito outbreaks in wetlands. Ecol. Lett. 2003, 6, 1017–1024. [Google Scholar] [CrossRef]
  61. Ezenwa, V.O.; Milheim, L.E.; Coffey, M.F.; Godsey, M.S.; King, R.J.; Guptill, S.C. Land cover variation and West Nile virus prevalence: Patterns, processes, and implications for disease control. Vector Borne Zoonotic Dis. 2007, 7, 173–180. [Google Scholar] [CrossRef]
  62. Wimberly, M.C.; Hildreth, M.B.; Boyte, S.P.; Lindquist, E.; Kightlinger, L. Ecological niche of the 2003 West Nile virus epidemic in the northern Great Plains of the United States. PLoS One 2008, 3, e3744. [Google Scholar] [CrossRef]
  63. Ezenwa, V.O.; Godsey, M.S.; King, R.J.; Guptill, S.C. Avian diversity and West Nile virus: Testing associations between biodiversity and infectious disease risk. Proc. R. Soc. Biol. Sci. Ser. 2006, 273, 109–117. [Google Scholar] [CrossRef]
  64. Shepherd, A.; McGinn, S. Assessment of climate change on the Canadian prairies from downscaled GCM data. Atmosphere-Ocean 2003, 41, 301–316. [Google Scholar] [CrossRef]

Share and Cite

MDPI and ACS Style

Chen, C.-C.; Epp, T.; Jenkins, E.; Waldner, C.; Curry, P.S.; Soos, C. Modeling Monthly Variation of Culex tarsalis (Diptera: Culicidae) Abundance and West Nile Virus Infection Rate in the Canadian Prairies. Int. J. Environ. Res. Public Health 2013, 10, 3033-3051. https://doi.org/10.3390/ijerph10073033

AMA Style

Chen C-C, Epp T, Jenkins E, Waldner C, Curry PS, Soos C. Modeling Monthly Variation of Culex tarsalis (Diptera: Culicidae) Abundance and West Nile Virus Infection Rate in the Canadian Prairies. International Journal of Environmental Research and Public Health. 2013; 10(7):3033-3051. https://doi.org/10.3390/ijerph10073033

Chicago/Turabian Style

Chen, Chen-Chih, Tasha Epp, Emily Jenkins, Cheryl Waldner, Philip S. Curry, and Catherine Soos. 2013. "Modeling Monthly Variation of Culex tarsalis (Diptera: Culicidae) Abundance and West Nile Virus Infection Rate in the Canadian Prairies" International Journal of Environmental Research and Public Health 10, no. 7: 3033-3051. https://doi.org/10.3390/ijerph10073033

APA Style

Chen, C. -C., Epp, T., Jenkins, E., Waldner, C., Curry, P. S., & Soos, C. (2013). Modeling Monthly Variation of Culex tarsalis (Diptera: Culicidae) Abundance and West Nile Virus Infection Rate in the Canadian Prairies. International Journal of Environmental Research and Public Health, 10(7), 3033-3051. https://doi.org/10.3390/ijerph10073033

Article Metrics

Back to TopTop