key: cord-0295350-dlehspaq
authors: Springborn, M. R.; Weill, J. A.; Lips, K. R.; Ibanez, R.; Ghosh, A.
title: Amphibian Collapses Exacerbated Malaria Outbreaks in Central America
date: 2020-12-09
journal: nan
DOI: 10.1101/2020.12.07.20245613
sha: d93dfbcb6fc0d566b4f937d464d0af61ac2feecd
doc_id: 295350
cord_uid: dlehspaq

Ecosystems play an important role in supporting human welfare, including regulating the transmission of infectious diseases. Many of these services are not fully-appreciated due to complex environmental dynamics and lack of baseline data. Multicontinental amphibian decline due to the fungal pathogen Batrachochytrium dendrobatidis (Bd) provides a stark example. Even though amphibians are known to affect natural food webs---including mosquitoes that transmit human diseases---the human health impacts connected to their massive decline have received little attention. Here we show a causal link between a wave of Bd-driven collapse of amphibians in Central America and increased human malaria incidence. At the canton-level in Costa Rica and district-level in Panama, expected malaria incidence increased for eight years subsequent to amphibian losses, peaking at an additional 1.0 cases per 1,000 population (CPK). The increase is substantial in comparison to annual incidence levels from outbreaks in these countries, which peaked at 1.1-1.5 CPK during our period of study from 1976-2016. This pattern holds across multiple alternative approaches to the estimation model. This previously unidentified impact of biodiversity loss illustrates the often hidden human welfare costs of conservation failures. These findings also show the importance of mitigating international trade-driven spread of similar emergent pathogens like Batrachochytrium salamandrivorans.

two countries are consistent with a lagged impact of amphibian decline, this correlation does 59 not establish causality on its own. (13, 14) . Deforestation in particular has received increased attention in recent years and is hy-65 pothesized to operate through changes to the physical environment, malarial mosquito biology, 66 and human exposure (15). While most have found that deforestation is associated with increased 67 malaria incidence (16), this result does not hold across all regions and study designs (17). How-68 ever, linkages between malarial dynamics and ecosystem disruption by invasive species has not 69 been previously studied, aside from well-known linkages to invasive mosquito vectors. 70 Methods 71 We use a multiple regression model to estimate the causal impact of Bd-driven amphibian 72 decline on malaria incidence at the canton level in Costa Rica and distrito level in Panama, 73 5 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint hereafter referred to as "county level". We exploit variation in outcomes for units (counties) 74 that experienced "treatment" (amphibian decline) at different times, which is a difference-in- 75 difference, event-study design (18). This approach takes advantage of the staggered treatment 76 of counties, leveraging differences in malaria outcomes over time between administrative units 77 that have and have not been treated with Bd. For the regression model we use an event-study framework to characterize the impact of 90 Bd-driven amphibian decline on per capita incidence of malaria, while controlling for other 91 potential drivers. This approach is standard in the econometric literature in cases where multiple 92 units, such as states or counties, receive the same "treatment" at different times (19-21). We 93 specified the model as to the event is K c,1992 = 1992 − 1990 = 2; county c has completed two years of treatment and is 99 entering its third. Our main regressor of interest is τ ck , which is a dummy variable equal to one 100 if county c is k years away from the initial treatment event: τ ck = 1{K ct = k}. We focus on 101 the five relative years before the event and 10 years following: k ∈ {−5, −4, ...9, 10}. We also 102 included a single dummy for all relative years before this, denoted by k = −6 − and another for 103 all relative years after, k = 11 + . Allowing the coefficient γ k to vary for each relative year (K ct ) 104 in this way facilitates flexible and dynamic treatment effects. Because we imposed γ −1 = 0 to study framework is to confirm the absence of a pre-trend, i.e. a directional trend in the effect 132 8 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint on malaria cases of the year relative to amphibian decline before treatment at year 0. Simply 133 put, we would not expect to see systematic movement in malaria incidence in years before 134 amphibian decline begins. We confirmed lack of such a pre-trend: none of the coefficients for 135 k < 0 are significantly different from zero. 136 Overall, we estimated a significant increase in malaria cases due to the onset of amphibian 137 decline, an effect that starts gradually, plateaus after 3 years, and starts to attenuate after 8 years.

The first year of amphibian decline (k = 0) reflects partial treatment for most counties since Bd-139 saturation of a county takes a median of 1.1 years and spread may arrive anytime in a calendar 140 year. Starting in year k = 1, amphibian decline is associated with a statistically significant 141 increase in malaria cases. We estimate that this average effect reaches a relative plateau by year 142 k = 3 and stays relatively constant for six years. For one year in this range (k = 6) the effect is 143 not significantly different from zero. This is not due to a decline in the coefficient but rather to 144 an increase in the standard error due to an increase in residuals, i.e., additional noise. Starting In Table 1 we present the full set of regression estimates for the preferred specification (dis-154 cussed above) in column 1, alongside estimates for three alternative specifications (columns 155 2-4) to check robustness (discussed further below). In the table, regression coefficients are pre-156 9 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 9, 2020. ; Fig. 3 . Estimated effect on malaria cases per 1,000 population (vertical axis) of year k (horizontal axis) relative to Bd-driven date of amphibian decline (DoD). Shading represents 90% confidence intervals.

sented with standard errors in parentheses, which are clustered at the county level. We clustered 157 due to the sampling design (we are inferring something about the larger population based on 158 data sampled at the county level) and our quasi-experimental design ("treatment" occurs at the 159 county level). Overall, we found that our key qualitative results discussed above hold across an 160 array of robustness checks.

In the table, our independent variables (rows) start with two ground cover measures and 162 three weather measures. Next are the key coefficients of interest,γ k , representing the estimated 163 effect of relative years before a county's DoD (k < 0) and after (k ≥ 0). These coefficients 164 for our preferred specification for k = −5, −4, ..., 10 are plotted in Fig. 2 and discussed above. 165 We excluded k = −1 so that the rest of these coefficients are interpreted as effects relative to 166 this year just before a county's DoD. For our first alternative specification in column (2) we 167 augmented the data set with regions of Panama excluded in our preferred specification due to 168 data limitations as described in the SI Appendix. ii In column (3) we considered an alternative 169 10 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint using the minimum date reflecting initial arrival to the county border, we considered the average 171 DoD for the county. In column (4) we conducted weighted least squares regression where the 172 weights were given by county-level population. This heightened emphasis on observations 173 from high-population counties is motivated by the conjecture that malaria incidence measures 174 from such counties are less subject to measurement error given the larger number of "samples" year appears to stem from an uptick in noise-the coefficient is within the range of surrounding 183 years while the standard error is elevated. 184 We found that the same general pattern holds across all specifications: we fail to find a pre- In results presented in the SI Appendix, we also examined the sensitivity of estimates to an 194 11 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; Table 1 . Estimates for the regression model specified in Equation for the preferred specification (column 1) and alternatives for robustness checks. Standard errors (clustered at the county level) are presented in parentheses. 12 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint alternative method used for estimating the gridded DoD values, specifically a thin-plate spline 195 (TPS) method for direct interpolation of the DoD data over space. Although this method ignores 196 the implications for spread of Central America's irregular coastline, the results are also broadly 197 consistent with our preferred specification.

Additional Drivers of Malaria 199 We also found that decreasing tree cover is associated with a statistically significant increase in Overall we provide novel causal evidence that pathogen-driven amphibian decline can play a 213 significant role in increasing incidence of vector-borne disease. Our results also contribute to CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint new calamities, like international spread of an emergent and closely related pathogen Batra-218 chochytrium salamandrivorans through incompletely regulated live species trade (4).

While results were robust to several alternative specifications, we were unable to examine 220 whether outcomes for malaria held for other diseases. Showing that our results held for other 221 vector-borne diseases (e.g. dengue and leishmaniasis) would have provided additional support 222 for the mechanism we propose. Showing that our results failed to hold for non-vector-borne 223 illnesses like influenza would have provided additional support for the argument that the effect 224 we identify is specific to vector-borne diseases and not a general disease effect. We attempted 225 to obtain these disease data sets from the national ministries of health in both countries but they 226 were not available for our period of study at the county-level needed.

From the data we were able to obtain from the Panamanian Ministry of Health, we found that CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint available at the county level and such investment is endogenous with malaria cases (the outcome 244 variable of interest). CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint

The authors declare no competing interests.

after widespread amphibian loss. Science 367, 814-816 (2020). CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. 18 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint

Below we describe data sources and methods for constructing our panel data set spanning each 330 year from 1976-2016 at the canton-level in Costa Rica and distrito-level in Panama, hereafter 331 referred to as "county-level". This core data set includes variables presented in Table S1 as well 332 as the Bd-driven date of amphibian decline (DoD) over space as presented in Fig. 2 . Because 333 Panama experienced multiple distrito boundary changes during our study period, we aggregated 334 certain distritos for spatial consistency. We omitted a small set of distritos and regions that 335 present aggregation problems or are islands. These steps are described in detail further below. 1 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint 366 relative split between short vegetation and bare ground is the same in each county as in 1982. 

To provide a sense of the magnitude of the impact of deforestation on malaria incidence, 376 we first identified a plausible change in the tree cover variable. We followed the approach of 377 (4) and residualized the tree cover variable with respect to both county and year fixed effects. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. Because evaluation of the residuals resulting from a given set of parameter estimates is not 410 a direct function calculation, we used an iterative algorithm for optimization. We started with 411 a guess for the parameter vector of spread rates, which is informed by the observed physical 412 and temporal distances between data points. iv Then, in each iteration of the algorithm, one of 413 4 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint the parameters was randomly selected to be perturbed, both up and down, by a fixed adjustment 414 factor (e.g. 5%). Then, from the three candidate models (unperturbed, selected parameter ad-415 justed up, selected parameter adjusted down) the one with the lowest SSR is retained and used 416 as the basis for the next iteration. When improvements in the SSR cease, the adjustment factor is 417 halved (e.g. from 5% to 2.5%) to further refine the solution. We performed 200 iterations of the 418 algorithm since after 175 iterations, no further improvements in SSR were typically achieved.

To guard against this stochastic algorithm identifying a local and not a global optimum, the 420 best-performing model is selected from the set resulting from running the algorithm 100 times. 

where x i1 and x i2 are the longitude and latitude coordinates at point i, and y i is the observed 435 5 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint allowed and the optimization problem boils down to ordinary least squares.

In practice, λ can be fixed at a given value or selected through a cross-validation method. 439 We selected λ using generalized cross validation and then further assessed the robustness of our 

The TPS approach is appealing in its simplicity as a one-step fitting procedure that directly 445 extrapolates observed DoD values over space. However, the pattern in Fig. S1 illustrates a key 446 drawback of this method: it ignores realities of irregular landscape shapes, in essence allowing 447 for spread unimpeded across ocean, equivalently to land. 6 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint Our "county-level" spatial unit of analysis corresponds to cantons in Costa Rica and distritos CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint 7 times the median duration, we exclude it from the preferred specification.

The second region we omit from our preferred specification encompasses eastern Panama is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted December 9, 2020. ; Fig. S2 . Estimated effect on malaria cases per 1,000 population (vertical axis) of the k-th year relative to the DoD or Bd-driven date of amphibian decline (horizontal axis) where DoD is estimated using a thinplate spline (TPS). Lines depict results under a range of TPS interpolations, which vary in the penalty parameter, λ. Shading represents 90% confidence intervals.

In main text Fig. 1 we show that during our period of study malaria cases spiked in Costa In Costa Rica, relative to the 1980s, we see a substantial uptick in total funding (govern-504 ment and external) starting 1993 followed by massive increases in 1998 and again in 2008.

For Costa Rica, these funding pulses are correlated first with the onset of the spike in cases 506 and then eventual decline. While spraying houses for mosquitoes increases somewhat in the 507 9 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint iii To construct these variables, Hijmans (5) obtained these monthly climate variables at a 0.5 • spatial resolution 520 from the Climate Research Unit data website (12, 13) and downscaled to the finer spatial resolution using the delta 521 method (14). 522 iv An initial estimate of the rate of spread (meters/week) at each of the DoD data points is calculated as follows. the physical distance. The observed rate of spread from from i − 1 to i, or "approach rate" is given by d i /t i . 526 We specified the rate of spread at point i as a weighted combination of the approach rate and departure rate, CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint nian amphibian species, habitats and elevations during epizootic and enzootic stages. Dis- 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted December 9, 2020. ; https://doi.org/10.1101/2020.12.07.20245613 doi: medRxiv preprint

Global land change from 1982 to 2016

Forest area in costa rica: a comparative study of tropical 541 forest cover estimates over time

Rates of deforestation in the humid tropics: estimates and measurements

Improving the Interpretation of Fixed Effects Regression Results

High resolution monthly climate variables

Riding the wave: Reconciling the 550 roles of disease and climate change in amphibian declines

The amphibian chytrid fungus along an altitudinal 552 transect before the first reported declines in costa rica

Batrachochytrium dendrobatidis infection patterns among panama