Mosquito-vertebrate Interaction Networks in the Magdalena Medio Region of Colombia

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. Mosquito genomic DNA/RNA was extracted for vertebrate host identication 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 aviviruses 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 identied 26 vertebrate species, the majority of which were mammals (N=16; 61.5%). No avivirus infection was conrmed. 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 specic.


Introduction
Transmission of mosquito-borne pathogens depends on the vertebrate host selection of mosquitoes when seeking bloodmeals for their egg development. In ecosystems with anthropogenic disturbances, de ning vector feeding patterns (specialist vs opportunistic), is key to assessing potential changes in pathogen transmission cycles and the risk of pathogen spillover [1] [2][3] [4]. Some mosquito species are generalists and often show opportunistic feeding behavior such as those in the Culex genus, which feed on a diverse range of species [5][6]. Others are specialists, such as Aedes aegypti who feeds primarily on humans [7]. In the sylvatic-urban interface, humans, mosquitoes, and wild animals co-occur, increasing the risk of zoonotic disease spillover to humans [8]. Additionally, habitat disturbance by humans may facilitate the dispersion of anthropophilic mosquito species, changing vector-host interactions and driving increased contact between human-biting mosquitoes and wild reservoirs of zoonotic pathogens [9]. Opportunistic mosquito species can act as bridge vectors between wild and urban transmission cycles [10], increasing the risk of spillover from sylvatic to urban environments and vice versa [11]. Thus, understanding the interactions among mosquitoes and the vertebrates on which they feed is of great importance [1][2] [3]. Mosquito-borne viruses associated with sylvatic and urban transmission cycles, constitute an important challenge for the design of prevention and control strategies because they usually involve a diversity of vector and reservoir species in complex multi-pathogen multi-host dynamics [11] [12]. For example, West Nile virus (WNV) can replicate in a great variety of Culicine vectors and vertebrate hosts and can affect hundreds of species from various classes [13] [14].
In Colombia, many pathogenic arboviruses are enzootic/endemic with many mosquito species potentially involved as vectors [15] [16] [17]. However, data from ecological studies on host selection by those mosquitoes of medical importance are scarce [6]. This study aims to evaluate the patterns of host selection by mosquitoes, in relation to the availability of vertebrates using network analysis.

Results
Mosquito diversity: From six eld campaigns, 292 blood-engorged female mosquitoes, representing eight genera and 20 species, were collected ( Table 1). The most abundant species was Aedes fulvus (15.8%), and the most diverse genus collected was Culex, with six species identi ed, followed by Uranotaenia with four species. By analyzing the obtained Cytochrome c Oxidase Subunit I (COI) sequences, we con rmed the identity of 10 species, supporting the morphological identi cation. Genetic distances, as determined by the Kimura-2-parameter (K2P) for intraspecies values, were under 2.4% for all species (see Supplementary Fig. S1). The maximum likelihood (ML) tree, built using the K2P distances between species, is shown in Fig. 2. Con rmation of species identity through barcodes was achieved for nine individuals with genetic identities to reference sequences above 97.6%. Reference sequences for Culex adamesi were not available in any published database, so we contributed the rst sequences for this species. Also, species belonging to the genus Uranotaenia could not be morphologically identi ed and our barcodes showed only 92% of homology with Uranotaenia sapphirina, the closest species in the barcode database. within the pools were positive for the detection of the avivirus NS5 gene. Retesting these 25 RNA samples for detection of the Zika virus envelope gene yielded negative results. The cDNA obtained in each reaction was revealed by agarose gel electrophoresis, but we did not obtain unique bands or viral sequences through Sanger sequencing in any of the samples. Sequencing the amplicons of the NS5 targets failed to detect known aviviral sequences.
Network analysis: The mosquito-vertebrate network comprised 41 different species interactions (Fig. 3a). Modularity was high (Q = 0.59, z = 13.4, p < 0.05), indicating signi cant specialization of mosquito choice of vertebrates utilized for bloodmeals. We obtained eight modules after 50 iterations (Fig. 4). 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 ve 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.
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 con rmed that mosquito-feeding behavior is not signi cantly 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. 3a).  [20]. The use of mitochondrial COI gene has become a great help in delimiting species within groups displaying ambiguous traits [21]; we contributed the COI sequences for 14 mosquito species from an unsampled ecoregion. The e ciency of species identi cation was supported by the barcoding gap analysis, which implies that if a gap exists, a cut-off value should be de ned allowing a clear distinction between the interspeci c and intraspeci c distances. The K2P intraspeci c and interspeci c values obtained here are similar to values found in other studies that include related species [22].
We identi ed 26 species of vertebrates among 154 mosquito bloodmeals, and found that only 10 of these species were observed during eldwork. Four species P. opposum, M. caucae, D. novemcinctus, and Z. brevicauda, were only detected trhough 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 identi ed 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 ndings, considering that despite the abundance of samples obtained, we did not nd 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 ndings, Mitchell et al. (1987) identi ed Mansonia titillans frequently feeding on domestic species including Phasianidae birds and Bovidae [30]. 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 di cult to monitor, limiting the understanding of innate behaviors of the vector [32] [33] [34] [35]. One very speci c 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 rst 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.
In our study, in most of the mosquito species in which human blood was identi ed, 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 (signi cant 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 antagonistic-relationship 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 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 de ned 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 Betweness 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 avivirus detection, we expected to nd 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 a 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 avivirus, we didn't detect any known virus. Some of the "false positives" may have been a combination of avivirus 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 eld of disease ecology by identifying mosquito-vertebrate interactions including species with medical importance and establishing well-de ned 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. 1). The area is a ooded 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 avivirus 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 identi cation: 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 10m. All traps remained active from 18:00 hrs to 6:00 hrs during ve 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 hrs of aspiration. After collection, we anesthetized 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 identi ed mosquito species according to external morphology using available taxonomic keys [53] [54], and species identity con rmation 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 ampli ed from a 658-bp region of the mitochondrial COI gene [55]. The PCR mixture was prepared to a nal volume of 25 µL which contained 12.5 µL of 2x GoTaq Green Master Mix (Promega), 10 µM of primers LCO1490 and HCO2198 and 5 µL of DNA. To determine the band size, ampli cation products were visualized on 2% agarose gels and stained with Sybr Safe (Invitrogen). The ampli cation products were sequenced bidirectionally 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 nd values of intra and interspeci c 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 eld trips to the study region. An additional eld trip of ve days for bird survey was conducted (June, 2017); during this survey, we conducted ten point counts for 20 minutes each, distanced by 250 m apart, during the peaks of activity (close to sunrise and sunset).
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 ampli cation 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 2x GoTaq Green Master Mix® (Promega), 5 µM of each primer and 2.5 µL of DNA. The ampli cation 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, con rmation with Blastx. Samples were identi ed accepting homology values higher than 96%. Network analysis: To evaluate if mosquitoes show speci c 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 identi cation in columns. In the matrix, the values used for species interactions corresponded to the number of identi ed 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 signi cance of Q can be assessed with the z-score. A network with a z-score above two is considered signi cantly modular. Complementary specialization occurs when a mosquito, or group of mosquitoes, feeds selectively on a speci c 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, To identify the ecological traits at vector species level within the network, species specialization, speci city and strength indexes were also calculated. Related to the Shannon diversity index, species specialization (d) re ects the number of vertebrate species utilized for bloodmeals, and ranges from 0 (generalist) to 1 (specialist) [68]. It is calculated with the coe cient of variation of interactions (COV), which is a measure of speci city as described [69]. COV values range between 0 and 1 as minimum and maximum variation, indicating low to high speci city. 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.