Viral diversity of bat communities in human-dominated landscapes in Mexico

Using integrative epidemiologic techniques, we studied the changing relationships (beta and phylogenetic beta diversity) of multihost systems and virus associations in bat communities in fragmented landscapes from Chiapas, Campeche and Greater Mexico City. We combined computing applications, molecular detection, and nucleotide sequencing of coronaviruses, hantaviruses, paramyxoviruses and pegiviruses with ecological and phylogenetic analyses. A total of 22 viruses were discovered in 1,067 samples from 42 bat species, representing an estimated 78% of all viral richness in the system. Based on 17 virus genotypes discovered with an equal sampling effort, a total viral richness of 23 genotypes was estimated using a Chao2 statistic model. Using a residual model, we categorized host species and habitat types that are prone to harboring higher viral richness. Positive relationships were found between phylogenetic host diversity and both viral diversity (r = 0.41, p < 0.05) and viral richness (r = 0.51, p < 0.05). The beta diversity (the rate of change) of viral communities was explained by host beta diversity (r = 0.86, p < 0.05). To understand the change in viral and host communities, we partitioned beta diversity in nestedness (species loss) and turnover (compositional dissimilarity) components. In Chiapas, the host beta diversity was explained by the nestedness of species composition, while the phylogenetic host diversity was explained by turnover of the host lineages. Campeche showed a high phylogenetic host nestedness and low host turnover. Beta-diversity and beta-phylogenetic diversity indicated that patterns of local species assemblages and regional abiotic features in human-dominated landscapes are significant drivers of viral community composition. Our study represents the first effort in Mexico to study the relationship between viral diversity in bat communities in modified landscapes to understand host-virus relationships.


Introduction
Land use change appears to be the primary mechanism driving zoonotic diseases (Patz et al., 2004).The expansion of agricultural production and urbanization has simultaneously modified ecosystem structure and function, community structure (species assemblages), patterns of species distribution and biodiversity (Christian et al., 2009;Gibbs et al., 2009).These modified systems have produced suitable environments for multi-species interactions, particularly hosts, vectors, and/ or pathogens (McMichael, 2004;Rivard et al., 2007;Jones et al., 2013;Rubio et al., 2014).
To properly understand complex interactions in multihost systems, several ecological and phylogenetic tools have been used.In disease ecology, diversity indices have been used to correlate the number of species (richness) and the relative abundance of species present in a given community (alpha diversity) (Suzán et al., 2009) and in microbiome systems (Anthony et al., 2013a;Olson et al., 2014) with disease prevalence.Diversity indices have also been used to evaluate changes in host-parasite composition in host communities at the local, regional, and biogeographic scales (alpha, beta gamma diversity) (Svensson-Coelho and Ricklefs, 2011; Scordato and Kardish, 2014).From an evolutionary perspective, host and pathogen phylogenetic relationships have been studied, and diversity indices have been incorporated to measure changes in host species community assemblages through environmental gradients (Webb et al., 2002;Helmus et al., 2007).These phylogenetic methods offer additional dimensions to explore host-parasite interactions over time, such as host specificity, host-parasite co-evolution, host switching events, and phylogenetic barriers preventing pathogen transmission (Legendre et al., 2002;Streicker et al., 2010;Poulin et al., 2011).The study of ecological and phylogenetic interactions between host-pathogen systems integrates the role of environmental influences on host and pathogen distributions across time and spatial scales and across different levels of biological organization beyond taxonomic levels (Hawley and Altizer, 2011) In this study, we examined the relationship between host diversity and the diversity of four viral taxa in bats from human-dominated landscapes in Mexico.Two hypotheses were tested related to the effect of host species and host phylogenetic diversities on viral diversity and the influence of habitat type on the composition of host and viral communities.First, we hypothesized that (1) host communities with high species and phylogenetic diversities will support high values of viral diversity; (2) changes in host and viral community composition across a habitat type will be reflected in high values of beta diversity and phylogenetic-beta diversity.

Sample collection
Bats were captured at three different sites in Mexico: the Reserva de la Biósfera Montes Azules (RBMA) in Chiapas, the Reserva de la Biósfera Calakmul (RBC) in Campeche, and Greater Mexico City (GMC) including the Distrito Federal and Metropolitan Area.The first two sites, located in southeastern Mexico, represent regions of high species diversity and are characterized by large tracts of continuous primary vegetation, while the GMC site is highly urbanized with vegetation patches.A high evergreen forest characterizes RBMA, while RBC is dominated by tropical semi-deciduous forest; both regions have high anthropogenic pressure.In RBMA and RBC, bats were collected from three different habitat types: 'Forested' (Fd), where signs of human impact are largely absent and the original vegetation persists; 'Fragmented' (F), where areas of primary vegetation are interspersed with agricultural/ rangeland; and 'Disturbed' (D), the transition zone between areas of secondary vegetation and agricultural/rangeland or urban areas.In the GMC sites, bats were captured in 'Urban' (U), human-dominated areas and 'Fragmented' habitats.We used 5 mist-nets (each 9 x 3 m wide) that were opened at dusk and remained open for four consecutive hours.Each habitat was sampled once in six months.Bats were identified using a field guide (Medellín et al., 2008).The minimal distance in RBMA was 2 km, while in RBC it was 10 km.A mantel test was performed to ensure site independence due to the geographic distance (RBMA; r = 0.55, p =0.01; RBC r = 0.57 p = 0.006).Oral and rectal swabs and, when possible, blood samples were collected from each animal.Samples were collected in lysis buffer and preserved at -80°C until transfer to the Center for Infection and Immunity, Columbia University, New York for viral screening.

Virus discovery
A total of 1,067 samples from 608 individuals representing 42 bat species were tested for the five viral families/genera (Table S1).Total nucleic acid was extracted from all samples using the EasyMag® (bioMérieux, Inc Darham, NC, USA.) platform, and cDNA synthesis performed using SuperScript® III first strand synthesis supermix (Invitrogen), all according to the manufacturer's instructions.Viral discovery was performed using broadly reactive consensus PCR primers, targeting the L-Segment for hantavirus (HTV) detection (Klempa et al., 2006) and the polymerase (pol) gene for paramyxovirus (PMV) detection (Tong et al., 2008).PCR products of the expected size were cloned into the StrataClone™ PCR cloning vector and sequenced using standard M13R primers.CoV, hepacivirus (HPV) and pegivirus (PGV) detections have been previously reported, and these viral sequences were detected in the same 1,067 samples (Anthony et al., 2013b;Quan et al., 2013).

Estimates and Completeness of Viral Richness
We evaluated our sampling effort (the number of samples tested for a given virus), using two methods: by producing rarefaction and extrapolation curves and by calculating the values of the residuals of the linear regression between viral richness within a host and sampling effort by host.Rarefaction and extrapolation curves are statistical techniques to estimate the number of species for a given number of samples (Magurran, 2004;Chao and Jost, 2012), allowing the evaluation of the sampling effort and estimation of the number of host samples required to obtain a viral richness value with 95% confidence (Chao et al., 2014).We evaluated the viral richness, defined as unique viruses discovered in the 1,067 samples, by constructing sample size-based rarefaction and extrapolation curves using a three-fold original sample effort (3,201 samples) (Chao et al., 2014) with the R iNEXT library (Hsieh et al., 2013).The same methodology was used to explore viral completeness by habitat type.For this purpose, we only considered samples with the same number of PCR screenings (CoVs, PMVs and HTVs).To identify host species associated with higher viral richness, we used a methodology proposed by Herbreteau, 2012.We calculated residual values from the linear regression of the logarithm of viral richness and the logarithm of sampling effort for each species and at each disturbance level.These data were logarithmically transformed to stabilize the variance.Host species or disturbance levels with positive or negative residual values were identified as host species with more or less viral richness than expected by the regression model (Herbreteau et al., 2012).

Host and Viral Diversity
To study regional host and viral alpha diversities, abundance matrixes (host and virus) were constructed, where the rows were disturbance level and the columns were (i) host species and (ii) viruses discovered in each disturbance level.Using the R vegan library (Oksanen et al., 2013), a Shannon-Wiener diversity index (Shannon, 1948) was calculated for each matrix.Values ranged from 0, when there is only one species present, to 1, when all species are equally represented in the sample (Magurran, 2004).

Phylogenetic Diversity and Host Specificity
The mammalian super tree (Bininda-Emonds et al., 2007) was used to calculate the phylogenetic diversity (PD) of host communities using the R Picante library (Kembel et al., 2010).The PD was measured by calculating the sum of the total branch length of the host species phylogeny sampled in each habitat type (Faith, 1992).Because the data on host taxonomic diversity were not normally distributed, only phylogenetic analyses were performed.The relationships of viral richness and viral diversity with host PD was explored using a linear model.To quantify host-viral taxonomic associations, we used a modified index of host specificity proposed by Poulin and Mouillot, 2003, that measures the PD of host communities associated with each virus.Viruses with high values of host specificity exhibit plasticity to infect a wide range of hosts, while viruses with lower values are restricted to a few closely related host species (Poulin and Mouillot, 2003;Poulin et al., 2011).

Beta Diversity and Phylogenetic Beta Diversity
A Pearson correlation test was performed to explore the relationship between host BD and PBD with changing compositions of viral communities by habitat type calculated by the Sorensen index.To evaluate the change in composition of viral and host communities (BD) within regions, we used measures of beta-diversity: spa-tial turnover (ß SIM ) and nestedness (ß SNE ) components (Baselga, 2010).Spatial turnover (ß SIM ) measures the replacement of species by other species due to environmental factors or spatial isolation, such as by habitat fragmentation (Calderón-Patrón et al., 2012).The nestedness component (ß SNE ) measures whether sites with smaller number of species are subsets of richer sites (Ulrich et al., 2009).These components were calculated for taxonomic and phylogenetic beta-diversity analyses.Phylogenetic beta-diversity (PBD) measures how phylogenetic relatedness changes across space in the same manner that BD measures how species composition changes across space (Graham and Fine, 2008).The PBD between disturbance levels was obtained using the inverse of the PhyloSor index (Bryant et al., 2008).This index represents shared branches between communities from two sites.Values range from 0 when no species are shared to 1 when all species in the two locations are the same.The methodology of Leprieur (2012) was used to calculate the phylogenetic turnover (Pß SIM ) and to measure the phylogenetic dissimilarity, nested patterns of species assemblages (Pß SNE ).The functions beta.multi for BD and phylo.beta.multifor PBD from the R betapart library were also applied to calculate the influence of each component on host and viral community composition by habitat type in each region (Baselga and Orme, 2013).

Estimates and Completeness of Viral Richness
Based on the 17 virus genotypes discovered, we estimated a maximum richness of 23 genotypes using a Chao2 statistic model (Chao and Jost, 2012).The sampling effort of 1,067 samples represents a completeness of 81% in relation to the estimated viral richness.The rarefaction sample coverage function estimates 97% completeness with a sample size of 3,201 (three-fold sample size) (Fig. 2).The comparison between habitats showed the highest value of completeness (53%) in forested habitat, followed by disturbed (43%) and fragmented (15%) habitats.The estimates of viral richness with a three-fold original sample effort by habitat were 24 from fragmented habitats, 19 from forested habitats and 10 from disturbed habitats (Table 1).

Host and Viral Diversity
RBMA was the region with the most virus genotypes discovered (12), followed by RBC (11) and GMC (3) (Fig. 4).As we expected, RBMA was the most diverse region in terms of host species and host phylogeny, with Fd being the most diverse area in both diversity scales (H = 2.79, PD = 484.8).Interestingly, D presented high values of both species and phylogenetic diversities compared to F (Table 2).A

Phylogenetic Host Specificity
Of 22 viruses discovered, only six were associated with more than two host species, and 16 viruses were detected in one host species (Table S3).The flavivirus PgV was the virus with the lowest value of phylogenetic host specificity (114.5) and was detected in three species from three different Phyllostomidae subfamilies; Stenodermatinae and Phyllostominae, both from forested sites, and Glossophaginae, from disturbed habitats.Hanta 2 (87.8) was detected in two bat species from different habitats; Carollia sowelli (disturbed) and Trachops cirrhosus (forested), while the CoV MexCoV 5b was found in three bat species, two of which belong to the same genus; Artibeus jamaicensis, A. lituratus (forested) and C. sowelli (fragmented).MexCoV 11a (68.3),MexCoV 11b (68.3) and MexCoV 1 (65.6) were detected in two bat species from the same genus (Figure 1A).

Beta Diversity
We found that 86% (p < 0.05) of the beta diversity of viral communities was explained by the turnover of host species between habitats.The same trend was observed with host PBD and virus beta diversity but was not statistically significant (r = 0.79, p > 0.05).For RBMA, the community dissimilarity in host species composition through the habitat gradient was relatively high (ß SOR = 0.49) and is explained by the nestedness component (ß SNE = 0.57), as most of the bat species sampled in low-richness sites are also contained in the richest site (Fd= 26 species).PBD showed different behavior; although the value of overall PBD is similar (Pß SOR = 0.42), the phylogenetic composition change was partially explained by the turnover component (Pß SIM = 0.25) due to the low phylogenetic relationship contained in the species subset replaced by nestedness.The overall beta diversity in RBC was relatively low (ß SOR = 0.37), suggesting that the habitat type has a low influence on the species assemblage.The observed pattern in RBC was different from RBMA, while species composition dissimilarity was moderately driven by species replacement (Pß SIM = 0.22).In contrast to RBMA, the PBD analysis indicated a replacement dominated by highly phylogenetically related species (Pß SNE = 0.79).The interpretation of the pattern observed in GMC is limited, as only two communities were sampled and the nestedness component is impossible to calculate.In RBMA and RBC, we observed high overall values of beta diversity in viral communities (Table S3), and both were explained by the turnover component (ß SIM =0.66), suggesting that changes in habitat quality drive a high replacement of virus species regardless of host composition.

Discussion
In this study, we evaluated the relationship between bat diversity and the diversity of 4 medically important viral families within an environmental gradient in  As hypothesized, significant positive correlation of phylogenetic diversity with viral richness and viral diversity was detected, supporting both the habitat heterogeneity hypothesis (Lawton, 1983), which proposes a strong relationship between environmental diversity (in this case the phylogenetic host diversity) and biological diversity (viral diversity), and the keystone structure hypothesis (Tews, 2004), which refers to "keystone structures" as host species that provide distinct resources, which may be linked to different viral species.Considering bat hosts as habitats, the keystone structures can be categorized in characteristics that relate to reproduction and abundance (reproductive potential, longevity, trophic guild, abundance and adult mass), transmission potential (home range, diet breath and roost size) and phylogenetic relationships (phylogenetic distance and phylogenetic distinctiveness).Our findings suggest that factors such as fragmentation and habitat loss drive species assemblages, resulting in areas of greater risk for zoonotic disease emergence, as proposed by (Gay et al., 2014;Kamiya et al., 2014;Rubio et al., 2014).Future studies are required to identify which host traits determine the viral community assemblages, but evidently greater viral diversity does not imply greater health risks; in fact, the correlation between viral and host diversity suggests that only under The phylogenetic host specificity calculated in this study not only reflects the number of bat species infected by a single virus but also helps to explore the phylogenetic relationships among these species.Our findings show that few viral species possess high host plasticity, such as PgV and HTV 2, while most viruses detected showed high phylogenetic specificity.Due to the quality of our data, we cannot conclude that viruses found in only one bat species are strictly exclusive to them, and further studies are needed.However, the coronavirus family showed high host phylogenetic specificity at the genus level, as shown in Anthony et al., 2013b).
Changes in viral community composition through the evaluated anthropogenic landscapes showed a strong dependence on host species turnover; however, this relationship was not statistically significant when host PBD was considered.In general, different patterns between regions and between diversity levels were reported.For instance, RBMA is characterized by a nestedness process in host community composition, and changes reflect of a high PD host pool, while in RBC, the beta diversity values were relatively low due to low PD host diversity, and the PBD was explained by the nestedness component.As hypothesized, we found high values of beta diversity in viral communities, supporting the hypothesis of perturbation, where land use change modifies parasite dynamics in multihost systems by cross-species shifting in parasite transmission (Murray and Daszak, 2013).It has been widely proposed that habitat transformation drives the exposure of novel hosts to a rich pool of parasites, especially in high-diversity regions, influencing the cross-species transmission rate (Lloyd- Smith et al., 2009;Brearley et al., 2013;Murray and Daszak, 2013).
The use of beta analyses at both the taxonomic and phylogenetic scales has provided a useful tool to understand whether host species or environmental filters can determine parasite composition (Svensson-Coelho and Ricklefs, 2011; Scordato and Kardish, 2014).We could not demonstrate that host phylogeny determines the composition of viral communities due to spatial scale limitations in our study.Further studies are necessary to test the correlation between host phylogeographi-

Figure 1 .
Figure 1.Bipartite graph of 22 virus genotypes discovered in 17 bat species.A. Viral richness associated with bats.The width of the green boxes represents viral and positive hosts abundances.B. Viral richness associated with habitat type.

Figure 2 .
Figure 2. Rarefaction and extrapolation sampling curve based on 1,067 samples.Orange line represents accumulation curve of virus genotypes over samples tested.Green point, viral richness = 22 in 1,067 samples analyzed for virus tested.Solid green line: rarefaction curve, green dashed line: extrapolation-sampling curve.The numbers of samples needed to obtain the completeness percentages of 85, 90, 95 and 97% are presented.

Figure 3 .
Figure 3. Distribution of residual values from the linear relations (A) between host viral richness and sampling effort and (B) between habitat type viral richness and sampling effort.Host species and habitat type were reordered by residuals of the sampling effort regression.

Figure 4 .
Figure 4. Viral richness in the three regions of study: (A) Greater Mexico City, (B) Montes Azules and (C) Calakmul.Number of virus genotypes discovered per region is shown.

Table 1 .
Estimates of viral richness and sample completeness by habitat type.