Evidence for family-level variation of phenotypic traits in response to temperature of Brazilian Nyssorhynchus darlingi
Parasites & Vectors volume 13, Article number: 55 (2020)
Nyssorhynchus darlingi (also known as Anopheles darlingi) is the primary malaria vector in the Amazon River Basin. In Brazil, analysis of single nucleotide polymorphisms (SNPs) previously detected three major population clusters, and a common garden experiment in a laboratory setting revealed significant population variation in life history traits. Increasing temperatures and local level variation can affect life history traits, i.e. adult longevity, that alter vectorial capacity with implications for malaria transmission in Ny. darlingi.
We investigated the population structure of Ny. darlingi from 7 localities across Brazil utilizing SNPs and compared them to a comprehensive Ny. darlingi catalog. To test the effects of local level variation on life history traits, we reared F1 progeny from the 7 localities at three constant temperatures (20, 24 and 28 °C), measuring key life history traits (larval development, food-starved adult lifespan, adult size and daily survival).
Using nextRAD genotyping-by-sequencing, 93 of the field-collected Ny. darlingi were genotyped at 33,759 loci. Results revealed three populations (K = 3), congruent with major biomes (Amazonia, Cerrado and Mata Atlântica), with greater FST values between biomes than within. In the life history experiments, increasing temperature reduced larval development time, adult lifespan, and wing length in all localities. The variation of family responses for all traits within four localities of the Amazonia biome was significant (ANOVA, P < 0.05). Individual families within localities revealed a range of responses as temperature increased, for larval development, adult lifespan, wing length and survival time.
SNP analysis of several Brazilian localities provided results in support of a previous study wherein populations of Ny. darlingi were clustered by three major Brazilian biomes. Our laboratory results of temperature effects demonstrated that population variation in life history traits of Ny. darlingi exists at the local level, supporting previous research demonstrating the high plasticity of this species. Understanding this plasticity and inherent variation between families of Ny. darlingi at the local level should be considered when deploying intervention strategies and may improve the likelihood of successful malaria elimination in South America.
Malaria has made a comeback in Latin America in the last few years despite a recent period of decline from 2000 to 2014 [1, 2]. The Americas is the only region to have an increase in malaria mortality in 2017 compared to 2010, with a greater number of malaria cases reported in Venezuela, Brazil, and Nicaragua during this period . The main vector and driver of this disease in South America is Nysorrhynchus darlingi (also known as Anopheles darlingi ), that exhibits significant geographical variation in behavior [4, 5] and in phenotypic plasticity [6,7,8]. This species has a natural infection rate by Plasmodium of up to 20% [5, 9], though a more common rate is 1–5% [10,11,12]. Adult vector traits relevant to disease transmission, such as adult life span and body size , can vary between populations [6, 7] and are influenced by larval development conditions such as food quantity  and temperature . Globally, temperatures are projected to rise between 1–4 °C due to climate change . Whereas even small changes in temperature may reduce vectorial capacity , the effects of temperature are not uniform across Ny. darlingi populations . To be successful, future interventions in this region require a better understanding of this vector in a changing environment. Here we assess levels of genetic and phenotypic differentiation among Ny. darlingi populations from Brazil.
The geographical distribution of Nyssorhynchus darlingi includes diverse South American biomes  and is associated with a range of larval habitat types, including natural breeding sites with clean, shaded water and aquatic vegetation near human dwellings , as well as anthropogenic habitats, such as fish ponds  and dams . Habitat modification, e.g. deforestation, was linked to Ny. darlingi presence in Peruvian breeding habitats  and has been positively correlated with malaria cases in Brazil . A mathematical model using field-collected data found that the high biting rate and susceptibility to Plasmodium of Ny. darlingi in the Brazilian Amazon led to a high basic reproductive rate (R0) of malaria (mainly caused by Plasmodium vivax) . The heterogeneity in distribution, vector competence and vectorial capacity of Ny. darlingi presents a major challenge to malaria elimination.
Research on the effects of juvenile stages on adult traits has increased the understanding of developmental trade-offs. Changes in life history traits, such as body size and adult survival, can modify vectorial capacity  and directly impact malaria transmission. In a theoretical climate risk model, the inclusion of the effects of temperature during the full life-cycle, such as juvenile development rate and mortality, revealed that mosquito populations are more sensitive to changes in temperature than adult data alone would indicate . A study of full-sib F1 progeny from field collected An. coluzzii (formerly known as An. gambiae M form , Burkina Faso) found adult longevity increased with adult body size but decreased with longer larval development . Population differentiation for both larval and adult life history traits of Ny. darlingi has been reported at the regional level in Brazil [6, 7], but variation at smaller scales is unstudied.
Average global temperatures are projected to increase 1–4 °C over the next 100 years because of climate change  and tropical insect populations are anticipated to be more negatively affected compared to those in temperate regions . Exotherm development is very sensitive to temperature, which can affect traits relevant to disease transmission, such as body size and adult fitness [15, 25]. Laboratory rearing of An. gambiae suggested an upper thermal limit of 31 °C and complete larvae mortality at 35 °C, with increasing temperatures reducing adult body size and egg production [30, 31]. A malaria model predicted optimal transmission at 25 °C and was validated by an independent malaria transmission data set for An. gambiae (s.l.) and P. falciparum . Modeled parasite and mosquito development rate reached a peak at 30 °C, in contrast, vector competence and vector survivorship peaked at 25 °C.
The analysis of single nucleotide polymorphisms (SNPs) in Anopheles has shed light on population structure [33, 34] and phenotypes [35, 36]. However, results of tests using SNP data to identify population structure of Ny. darlingi in South America have been mixed. Analysis of Ny. darlingi collected from 12 states across Brazil detected three genetic clusters  associated with major biogeographical regions. In contrast, analysis of specimens from three sites within a single biome (Amazonia), between 60–700 km apart, detected significant population divergence at a regional scale , although a subsequent analysis of two of these sites (60 km apart, new vs old settlement) to test for local differentiation in biting behavior found no significant genetic variation . Despite similar methods used in these studies, comparisons between datasets are difficult given the variation in identified loci.
The aim of this study was to investigate local level variation in population structure and in life history traits of Ny. darlingi using a common garden experiment approach to address the following questions: (i) What is the scale of genetic differentiation among populations of Ny darlingi? and (ii) Is there evidence of small-scale variation in life history traits and in plastic responses to temperature variation? Our research combined broad-scale population genetic assays with empirical data from a common garden experiment. We investigated the effects of variation in rearing temperature on a major Neotropical vector, allowing us to assess the extent to which population differences in life history traits were due to environment (temperature), genetics, or both.
Here, we identified molecular genetic variation across biomes, significant phenotypic and genetic variation in life history traits, as well as within-population genetic variation for plasticity of Ny. darlingi. This variation could help tailor current intervention efforts, such as long-lasting insecticide nets (LLINs), indoor residual spraying (IRS) and larval source management (LSM), to regional and local scales for maximum efficiency and malaria elimination.
In this study, we first analyzed the population structure of Brazilian Ny. darlingi with mosquitoes collected within the same year, over a range of seven locality sites. In order to increase our chances of identifying fine-scale population genetic structure, we created a catalog incorporating sequences of Ny. darlingi from Peru and Brazil and used this to reexamine the population structure of Brazilian Ny. darlingi and to investigate the possibility of fine-scale differentiation within three Brazilian biomes. We then investigated the among- and within-locality variation in life history traits such as larval development time, adult lifespan and body size, among these seven localities across Brazil. A common garden experiment of the mosquito populations from the seven localities was conducted in three constant temperature environmental chambers. Mosquitoes were observed from egg hatch to adult death and life history traits recorded.
Study area and field collections
Adult female Ny. darlingi mosquitoes were collected from 7 localities across Brazil (Table 1, Fig. 1), spanning four states and 3 biomes. Details of the collection site criteria for paired sites are found in . Mosquitoes were collected in the evening for 5 hours (17:00–22:00 h) using barrier screens as described in Moreno et al.  for 1–5 days, depending on locality and successful collection of the target species, Ny. darlingi. Blood-fed female mosquitoes from barrier screens were morphologically identified as Ny. darlingi  and maintained individually in a humid box and provided ad libitum sucrose solution during transport to the laboratory in São Paulo, Brazil (Laboratório de Entomologia de Saúde Pública – Culicidae, Faculdade de Saúde Pública, Universidade de São Paulo).
The laboratory rearing for progeny of field-caught individuals was carried out as previously described . Briefly, eggs laid by individual females from each locality, referred to as families, were allowed to hatch and larvae were equally divided (n = 15) into each of three constant temperature environmental chambers (20, 24 and 28 ± 1 °C) (Additional file 1: Table S1) with a 12:12 h light:dark cycle and a relative humidity of 70–80%. Larvae were fed ad libitum and water was changed every other day; adults were provided only water after emergence. Larval, pupal and adult development was assessed daily. Mosquito specimens were maintained in these chambers until natural adult death, whereupon the left wing was collected for body size estimation.
DNA extraction and nextRAD preparation
A subset (n = 93) of the field-collected Ny. darlingi (used to create the families in the life history experiment) was genotyped to evaluate population structure (Table 1). Individuals were selected from each of the 7 localities (n = 12–14 per locality) based on (i) successful egg laying (with priority given to the dams of the families used in the life history research); (ii) complete wing data; and (iii) a DNA concentration between 2.87 and 16.2 mg/ml. Genomic DNA was extracted from all specimens using Qiagen DNeasy Blood and Tissue kit (Qiagen, Germantown, MD, USA) and concentrations were quantified with a Qubit Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). Individuals were sequenced using nextRAD genotyping-by-sequencing methods as described in Emerson et al.  (SNPSaurus, LLC, Eugene, OR, USA). Briefly, the genomic DNA was first fragmented using a Nextera reaction to ligate adapter sequences to the fragments. The fragments were then amplified with an 8 bp Nextera primer (5′-TGC AGG AG-3′), and the library was pooled and purified, with size selection between 350–500 bp. The resulting library was then sequenced, generating 150 bp reads on two lanes of an Illumina HiSeq 4000.
Nyssorhynchus darlingi catalog creation
All raw sequences were analyzed using STACKS v2.3b . Sequences of 24 representative field collected female Ny. darlingi were used to create a catalog using STACKS cstacks, permitting 4 mismatches between stacks and enabling gapped alignments. In order to generalize the catalog to be useful across projects, samples used in this catalog were from previous publications [37, 41, 42]; similar sequencing methods were employed. This catalog consisted of 13 individuals from this study (the life history localities) and additional Ny. darlingi collected between 2006–2016 from Brazil (n = 7, additional states: Pará, São Paulo, Acre, Espirito Santo, Mato Grosso ) and Peru (n = 4, Lupuna and Cahuide, Loreto Department ). The process radtags program was used to drop low quality sequence reads and ustacks aligned reads into stacks with the following parameters: minimum depth of coverage for stack creation was set to 3, maximum distance allowed between stacks set to 4, and the maximum distance allowed to align secondary reads to primary reads set to 6. This generated a master catalog from Ny. darlingi sequences using consistent loci which allows for parallels to be drawn from different research projects.
nextRAD data analysis
For the present study, the sequences of the 93 Ny. darlingi from the seven collection localities (Table 1) were processed with process radtags and ustacks program as described above, compared against the above described catalog, and then SNPs were called with the STACKS -sstacks, -tsv2bam and -gstacks programsʼ settings set to default. The STACKS populations program was used to select a single SNP from each locus found in at least 40% of individuals in the dataset, a threshold slightly more lenient than the 50% used in previous population structure research ; this modification resulted in a greater number of loci for comparison.
STRUCTURE analysis was run using StrAuto v1.0 , allowing for parallel computation. To test the hypothesis of distinctive sub-populations within the major Brazilian biomes, a Bayesian clustering analysis was performed using the STRUCTURE admixture model assuming correlated allele frequencies for 10 replicates each of K = 1 through 7, with a ‘burn-in’ of 50,000 generations and a Markov chain Monte Carlo (MCMC) chain of 500,000 generations. CLUMPAK  was used to average runs and visualize STRUCTURE results. Principal components analysis (PCA) was conducted to test the hypothesis that the reduction of variables to principal components would lead to population separation based on SNP variation congruent with the populations from the Bayesian analysis. We performed PCA with the STRUCTURE file comparing different levels of population (biome, state, locality) in R (v. 3.6.0) using the ade4 package v.1.7.13  via the dudi.pca() function and visualized with factoextra package v.1.0.5  fviz_pca_ind() function. To partition the genetic variation into clusters, and confirm the optimal cluster number, we employed discriminant analysis of principal components (DAPC)  with the R package adegenet v.2.1.1 . A hierarchical analysis of molecular variance (AMOVA), with individuals grouped by locality within states, was calculated using the poppr.amova function in the R package poppr v.2.8.3 .
Life history trait analysis
All statistical analyses were conducted in R (v. 3.6.0) (Additional file 3: Dataset S2). A generalized linear model (GLM) was used to compare the effects of population (localities within state and families within localities) and temperature on life history traits. Genetic variation (populations or families), phenotypic plasticity (temperature levels) and genetic variation for plasticity (population/family-by-temperature interactions) were assessed for larval development, adult lifespan, and wing length with ANOVA (Type II). Comparison of localities within state was conducted on all states except for Rio de Janeiro because there was only one locality in that state (Fig. 1, Table 1). The Kaplan-Meier estimate of survival (time between larvae hatch and adult death) of individual families within each locality by temperature was visualized with the survival v.188.8.131.52  and survminer v.0.4.3  R packages.
Estimators of population differentiation: FST and PST
Trait-based data can be used to estimate the amount of genetic variance among populations (PST), which we expect to be comparable to calculated FST. In order to compare population genetic structure results from the sequenced P generation (field-collected females) and their laboratory-reared F1 progeny, FST and PST values were calculated, respectively. Pairwise FST values by locality (for the 7 localities tested in the present study) were calculated with the populations program from STACKS  using the sequenced 93 field Ny. darlingi. Life history data of the reared progeny were used to estimate PST, a phenotype-based analog for FST that measures the amount of genetic variation among populations relative to the total genetic variation , assuming that the proportion of phenotypic variance due to genetic effects is equivalent between and within populations. Pairwise PST values by locality for each life history trait were calculated with the Pstat R package v.1.2 . The ratio of PST to FST values is a useful proxy for estimating the strength of selection  on phenotypic traits. As FST is typically estimated from neutral loci, deviations of PST from FST can lead to inferences of selection: if PST > FST, directional selection can be inferred, conversely, if PST < FST, stabilizing selection is indicated.
Evidence of population genetic structure by major biome
There was an average of 3,891,842 (range: 359,767–6,636,895) sequences, or reads, per individual (n = 93) after quality filtering. The average number of reads per stack (or unique groups of matched reads) was 3,002,165 (range: 228,591–5,437,712) with an average number of 100,369 (range: 23,754–232,583) stacks per individual. The final SNP dataset included one biallelic SNP from each locus genotyped in at least 40% of the 93 individuals, for a total of 33,759 loci. The average coverage depth was 49X at each locus. Multiple values for K (K = 1–7) clusters were examined from STRUCTURE and STRUCTURE Harvester analyses (Additional file 4: Figure S1, Additional file 5: Figure S2). There was a dramatic dip in ΔK at K = 3 and the highest probability of K was at K = 3 (Additional file 4: Figure S1). The Bayesian information criterion (BIC) of the K-means clustering algorithm implemented in adegenet  in preparation for DAPC indicated K = 3 was the optimal number of clusters (Additional file 4: Figure S1b). Small amounts of admixture within the southern Amazonia biome (Rondônia state) samples were detected (Fig. 2a). Both the STRUCTURE (Fig. 2a, Additional file 5: Figure S2) and PCA (Fig. 2b) analyses identified three major clusters corresponding to the biome classifications.
Variation of life history traits between localities within states
Increasing temperature reduced larval development time, adult lifespan, and wing length within all states (Additional file 1: Table S2, Additional files 6–8: Figures S3-S5). There were significant genetic differences between localities within Amazonas (F(1, 968) = 52.0, P < 0.0001), Rondônia (F(1, 1049) = 15.3, P < 0.0001) and Tocantins (F(1, 332) = 6.7, P = 0.01) states for larval development time (Additional file 6: Figure S3), and only within Tocantins for adult lifespan (F(1, 332) = 4.57, P = 0.03) (Additional file 7: Figure S4) and wing length (F(1, 320) = 32.9, P < 0.0001) (Additional file 8: Figure S5).
The two localities within Amazonas State had significantly different larval development time at 20 °C (t(968) = 3.77, P < 0.0001), whereas the two localities in Rondônia had significantly different larval development time at 20 °C (t(1049) = 5.23, P < 0.0001) and 28 °C (t(1049) = − 3.41, P < 0.0001). Only the localities in Amazonas State had significantly different adult lifespan (t(968) = − 2.05, P = 0.04) and wing length (t(940) = 2.44, P < 0.0001) at 24 °C. In contrast, localities in Tocantins had significantly different wing lengths at 20 °C (t(320) = 2.44, P = 0.02) (Additional file 1: Table S2).
Within-population genetic variation for traits and their plastic responses
There was significant genetic variation among families within populations for larval development time (ARS: F(9, 369) = 6.71, P < 0.0001; APR: F(14, 530) = 5.48, P < 0.0001; RPV: F(14, 461) = 4.61, P < 0.0001; RMO: F(14, 504) = 1.77, P = 0.04), adult lifespan (ARS: F(9, 369) = 3.74, P < 0.0001; APR: F(14, 530) = 4.88, P < 0.0001; RPV: F(14, 461) = 3.94, P < 0.0001; RMO: F(14, 504) = 3.96, P < 0.0001) and wing length (ARS: F(9, 362) = 2.07, P = 0.03; APR: F(14, 509) = 2.66, P < 0.0001; RPV: F(14, 449) = 9.03, P < 0.0001; RMO: F(14, 490) = 6.72, P < 0.0001) for both localities of Amazonas and Rondônia states (Figs. 3, 4, 5), and for adult lifespan (SJU: F(10, 254) = 2.30, P = 0.01) and wing length (SJU: F(10, 247) = 3.75, P < 0.0001) in the southern high latitude population. Populations from Tocantins showed little genetic variation among families, with the exception of wing length for TLC (F(7, 214) = 6.82, P < 0.0001).
All traits from all populations showed significant responses to temperature - increasing temperature reduced larval development time, adult lifespan, and wing length in all localities (Figs. 3, 4, 5, Additional file 1: Table S2). The genotype-by-environment term was significant (or nearly so) for all three traits in multiple populations, indicating significant genetic variation among families for the response to temperature (different slopes among families; Figs. 3, 4, 5). Median survival times (larval hatch to adult death) of families within each locality were highest at 20 °C and decreased with increasing temperature. Family survival was significantly different within each locality at each temperature (Additional file 9: Figure S6).
Estimators of population differentiation: AMOVA, F ST and P ST
The hierarchical AMOVA (Table 2) revealed highly significant levels of variation at each level (P < 0.001), and the genetic variation was mainly accounted for within individuals (71.8%), followed by among individuals (19.1%) and between states (9.1%). Pairwise FST between localities of the field collected females (P generation) ranged from 0.045 to 0.183, with the lowest FST values between paired sites within the same state, as expected (Table 3). The calculated pairwise PST values across all lab-reared progeny (F1 generation) by locality were generally greater than FST for larval development (range: 0.660–0.995) (Table 4), adult lifespan (range: 0.0004–0.972) and wing length (range: 0.211–0.994) (Table 5), with few exceptions. The lowest pairwise PST values for these three traits were between localities within the same state (Tables 4, 5). Comparison of PST values for progeny life history traits with parental FST values reveals that PST values are nearly all substantially greater, an indication that directional selection is responsible for some of the genetic differentiation of life history traits among regions.
Using a representative sample of fine-scale SNP data, and an extensive common-garden experiment, we found strong evidence of among-population genetic differentiation as well as within-population variation in plasticity in major life history traits. In particular, our SNP analysis found population differentiation according to major biome designation rather than at a local scale, similar to previous findings . We also found evidence of variation in plasticity in life history traits at local scales. The genetic structure of Ny. darlingi populations as well as the local-level variation for plasticity of this vector has major implications for the future of malaria elimination in South America. A recent report from the World Health Organization highlighted the importance of novel and local approaches as vital for malaria elimination  and our findings indicate that the variation at the local level could enable some populations to potentially tolerate changing temperature (i.e. Amazonia biome) and for potential increased local transmission in southern populations (i.e. Rio de Janeiro).
In partial support of previous findings of differentiation according to biome and physical barriers , populations analyzed in the present study clustered by major biome. A new finding was the detection of low admixture between the two states at the same latitude (Rondônia, Tocantins) (Fig. 2). Our study revealed geographical divisions, with pairwise FST between localities (34–120 km apart) lowest within (range: 0.046–0.070) compared to between biomes (range: 0.081–0.183) (Table 3) indicating weaker genetic differentiation at smaller geographical scales. These findings suggest that biome boundaries may represent strong barriers to gene flow in Ny. darlingi. There is limited evidence from SNP analysis of significant microgeographical genetic differentiation of Ny. darlingi from western Amazonian Brazil related to varying levels of deforestation among municipalities (that ranged from 60–1600 km apart); one study detected significant differentiation when comparing an older, highly deforested agricultural settlement 60 km from a newly settled one that retained high forest cover levels , and another detected low and non-significant variation when comparing multiple deforestation levels among several Brazilian Amazon settlements .
A study of the closely related species An. gambiae (s.s.) and An. coluzzii, collected from 15 sites across Africa and tested with over 50 million SNPs, revealed clustering by geographical region rather than species, and as predicted, lower FST values within rather than between biomes . Population structure of Anopheles species that have been analyzed is strongly affected by geographical division; such demarcations may break down in the future as the integrity of biomes is eroded by deforestation and climate change. Specifically, in South America, under a high CO2 emission model, there could be substantial reduction (3%) in tropical forest area in South America in the next 10 years and up to 18% by 2100 . As Ny. darlingi is primarily associated with forested areas, its range  and population structure will likely be altered. Variation at the individual level was high in this study (72%) and suggests potential for adaptation.
Our study extends previous findings of regional variation in life history traits of Ny. darlingi  to the detection of significant genetic variation within localities. Significant genetic variation between families was found consistently within populations within the Amazonia biome (larval development time, adult lifespan and adult body size) and families from the Mata Atlântica exhibited significant genetic variation for adult lifespan and body size. These populations (Amazonas and Mata Atlântica) may have greater adaptive potential to increase their resistance to changing temperature given their responses in the laboratory experiment. PST values  of localities within each state were the lowest for all three life history traits (larval development time, adult lifespan and wing length), as expected. Because PST values are nearly uniformly substantially greater than the parental FST values, we infer that there is directional selection driving genetic differentiation of life history traits among regions. Combined with our evidence that there is standing genetic variation for performance at different temperatures within populations, future selection could favor phenotypes that tolerate increased temperatures.
The environment substantially influences mosquito vector traits. In our study, increased temperature reduced larval development time in all populations with the magnitude of reduction in adult lifespan and body size dependent on population. The dramatic differences in larval development time however, were not a clear predictor of adult longevity due to differences between populations. Our data show that the relationship between larval conditions and adult traits is not linear but rather complex. Temperature of larval and adult environments had significant effects on An. gambiae (s.s.) development: increased larval rearing temperature (23–31 °C) resulted in smaller larvae and adults whereas increased adult temperatures reduced the proportion of egg hatch . The effects of temperature can reduce mosquito population size over time, with smaller individuals laying fewer eggs. On the other hand, smaller increases in temperature can increase mosquito population size, with a field study of An. gambiae (s.s.) in Kenya revealing greater fecundity and vectorial capacity of mosquitoes placed in homes that were 0.7–1.2 °C warmer compared with control homes .
Vector control interventions need to consider variation in life history traits , behaviors , and habitats . At a high nutrition diet, low temperature treatment Anopheles were found to be larger and more likely to survive exposure to a LC50 dose of permethrin . Data from our study suggest that larger doses of permethrin would be required in southern compared to northern populations of Ny. darlingi. Interventions such as long-lasting insecticidal nets (LLINs) are highly effective and target adult mosquitoes that are mainly endophagic and endophilic. The biting behavior of Ny. darlingi is variable [41, 61], compromising the efficacy of IRS or LLINs. Field studies of Ny. darlingi reveal endophagy and exophagy at different times throughout the night , and there is no evidence for a genetic basis of these behaviors . Ivermectin treatment of cattle was shown to reduce An. arabiensis fecundity by nearly 60% after deployment of LLINs compared to LLINs alone, supporting the use of combination interventions to help achieve population elimination .
The plasticity of Ny. darlingi, including biting behavior [5, 41], host [12, 63] and breeding site [19,20,21] preferences, coupled with the potential of families within certain populations to withstand changing environments, help explain its status as the major malaria vector in South America. The present study contributes to the growing body of evidence of high levels of plasticity in Ny. darlingi, and significantly, presents evidence for genetic variation in plasticity within populations.
A potential study limitation was that we were unable to collect Ny. darlingi from a second locality in Rio de Janeiro although we had previous evidence of its presence . This limited the comparison of life history trait responses from paired localities within Rio de Janeiro State and reduced our ability to adequately test population structure within this biome. While the FST values were calculated from field collected mosquitoes and PST values from laboratory-reared progeny, it is unlikely that there would accrue significant genetic variation between parent and offspring in one generation. To date, Ny. darlingi has not been tested for polyandry, although An. gambiae exhibits low polyandry (12%)  whereas nearly 25% of An. arabiensis females had been multiply inseminated . We treated individuals within families as full-siblings assuming that polyandry did not contribute significantly to the observed variation between families.
The study design of our laboratory experiment to observe life history traits at a constant temperature throughout the mosquito life cycle was somewhat limited by space and resources. The temperatures in this experiment were chosen to avoid extremes that can lead to excessive mortality that would limit comparisons , and reflect averages at specific latitudes; they may not reflect specific microclimates for each locality. The temperature range (8 °C) we used may not reflect actual temperatures projected for Brazil under climate change [16, 66, 67]. Research with established laboratory colonies has also shown that fluctuating temperatures may more accurately reflect the natural environment, and affect life history traits differently compared to constant temperatures .
Our treatment of adult mosquitoes (providing only water, no food) deviates from the natural adult environment which involves sugar feeding and potential blood meals for females as well as sugar-feeding for males. On the other hand, the average adult longevity in our study was 3.09 days compared to field data of daily survival rates that detected between 3.73 and 23.9 days for adult females from two Peruvian sites . Our research did not investigate variation in biting behavior, fecundity and susceptibility to Plasmodium that can be affected by temperature and Ny. darlingi population specificity. The establishment and maintenance of a laboratory colony in Peru  and Brazil  as well as successful Plasmodium infection of colony mosquitoes  would facilitate investigations of this variation over generations and between populations.
This study identified the population structure and degree of genetic variation and phenotypic plasticity of Ny. darlingi in Brazil. The genomic signatures indicate that genetic divergence occurs at the level of biomes, with phenotypic traits varying more than molecular markers indicating a role for natural selection via climate or vegetation structure in driving differentiation. A key result is our finding that there is genetic variation for both life history traits and their plastic responses within populations to temperature, indicating future adaptive capacity to changes in temperature. Future research that further quantifies the effects of environment and population on life history traits relevant to transmission will be vital for predicting variation in transmission potential and informing modeling efforts.
Availability of data and materials
Raw Illumina sequences were deposited in the NCBI Sequence Read Archive (SRA; BioProject ID PRJNA576174). All other relevant data are within the manuscript and its additional files.
analysis of molecular variance
analysis of variance
Manaus-Brasilierinho, Amazonas State, Brazil
Ramal Novo Horizonte, Amazonas State, Brazil
Bayesian information criterion
discriminant analysis of principal components
- F1 :
offspring generation 1
- F ST :
generalized linear model
long-lasting insecticidal nets
Markov chain Monte Carlo
nextera-tagmented, reductively-amplified DNA
principle components analysis
- P ST :
degree of phenotypic differentiation
- R0 :
basic reproductive rate
Machadinho dʼOeste, Rondonia State, Brazil
Porto Velho, Rondonia State, Brazil
Lake Juturnaiba, Rio de Janeiro State, Brazil
single nucleotide polymorphism
Lagoa da Confusão, Tocantins State, Brazil
Porto Nacional, Tocantins State, Brazil
Recht J, Siqueira AM, Monteiro WM, Herrera SM, Herrera S, Lacerda MVG. Malaria in Brazil, Colombia, Peru and Venezuela: current challenges in malaria control and elimination. Malar J. 2017;16:273.
WHO. World malaria report 2018. Geneva: World Health Organization; 2018. https://www.who.int/malaria/publications/world-malaria-report-2018/en/
Foster PG, Oliveira de TMP, Bergo ES, Conn JE, Sant’Ana DC, Nagaki SS, et al. Phylogeny of Anophelinae using mitochondrial protein coding genes. R Soc Open Sci. 2017;4:170758.
Campos M, Alonso DP, Conn JE, Vinetz JM, Emerson KJ, Ribolla PEM. Genetic diversity of Nyssorhynchus (Anopheles) darlingi related to biting behavior in western Amazon. Parasit Vectors. 2019;12:242.
Prado CC, Alvarado-Cabrera LA, Camargo-Ayala PA, Garzón-Ospina D, Camargo M, Soto-De León SC, et al. Behavior and abundance of Anopheles darlingi in communities living in the Colombian Amazon riverside. PLoS ONE. 2019;14:e0213335.
Motoki MT, Suesdek L, Bergo ES, Sallum MAM. Wing geometry of Anopheles darlingi Root (Diptera: Culicidae) in five major Brazilian ecoregions. Infect Genet Evol. 2012;12:6.
Chu VM, Sallum MAM, Moore TE, Lainhart W, Schlichting CD, Conn JE. Regional variation in life history traits and plastic responses to temperature of the major malaria vector Nyssorhynchus darlingi in Brazil. Sci Rep. 2019;9:5356.
Voorham J. Intra-population plasticity of Anopheles darlingi’s (Diptera, Culicidae) biting activity patterns in the state of Amapá. Brazil. Rev Saude Publica. 2002;36:1.
Da Rocha JAM, De Oliveira SB, Póvoa MM, Moreira LA, Krettli AU. Malaria vectors in areas of Plasmodium falciparum epidemic transmission in the Amazon Region. Brazil. Am J Trop Med Hyg. 2008;78:6.
Lainhart W, Bickersmith SA, Nadler KJ, Moreno M, Saavedra MP, Chu VM, et al. Evidence for temporal population replacement and the signature of ecological adaptation in a major Neotropical malaria vector in Amazonian Peru. Malar J. 2015;14:375.
Galardo AKR, Arruda M, D’Almeida Couto AAR, Wirtz RA, Lounibos PL, Zimmerman RH. Malaria vector incrimination in three rural riverine villages in the Brazilian amazon. Am J Trop Med Hyg. 2007;76:3.
Moreno M, Saavedra MP, Bickersmith SA, Prussing C, Michalski A, Tong Rios C, et al. Intensive trapping of blood-fed Anopheles darlingi in Amazonian Peru reveals unexpectedly high proportions of avian blood-meals. PLoS Negl Trop Dis. 2017;11:1–19.
Phasomkusolsil S, Pantuwattana K, Tawong J, Khongtak W, Kertmanee Y, Monkanna N, et al. The relationship between wing length, blood meal volume, and fecundity for seven colonies of Anopheles species housed at the Armed Forces Research Institute of Medical Sciences, Bangkok, Thailand. Acta Trop. 2015;152:220–7.
Shapiro LLM, Murdock CC, Jacobs GR, Thomas RJ, Thomas MB. Larval food quantity affects the capacity of adult mosquitoes to transmit human malaria. Proc R Soc B Biol Sci. 2016;283:20160298.
Barreaux AMG, Stone CM, Barreaux P, Koella JC. The relationship between size and longevity of the malaria vector Anopheles gambiae (s.s.) depends on the larval environment. Parasit Vectors. 2018;11:485.
IPCC 2013. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge: Cambridge University Press; 2013.
Murdock CC, Sternberg ED, Thomas MB. Malaria transmission potential could be reduced with current and future climate change. Sci Rep. 2016;6:27771.
Sinka ME, Rubio-Palis Y, Manguin S, Patil AP, Temperley WH, Gething PW, et al. The dominant Anopheles vectors of human malaria in the Americas: occurrence data, distribution maps and bionomic précis. Parasit Vectors. 2010;3:72.
Barros FSM, Arruda ME, Gurgel HC, Honório NA. Spatial clustering and longitudinal variation of Anopheles darlingi (Diptera: Culicidae) larvae in a river of the Amazon: the importance of the forest fringe and of obstructions to flow in frontier malaria. Bull Entomol Res. 2011;101:643–58.
Conde M, Pareja PX, Orjuela LI, Ahumada ML, Durán S, Jara JA. Larval habitat characteristics of the main malaria vectors in the most endemic regions of Colombia: potential implications for larval control. Malar J. 2015;14:476.
Arcos AN, Ferreira FA, Cunha HB, Tadei WP. Characterization of artificial larval habitats of Anopheles darlingi (Diptera: Culicidae) in the Brazilian Central Amazon. Rev Bras Entomol. 2018;62:267–74.
Vittor A, Gilman RH, Tielsch J, Glass G, Shields T, Pinedo-Cancino V, et al. Linking deforestation to malaria in the Amazon: characterization of the breeding habitat of the principle malaria vector, Anopheles darlingi. Am J Trop Med Hyg. 2009;81:5–12.
Chaves LSM, Conn JE, López RVM, Sallum MAM. Abundance of impacted forest patches less than 5 km2 is a key driver of the incidence of malaria in Amazonian Brazil. Sci Rep. 2018;8:7077.
Sallum MAM, Conn JE, Bergo ES, Laporta GZ, Chaves LSM, Bickersmith SA, et al. Vector competence, vectorial capacity of Nyssorhynchus darlingi and the basic reproduction number of Plasmodium vivax in agricultural settlements in the Amazonian Region of Brazil. Malar J. 2019;18:117.
Moller-Jacobs LL, Murdock CC, Thomas MB. Capacity of mosquitoes to transmit malaria depends on larval environment. Parasit Vectors. 2014;7:593.
Beck-Johnson LM, Nelson WA, Paaijmans KP, Read AF, Thomas MB, Bjørnstad ON. The effect of temperature on Anopheles mosquito population dynamics and the potential for malaria transmission. PLoS One. 2013;8:11.
Coetzee M, Hunt RH, Wilkerson R, Torre AD, Coulibaly MB, Besansky NJ. Anopheles coluzzii and Anopheles amharicus, new members of the Anopheles gambiae complex. Zootaxa. 2013;3619:246–74.
Lehmann T, Dalton R, Kim EH, Dahl E, Diabate A, Dabire R, et al. Genetic contribution to variation in larval development time, adult size, and longevity of starved adults of Anopheles gambiae. Infect Genet Evol. 2006;6:410–6.
Deutsch CA, Tewksbury JJ, Huey RB, Sheldon KS, Ghalambor CK, Haak DC, et al. Impacts of climate warming on terrestrial ectotherms across latitude. Proc Natl Acad Sci USA. 2008;105:6668–72.
Christiansen-Jucht CD, Parham PE, Saddler A, Koella JC, Basáñez M-G. Larval and adult environmental temperatures influence the adult reproductive traits of Anopheles gambiae s.s. Parasit Vectors. 2015;8:456.
Christiansen-Jucht C, Parham PE, Saddler A, Koella JC, Basáñez MG. Temperature during larval development and adult maintenance influences the survival of Anopheles gambiae s.s. Parasit Vectors. 2014;7:489.
Mordecai EA, Paaijmans KP, Johnson LR, Balzer C, Ben-Horin T, de Moor E, et al. Optimal temperature for malaria transmission is dramatically lower than previously predicted. Ecol Lett. 2013;16:22–30.
Fouet C, Kamdem C, Gamez S, White BJ. Extensive genetic diversity among populations of the malaria mosquito Anopheles moucheti revealed by population genomics. Infect Genet Evol. 2017;48:27–33.
Wiltshire RM, Bergey CM, Kayondo JK, Birungi J, Mukwaya LG, Emrich SJ, et al. Reduced-representation sequencing identifies small effective population sizes of Anopheles gambiae in the north-western Lake Victoria basin, Uganda. Malar J. 2018;17:285.
Cheng C, White BJ, Kamdem C, Mockaitis K, Costantini C, Hahn MW, et al. Ecological genomics of Anopheles gambiae along a latitudinal cline: a population-resequencing approach. Genetics. 2012;190:1417–32.
Maliti DV, Marsden CD, Main BJ, Govella NJ, Yamasaki Y, Collier TC, et al. Investigating associations between biting time in the malaria vector Anopheles arabiensis Patton and single nucleotide polymorphisms in circadian clock genes: support for sub-structure among An. arabiensis in the Kilombero valley of Tanzan. Parasit Vectors. 2016;9:109.
Emerson KJ, Conn JE, Bergo ES, Randel MA, Sallum MAM. Brazilian Anopheles darlingi root (Diptera: Culicidae) clusters by major biogeographical region. PLoS One. 2015;10:7.
Campos M, Conn JE, Alonso DP, Vinetz JM, Emerson KJ, Ribolla PEM. Microgeographical structure in the major Neotropical malaria vector Anopheles darlingi using microsatellites and SNP markers. Parasit Vectors. 2017;10:76.
Linthicum K. A revision of the Argyritarsis section of the Nyssorhynchus of Anopheles (Diptera: Culicidae). Mosq Syst. 1988;20:2.
Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA. Stacks: an analysis tool set for population genomics. Mol Ecol. 2013;22:11.
Prussing C, Moreno M, Saavedra MP, Bickersmith SA, Gamboa D, Alava F, et al. Decreasing proportion of Anopheles darlingi biting outdoors between long-lasting insecticidal net distributions in peri-Iquitos, Amazonian Peru. Malar J. 2018;17:86.
Prussing C, Emerson KJ, Bickersmith SA, Sallum MAM, Conn JE. Minimal genetic differentiation of the malaria vector Nyssorhynchus darlingi associated with forest cover level in Amazonian Brazil. PLoS ONE. 2019;14:e0225005.
Chhatre VE, Emerson KJ. StrAuto: automation and parallelization of STRUCTURE analysis. BMC Bioinform. 2017;18:192.
Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA. CLUMPAK: a program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour. 2015;15:5.
Dray S, Dufour A-B. The ade4 package: Implementing the duality diagram for ecologists. J Stat Softw. 2015;22:4.
Kassambara, Alboukadel Mundt F. Factoextra: Extract and visualize the results of multivariate data analyses. 2017. https://cran.r-project.org/package=factoextra.
Jombart T. Adegenet: A R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–5.
Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010;11:94.
Kamvar ZN, Tabima JF, Grünwald NJ. Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ. 2014;2:e281.
Therneau TM. A package for survival analysis in S. 2015. https://cran.r-project.org/package=survival.
Kosinski AK and M. survminer: Drawing Survival Curves using “ggplot2”. 2017. https://cran.r-project.org/package=survminer.
Leinonen T, McCairns RJS, O’Hara RB, Merilä J. QST -FST comparisons: evolutionary and ecological insights from genomic heterogeneity. Nat Rev Genet. 2013;14:179–90.
Stephane BDS. Pstat: Assessing Pst Statistics. R package version 1.2. 2017. https://cran.r-project.org/package=Pstat.
Raeymaekers JAM, Van Houdt JKJ, Larmuseau MHD, Geldof S, Volckaert FAM. Divergent selection as revealed by PST and QTL-based FST in three-spined stickleback (Gasterosteus aculeatus) populations along a coastal-inland gradient. Mol Ecol. 2007;16:4.
The WHO Strategic Advisory Group. Malaria eradication: benefits, future scenarios and feasibility (Executive summary). 2019. https://www.who.int/publications-detail/strategic-advisory-group-malaria-eradication-executive-summary.
The Anopheles gambiae 1000 Genomes Consortium. Genetic diversity of the African malaria vector Anopheles gambiae. Nature. 2017;552:7683.
Salazar LF, Nobre CA, Oyama MD. Climate change consequences on the biome distribution in tropical South America. Geophys Res Lett. 2007;34:9.
Laporta GZ, Linton Y-M, Wilkerson RC, Bergo ES, Nagaki SS, Sant’Ana DC, et al. Malaria vectors in South America: current and future scenarios. Parasit Vectors. 2015;8:426.
Afrane YA, Little TJ, Lawson BW, Githeko AK, Yan G. Deforestation and vectorial capacity of Anopheles gambiae Giles mosquitoes in malaria transmission, Kenya. Emerg Infect Dis. 2008;14:1533–8.
Owusu HF, Chitnis N, Müller P. Insecticide susceptibility of Anopheles mosquitoes changes in response to variations in the larval environment. Sci Rep. 2017;7:3667.
Hiwat H, Bretas G. Ecology of Anopheles darlingi Root with respect to vector importance: a review. Parasit Vectors. 2011;4:177.
Ng’habi K, Viana M, Matthiopoulos J, Lyimo I, Killeen G, Ferguson H. Mesocosm experiments reveal the impact of control measures on malaria vector life history and population dynamics. Rev. 2018;8:13949.
Zimmerman RH, Galardo AKR, Lounibos PL, Arruda M, Wirtz RA. Bloodmeal hosts of Anopheles species (Diptera: Culicidae) in a malaria-endemic area of the Brazilian Amazon. J Med Entomol. 2006;43:947–56.
Gomulski L. Polyandry in nulliparous Anopheles gambiae mosquitoes (Diptera: Culicidae). Bull Entomol Res. 1990;80:393–6.
Helinski ME, Hood RC, Knols BG. A stable isotope dual-labelling approach to detect multiple insemination in un-irradiated and irradiated Anopheles arabiensis mosquitoes. Parasit Vectors. 2008;1:9.
INMET. Climate data, National Institute of Meteorology; 2015. http://www.inmet.gov.br/portal/ Accessed 6 Jun 2015.
SEDAM. Climate data from the Government of the State of Rondônia State Secretariat for Environmental Development. http://www.sedam.ro.gov.br/ Accessed 8 Apr 2017.
Davies C, Coetzee M, Lyons CL. Effect of stable and fluctuating temperatures on the life history traits of Anopheles arabiensis and An. quadriannulatus under conditions of inter- and intra-specific competition. Parasit Vectors. 2016;9:342.
Villarreal-Treviño C, Vásquez GM, López-Sifuentes VM, Escobedo-Vargas K, Huayanay-Repetto A, Linton Y-M, et al. Establishment of a free-mating, long-standing and highly productive laboratory colony of Anopheles darlingi from the Peruvian Amazon. Malar J. 2015;14:227.
Araujo S, Andrade AO, Akira N, Pereira B, Costa S. Brazil’s first free-mating laboratory colony of Nyssorhynchus darlingi. Rev Soc Bras Med Trop. 2019;52:e20190159.
Moreno M, Tong C, Guzmán M, Chuquiyauri R, Llanos-Cuentas A, Rodriguez H, et al. Infection of laboratory-colonized Anopheles darlingi mosquitoes by Plasmodium vivax. Am J Trop Med Hyg. 2014;90:4.
Google Earth Pro. Brazil, eye alt 247.82mi. 2018.
We express our sincere gratitude to Denise C. Sant’Ana, Eduardo S. Bergo, and Leonardo Chaves for assistance in field collections and to Caio Morreira for laboratory rearing guidance. We thank D. Gamboa (Universidad Peruana Cayetano Heredia, Lima, Peru) and JM Vinetz (Yale University School of Medicine, New Haven, CT, USA) for specimens from Loreto state, Peru used in the creation of the catalog of Nyssorhynchus darlingi. We thank Sara Bickersmith for help with DNA extraction and sample storage. We greatly appreciate the co-operation and support of the local communities and health workers of Juturnaiba, Rio de Janeiro; Lagoa da Confusão, and Porto Nacional, Tocantins; Machadinho d’Oeste, and Porto Velho, Rondônia; Ramal Novo Horizonte, and Brasilierinho, Manaus, Amazonas. We are grateful to the Coleção Entomológica de Referência and Faculdade de Saúde Pública, Universidade de São Paulo for extensive use of laboratory facilities.
This work has been funded by the National Institutes of Health, USA (grant R01 AI110112 to JEC). MAMS is funded by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP, grant no. 2014/26229-7) and Conselho Nacional de Pesquisa (CNPq no. 301877/2016-5). The Boren Awards fellowship (2015) and the Biodefense and Emerging Infectious Disease training fellowship grant (T32AI05532901) provided partial funding support to VMC.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Number of larvae and families from each locality at each temperature. Table S2. Generalized linear model results for larval development, adult lifespan, and wing length. Table S3. Supplemental collection site data.
STRUCTURE file used for STRUCTURE, PCA, DAPC and FST analyses, with genotypes for 93 individuals at 33,759 loci.
Life history data used for ANOVA, PST, GLM and Kaplan–Meier analyses.
Estimation of the number of clusters in SNP dataset using STRUCTURE and Discriminant Analysis of Principal Components (DAPC). a STRUCTURE ΔK using Evanno method (left) and probability by K (right)  for K = 1–7. b Bayesian information criterion (BIC) for 1 to 20 clusters using K means clustering preparation for DAPC.
STRUCTURE plots for K = 1–7.
Average larval development time of localities within each state. Standard error bars and ANOVA results in each panel: G, genetic variation (locality); E, phenotypic variation (temperature); GEI, genotype-by-environment interaction (locality × temperature); *P < 0.05, **P < 0.01, ***P < 0.001.
Average adult lifespan of localities within each state. Standard error bars and ANOVA results in each panel: G, genetic variation (locality); E, phenotypic variation (temperature); GEI, genotype-by-environment interaction (locality × temperature); *P < 0.05, **P < 0.01, ***P < 0.001.
Average wing length of localities within each state. Standard error bars and ANOVA results in each panel: G, genetic variation (locality); E, phenotypic variation (temperature); GEI, genotype-by-environment interaction (locality × temperature); *P < 0.05, **P < 0.01, ***P < 0.001.
Kaplan-Meier survival curve of families within locality. Each family is uniquely color coded and consistent with all other figures’ color scheme for family.
About this article
Cite this article
Chu, V.M., Sallum, M.A.M., Moore, T.E. et al. Evidence for family-level variation of phenotypic traits in response to temperature of Brazilian Nyssorhynchus darlingi. Parasites Vectors 13, 55 (2020). https://doi.org/10.1186/s13071-020-3924-7