Host selection pattern and flavivirus screening of mosquitoes in a disturbed Colombian rainforest

Studies on the feeding behavior of hematophagous insects, particularly those of medical importance, are relevant for tracking possible pathogen transmission routes and identifying biases in the choice of vertebrates. We evaluated host selection of blood-feeding mosquitoes in a disturbed forest in the Magdalena Medio valley in Colombia from March 2017 to April 2018, after the introduction of Zika virus to the Americas from the 2015–2016 outbreak. We estimated vertebrate diversity and collected blood-engorged female mosquitoes. Genomic DNA/RNA was extracted from the mosquito’s abdomen for vertebrate host identification and pathogen detection. We performed conventional PCR and sequencing, using universal primers targeting vertebrate regions of the eukaryotic mitochondrial genome to determine bloodmeal host. Additionally, we tested for the presence of flaviviruses in all mosquito samples with RT-PCR. Based on the identity and quantity of detected bloodmeals, we performed mosquito-vertebrate interaction network analysis and estimated topology metrics. In total, we collected 292 engorged female mosquitoes representing 20 different species. Bloodmeal analyses identified 26 vertebrate species, the majority of which were mammals (N = 16; 61.5%). No flaviviruses of medical importance were detected from the samples. Although feeding patterns varied, network analyses showed a high degree of specialization by mosquitoes and revealed ecological and phylogenetic relationships among the host community. We conclude that host selection or preference by mosquitoes is species specific.

Vertebrate diversity. During fieldwork, we recorded 90 bird species belonging to 23 families using mist net collections and visual observations in three habitat types: pasture, forest understory, and river edge. Avian families with the greatest observed local diversity were Tyrannidae (13 species) followed by Ardeidae (egrets), Psittacidae (parrots), Icteridae (New World blackbirds, cowbirds, orioles), Picidae (woodpeckers) and Thraupidae (tanagers) with five species each. The bird species most frequently detected were Smooth-billed Ani (Crotophaga ani), Great Egret (Ardea alba), Ruddy Ground-Dove (Columbina talpacoti), Bicolored Wren (Camphylorynchus griseus) and Bare-faced Ibis (Phimosus infuscatus). Mammal populations were assessed based on captures of 164 individuals belonging to 21 species of nine families. Chiroptera was the order with highest richness with 15 species, followed by Primates (four species), Didelphimorphia (one species) and Rodentia (one species). Additionally, nine species of reptiles and eight of amphibians were recorded (Supplementary Table S1).  Network analysis. The mosquito-vertebrate network comprised 41 different species interactions (Fig. 2a).
Modularity was high (Q = 0.59, z = 13.4, p < 0.05), indicating significant specialization of mosquito choice of vertebrates utilized for bloodmeals. We obtained eight modules after 50 iterations (Fig. 3). Five of the modules were associated with a single mosquito species, meaning they showed a very different behavior compared to other species. The remaining modules grouped two, three or five species of mosquitoes. Regarding vertebrates, interestingly, we found one host module composed of domestic and synanthropic species, three exclusively of mammalian species, one of only reptile species and two only of amphibians. All species of marsupials (D. marsupialis, P. opposum, and M. caucae) were in the same module with the primate Alouatta seniculus and the rodent Zygodontomys brevicauda. Overall, vertebrates with similar traits in terms of habitat use, were grouped together. www.nature.com/scientificreports/ The specialization index at the community level (H2′) had a value of 0.55, indicating aggregated patterns of mosquito feeding behavior and preferences for certain vertebrates. In addition, the analysis of nestedness confirmed that mosquito-feeding behavior is not significantly nested (NODF = 6.84, p = 0.62). At the species level, the degree of mosquito specialization varied among species when comparing the d and COV indexes. The species Culex sp. 1 and Uranotaenia sp. 1 had the highest values for both indexes (d = 1, COV = 1) suggesting those were the species with the most exclusive interaction types and narrow host range. Anopheles neomaculipalpus and An. triannulatus had a low value of index d = 0.21, but high value for the COV index = 1 ( Table 2) as human blood was the only bloodmeal detected on these mosquitoes, but also in other species. Aedes fulvus, Ae. serratus Cx. adamesi, Cx. nigripalpus, and Cx. pedroi had the highest values for "species strength", denoting greater diversity in the detected interactions. Lastly, the value of betweenness centrality (BC) for host species varied from 0 to 0.39 (Fig. 2a).

Discussion
Using morphology and molecular analyses, we could identify 11 mosquito species previously involved in arbovirus transmission, out of 20 species of mosquitoes collected [15][16][17] , with the greatest number of species belonging to the Culex genus (six species). Species in the genus Culex, especially those belonging to the Melanoconion subgenus, are recognized as a highly diverse group with great taxonomic difficulties 18 , due to morphological www.nature.com/scientificreports/ similarity and the lack of recent taxonomic revisions of the Melanoconion subgenus 19,20 . The use of mitochondrial COI gene has become a great help in delimiting species within groups displaying ambiguous morphological traits 21 ; we contributed the COI sequences for 14 mosquito species from an unsampled ecoregion. The efficiency of species identification was supported by the barcoding gap analysis, which implies that if a gap exists, a cut-off value should be defined allowing a clear distinction between the interspecific and intraspecific distances. The K2P intraspecific and interspecific values obtained here are similar to values found in other studies that include related species 22 . We identified 26 species of vertebrates among 154 mosquito bloodmeals, and found that only 10 of these species were observed during fieldwork. Four species P. opposum, M. caucae, D. novemcinctus, and Z. brevicauda, were only detected through bloodmeal analyses and were expected at the study site based on available distribution records 23 . We only placed vertebrate traps at the ground level, so we missed arboreal mammals during sampling; future studies should include trapping at different levels in the forest to better characterize the vertebrate community 24 .
Human DNA was detected in 6.9% of the identified bloodmeal sequences in eight mosquito species, some previously known to feed on humans: Ae. serratus, Ps. albipes, Cx. nigripalpus, An. triannulatus and Cq. nigricans [25][26][27][28][29] . Studies in villages located in the Brazilian amazon forest found that Ae. fulvus fed on humans and peridomestic animals such as dogs, pigs, and cows 25 , while Silva et al. 26 found a greater tendency to ornithophily in this species. These data contrast our findings, considering that despite the abundance of samples obtained, we did not find evidence of Ae. fulvus feeding on humans or birds, and the selection of hosts for this species was dominated almost entirely by sylvatic mammalian species. Culex nigripalpus has been recognized as a species with broad host selection in laboratory-controlled experiments 28 , and in our study, it also showed a wide-ranging feeding behavior, mainly related to domestic animals. Similar to our findings, Mitchell et al. 30 identified Mansonia titillans frequently feeding on domestic species including Phasianidae birds and Bovidae. Preference of mosquitoes for vertebrates in nature can be modulated by many variables, including but not limited to relative abundance of vertebrate hosts, reproductive cycles, and host defensive behavior 31 . These ecological factors are difficult to monitor, limiting the understanding of innate behaviors of the vector [32][33][34][35] . One very specific behavior is that of the genus Uranotaenia which has been associated exclusively with ectothermic animals. We found Uranotaenia sp. 1 feeding on the amphibian Dendropsophus subocularis. The species U. lowii was one of the first mosquitoes in the Neotropics to show attraction towards the sound of an amphibian host 36 . In our study some Culex species were also attracted to ectothermic animals (frogs) as food source.
For most of the mosquito species in which human blood was identified, except An. triannulatus and An. neomaculipalpus, more than one host was found. Factors favoring humans as food sources and their epidemiological consequences merit further investigation. Frequent contact with a host leads to a strong vector-vertebrate interaction that improves the transmission of pathogens 37 . Mosquitoes with low specialization have special public health relevance, since they can disperse through a diverse range of ecosystems, becoming more likely to get infected with new pathogens, and raising the risk of spillover 10 . Conversely, a greater diversity of interacting vertebrate species favors a dilution effect, and this allows prevalence rates to remain low 38 .
Although the feeding patterns found in our study are varied, according to our analysis (significant modular topology Q) we conclude that they don't occur randomly, and there is a preference of mosquitoes for certain vertebrates. Among ecological networks, values of modularity and specialization tend to be higher in antagonisticrelationship networks if compared to other types of interactions such as mutualism 39,40 . The species involved in such antagonistic interactions tend to form enclosed modules within the interaction matrix 41 . This marked Table 2. Index values of specialization by mosquito species recognized in San Juan, Carare, Colombia. The Specialization index (d) and the Coefficient of interaction variation (COV) are estimated with scaled values between 0 to 1. www.nature.com/scientificreports/ compartmentalization within antagonistic networks (species grouped in discrete modules) occurs as a consequence of a series of pressures (evolutive, physiological, etc.), establishing restrictions over their interactions 40 . In a parasite-host modular network, the hosts of the same module are used by the same group of parasites, in our case by the same group of mosquitoes 42 . In our study site, for example, Aedes and Anopheles feed on a variety of sylvatic mammalian species, while only Uranotaenia and Culex fed on amphibians.
In modular networks, sub-communities defined by such modules are loosely connected to one another. Within this network topology, some nodes act as bridges connecting the different communities, facilitating perturbations to spread throughout the entire network 43 . The value of Betweeness Centrality calculated for each vertebrate species can be used to identify individuals that 'bridge' distant groups and may thus play an important role in spreading parasites among network sub-communities 44 . Among the detected vertebrates, chickens, common opossums and humans showed the highest values of betweenness centrality, denoting their potential role as bridge species between sub-communities. For instance the common opossum is a synanthropic species with home range on the urban-sylvatic interface, and host for a wide range of vector-borne pathogens 45,46 . The role of humans must also be considered for pathogen transmission from urban to sylvatic cycles, as has been shown in South America for imported pathogens such as YFV and dengue virus 1 .
Regarding flavivirus detection, we expected to find infection since the mosquito species collected are known to be involved in human arbovirus transmission in the Neotropics. For instance, Ae. serratus 47 and Mansonia titillans 48 have been implicated as potential secondary vectors of yellow fever in Brazil and have been found susceptible and potentially involved in Mayaro virus transmission. Also, Ae. fulvus and Cq. nigricans have been found infected with Madariaga virus 16,49 . However from the 25 samples that tested positive in the initial PCR screen for flavivirus, we didn't detect any known virus. Some of the "false positives" may have been a combination of flavivirus targets combined with false targets in the blood component. Mixed amplicons are generally uninterpretable by Sanger sequencing, and require either cloning or next generation sequencing methods to resolve the resulting sequence ambiguities.
In conclusion, our work provides an important contribution in the field of disease ecology by identifying mosquito-vertebrate interactions including species with medical importance and establishing well-defined patterns of host selection by Colombian mosquitoes.

Methods
Study site. Our sampling site was the locality of San Juan, Carare, in the department of Santander, Colombia (06°43′N, 74°09 ′W; 170 to 210 m.a.s.l) (Fig. 4). The area is a flooded tropical rainforest located between the central and eastern Andes, in the ecoregion of Magdalena Medio River Valley. The region is known to support alphavirus and flavivirus transmission with high mosquito diversity 50 and the presence of sylvatic fauna, particularly primates 51 . This area, as most of the Magdalena river valley, has been deforested by agriculture and livestock activities, generating a matrix of primary humid forest fragments with shrub-like areas mixed with crops 52 . We sampled during six campaigns at this site between March 2017 and April 2018.

Mosquito and vertebrate collection and identification.
We collected mosquitoes by establishing four transects, each with three CDC miniature light traps, three BG-Sentinel traps baited with BG-lure, and three resting traps, separated by 10 m. All traps remained active from 18:00 h to 6:00 h during five days in each campaign. Additionally, Prokopak ® aspirators were used at the undergrowth level during two consecutive hours in the evening, two days in each of the six campaigns for a total sampling effort of 24 h of aspiration. After collection, we sacrificed the insects using ethyl acetate and sorted them, keeping all the Culicidae. We stored the mosquitoes in 1.5-mL Eppendorf® tubes by trap-night collection in liquid nitrogen, until transport to the Research Center in Tropical Microbiology and Parasitology (CIMPAT for its acronym in Spanish) at the Universidad de Los Andes. In the laboratory, we selected female mosquitoes with blood-engorged abdomens for molecular analyses, and the remaining specimens were kept for a separate analysis. We identified mosquito species according to external morphology using available taxonomic keys 53,54 , and species identity confirmation was accomplished by DNA barcoding. High Pure PCR Template Preparation Kit (Roche) was used for extraction of nucleic acid following the manufacturer's protocol. The DNA barcode was amplified from a 658-bp region of the mitochondrial COI gene 55 . The PCR mixture was prepared to a final volume of 25 μL which contained 12.5 μL of 2 × GoTaq Green Master Mix (Promega), 10 μM of primers LCO1490 and HCO2198 and 5 μL of DNA. To determine the band size, amplification products were visualized on 2% agarose gels and stained with Sybr Safe (Invitrogen). The amplification products were sequenced bi-directionally at Gencore, Universidad de Los Andes. Raw sequences were assembled using DNAbaser and aligned with reference sequences from GenBank and the Barcode of Life database (BOLD). Pairwise nucleotide sequence divergence was estimated among all sequences using the K2P model implemented in MEGA X software. The platform ABGD software (Automatic Barcode Gap Discovery) was used to find values of intra and interspecific distance (barcode gap) for species delimitation 56 . Phylogenetic reconstruction was performed with PhyML 3.0 57 using the GTR + I + G nucleotide substitution models suggested by the SMART model selection and branch support was evaluated with the aLRT.
Additionally, to have an estimate of vertebrate species richness and relate it to blood sources found in mosquitoes, three transects each with 20 Sherman traps and 10 Tomahawk traps were set, for a total capture effort of 90 traps/night. Additionally, bats were surveyed using four mist nets, from 6:00 pm to 10:00 pm during eight consecutive nights during each of the six field trips to the study region. An additional field trip of five days for bird survey was conducted (June, 2017); during this survey, we conducted ten point counts for 20 min each, distanced by 250 m apart, during the peaks of activity (close to sunrise and sunset). www.nature.com/scientificreports/ Blood-meal analysis. To detect vertebrate blood sources, we removed the abdomen from each engorged female mosquito with blood content greater than 60% and extracted nucleic acid with DNA/RNA MiniPrep Viral kit (Zymo, Inc.) following the manufacturer's instructions. We used conventional PCR with universal primers for the amplification of a region of the mitochondrial cytochrome b (Cytb) gene of birds and mammals 58,59 . For samples that did not amplify, we used primers for amplifying a region of the mitochondrial COI gene 2 . PCR was performed in a reaction containing 12.5 μL of 2 × GoTaq Green Master Mix® (Promega), 5 μM of each primer and 2.5 μL of DNA. The amplification products were visualized on agarose gels, and positive samples were sequenced by capillary electrophoresis (the Sanger method) using an ABI-3500 genetic analyzer from Life Technologies (3500 Genetic Analyzer) at the Gencore Laboratory. We edited the raw sequences using the DNA Baser software (Heracle BioSoft SRL) and compared with sequences available in the GenBank (National Biotechnology Information Center) and Barcode of Life Data System (BOLD) databases through the BLAST platform using the Blastn algorithms, and when possible, confirmation with Blastx. Samples were identified accepting homology values higher than 96%.

Flavivirus detection.
To detect infected females we first tested pools of RNA derived from four individual mosquito abdomens (2 μl each), and when a positive pool was found, the individual RNA samples were retested individually (8 μl). This one-step RT-PCR assay utilized the Quantitect SYBR Green RT-PCR kit (Qiagen) following manufacturer's recommended protocol with primers PF1S and PF2R-bis, which recognize a 260 bp target in the NS5 gene of known flaviviruses 60 . Nucleic acid was extracted manually using the DNA/RNA Viral Mini-Prep kit (Zymo, Inc.). Amplification of a viral target sequence was considered successful when the fluorescence signal indicated a peak with a melting temperature between 75 and 85 °C. Prevalence was calculated as the number of positive samples over the number of processed samples (%). Flavivirus-positive females were evaluated for Zika virus infection in a second real-time RT-PCR following the protocol previously described 61 . Each reac-

Network analysis.
To evaluate if mosquitoes show specific blood-source selection, or if they feed randomly on available vertebrates, a network analysis was performed. First, an interaction-weighted matrix was built, using mosquito species in rows and vertebrates found through blood source identification in columns. In the matrix, the values used for species interactions corresponded to the number of identified bloodmeals derived from each vertebrate species for each mosquito species. For bipartite network visualization, we used the software GEPHI version 0.9 62 . The network was derived by using the force-directed algorithm FORCEATLAS2 and posterior manual editing. After creating the interaction network, values were calculated for quantitative modularity, complementary specialization and nestedness. Modularity recognizes unexpected species clusters, i.e. species interacting more frequently than expected by chance encounters, based on local availability of vertebrate species.
To calculate the modularity, Q, we used the QuanBiMo algorithm with "Beckett" method 63,64 . This algorithm uses quantitative interaction strengths across the community to portray the structure via maximizing weighted modularity. To incorporate stochasticity, we performed 50 iterations with 1,000,000 steps. We selected the command "metaComputeModules" to obtain the maximum possible value and set the modularity calculation. The statistical significance of Q can be assessed with the z-score. A network with a z-score above two is considered significantly modular. Complementary specialization occurs when a mosquito, or group of mosquitoes, feeds selectively on a specific vertebrate species (or group of vertebrates) and is measured using the H2 index. Lastly, the topology of the network at the community level was evaluated with the nestedness metric (NODF index). This index measures the degree of specialist mosquito species interacting with vertebrates that are also used by generalist mosquito species. We calculated the NODF index, with the function "nestednodf " included in the "Vegan" package, following methods previously described 65,66 . A value of 0 indicates non-nestedness, while 100 means complete nestedness 67 .
To identify the ecological traits at vector species level within the network, species specialization, specificity and strength indexes were also calculated. Related to the Shannon diversity index, species specialization (d) reflects the number of vertebrate species utilized for bloodmeals, and ranges from 0 (generalist) to 1 (specialist) 68 . It is calculated with the coefficient of variation of interactions (COV), which is a measure of specificity as described 69 . COV values range between 0 and 1 as minimum and maximum variation, indicating low to high specificity. Lastly, a strength metric (Species Strength) was used to evaluate the species richness (number of vertebrate species) utilized for bloodmeals by each mosquito species 67 , calculated as previously described 70 . To identify the most relevant vertebrate host species within the network, betweenness centrality (BC) was calculated. The BC index describes the importance of a node as a connector between different parts of the network. Nodes with BC > 0 connect areas of the network that would otherwise be sparsely or not connected at all. This BC index has been used to detect potential bridge hosts in other host-parasite networks 40,44 . All the network analyses were performed using the Bipartite package version 2.02 for R version 3.2.