Eukaryotic gut community of the bat Myotis arescens in anthropized landscapes in Chile

Sampling sites and bat fecal collection

The study areas are located in the Metropolitan and O’Higgins regions of central Chile. These regions are characterized by agricultural landscapes composed of a variety of crops, including horticultural species, cereals, alfalfa fields, orchards, vineyards, in addition to small remnants of native vegetation (sclerophyllous forest or shrublands), exotic tree plantations (Pinus sp. or Eucalyptus spp.), and urban or semiurban areas. According to the Köppen climate classification, the region has a temperate pluviseasonal Mediterranean climate (Peel, Finlayson & McMahon, 2007), with most rainfall concentrated in the winter season, June–August. The mean annual precipitation is approximately 360 mm, and the mean annual temperature of 15 °C (Luebert & Pliscoff, 2017). Fieldwork was approved by the National Agency for Research and Development, ANID (formerly CONICYT), under the sponsorship of the University of Chile (project number: 3160188).

To assess the gut eukaryotic community composition, we analyzed bat fecal samples, a widely validated proxy in microbiome studies (Thomas, Clark & Doré, 2015). Samples were collected from four M. arescens colonies, located in separated landscapes. Each colony inhabited human-made structures as roost sites, with colony sizes ranging from 40–100 individuals. These landscapes varied in their degrees of anthropogenic disturbance, defined as the percentage of native vegetation remaining within a 5 km radius of each colony, along with land cover or use types (Fig. 1). This buffer area was established considering documented foraging ranges of comparable insectivorous bat species with similar body sizes to those present in central Chile. Land cover classification was performed through manual digitization of 2017 Google Earth imagery at a 1:15,000 scale using QGIS software. Ten habitat categories were identified: riparian habitats, annual crops, orchards, water bodies, prairies, vineyards, rural areas, urban areas, hedgerows (linear vegetation barriers between fields), and native vegetation comprising sclerophyll forests and shrublands. This classification covered the dominant vegetation types in the study area. Spatial analysis was conducted using Patch Analyst 5.2 extension in ArcGIS 10.4.1. Detailed characteristics of the sampled landscapes are provided in Table 1.

Each colony was visited in March 2017. For each colony, three pooled fecal samples were collected, each pooled sample containing an average of 20 pellets. Feces were collected directly from the roost floor using sterile forceps. To ensure the freshness of bat fecal samples, we cleared the collection areas prior to the sampling period. This allowed us to confirm that any droppings found afterward were recently deposited. Samples were collected early in the morning, immediately following the bats’ overnight activity, further increasing the likelihood of obtaining fresh material. Only fresh samples were collected and all were stored in sterile Eppendorf tubes and kept at −20 °C (within 10 h of collection) to preserve the DNA prior to analysis (Hale et al., 2016).

DNA extraction and library preparation

DNA was extracted using the PowerFecal^® DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA) following the manufacturer’s instructions, with the exception that the DNA was eluted in 25 μL of Milli-Q water.

To amplify one portion of the 18S rRNA gene, the primer combination 7f forward and 570r reverse was used (Noyce et al., 2016). Library preparation followed the Illumina 16S Metagenomics Sequencing protocol, adapted for 18S rRNA gene. In addition to the specific primers, the sequence was supplemented with the design outlined by Fadrosh et al. (2014), which includes a linker sequence, an index sequence, and a heterogeneity spacer. Amplicons were sequenced using 600-cycle kits with paired-end technology on an Illumina MiSeq sequencer (Illumina, San Diego, CA, USA). DNA extraction, library preparation, and next-generation sequencing were carried out at AUSTRAL-omics core research facility of the Universidad Austral de Chile (Valdivia, Chile). The 18S metabarcoding dataset has been submitted to the European Nucleotide Archive (ENA) under the accession number PRJEB78431.

Bioinformatic and statistics analysis

Raw FASTQ files from 18S rRNA sequencing were processed using the DADA2 pipeline (Callahan et al., 2016), accessible at https://benjjneb.github.io/dada2/tutorial.html in R (Oksanen et al., 2013). This software performs fast and accurate inference of amplicon sequence variant (ASVs) with a single-nucleotide resolution, improving the precision of community composition analyses in metabarcoding studies.

Forward and reverse reads were truncated at 250 bp using the truncLen parameter. Sequences containing ambiguous bases (maxN > 0) or with expected errors exceeding 8 (maxEE = 8) were removed. A minimum quality threshold of Q = 0 (truncQ = 0) was applied, and PhiX control sequences were filtered out (rm.phix = TRUE). Filtered reads were used for ASV inference, and chimeric sequences were identified and removed using the removeBimeraDenovo() function with the consensus method.

Taxonomy assignment was performed using the SILVA v.132 database (Quast et al., 2013). ASVs assigned to the phyla Vertebrata, Arthropoda, or Phragmoplastophyta phyla were excluded from subsequent analyses. For alpha diversity analysis, the Shannon, Observed Richness and Simpson indices were calculated using the diversity function from the vegan package in R (Oksanen et al., 2013). To test for differences between sites, we first applied Kruskal-Wallis tests after assessing the data for normality using the Shapiro-Wilk test. In addition, we fitted generalized linear models (GLMs) to each diversity metric: a quasi-Poisson GLM was used for Observed richness to correct for overdispersion, and Gaussian GLMs were used for Shannon and Simpson indices. Rarefaction curves were generated using the ggrare function from the microbiomeSeq package in R. This function plots species richness (observed ASVs) as a function of sequencing depth across samples. Confidence intervals were displayed using standard error (se = TRUE), and samples were color-coded by sampling site.

To evaluate beta diversity and patterns of community composition across sampling sites, we performed non-metric multidimensional scaling (NMDS) analyses using two dissimilarity indices: Bray–Curtis (which considers relative abundance) and Jaccard (based on presence–absence data). For each index, we computed a dissimilarity matrix and performed NMDS using the metaMDS() function from the vegan package in R, with two dimensions (k = 2) and 100 random starts (trymax = 100) to ensure convergence on a stable solution. Stress values were used to assess the goodness of fit.

To identify environmental variables significantly associated with the NMDS ordinations, we applied the envfit() function from vegan, which fits environmental vectors onto the ordination space and uses permutation tests (n = 999) to assess significance. This analysis was applied to both Bray–Curtis and Jaccard ordinations.

Differences in community composition between sites were tested using permutational multivariate analysis of variance (PERMANOVA) via the adonis2() function, also from vegan, using both Bray–Curtis and Jaccard dissimilarity matrices. We used 999 permutations and evaluated the proportion of variance explained by the factor Site in each case.

The environmental matrix included site-level landscape characteristics recorded at each sampling location: riparian vegetation (Rip), natural vegetation (NaV), anthropogenic cover (AnC), orchards (Orch), water bodies (Wat), pastures (Pra), rural areas (Ru), urban cover (Urb), hedges (Hed), and vineyards (Vin). These variables were tested for correlations with gut eukaryotic community structure.

To summarize environmental gradients and reduce collinearity among landscape variables, we performed a principal component analysis (PCA). The first two principal components (PC1 and PC2) were used to visualize site separation and interpret dominant environmental gradients. These components were incorporated into vector fitting analyses (envfit), PERMANOVA models, and generalized linear models to assess their association with gut eukaryotic community diversity, using the vegan and stats packages in R.

A heatmap of relative abundance percentages by site was constructed using amp_heatmap, and boxplots were created with amp_boxplot, both from the ampvis2 package in R.

Phylogenetic analysis of Apicomplexa ASVs

Amplicon sequence variants (ASVs) assigned to the phylum Apicomplexa were extracted from the phyloseq object. Sequences were aligned with the AlignSeqs() function from the DECIPHER package, and the alignment was converted to a phyDat object. A preliminary phylogenetic tree was inferred using the Neighbor-Joining method and optimized under the GTR + I + Γ model with the phangorn package. The final tree was visualized with ggtree, and ASVs were labeled with genus-level taxonomic annotations. To illustrate ASV distribution across sampling sites, proportional colored circles were added to each tip, representing relative abundance and site of detection.

Raw data analysis and gut eukaryotic diversity

A total of 71,995 raw Illumina sequence reads were obtained from 12 samples. After quality filtering and chimera removal, 58,585 high-quality sequences remained. To minimize potential contamination from host or diet-derived DNA, ASVs assigned to the phyla Vertebrata, Arthropoda, and Phragmoplastophyta were excluded. The final dataset used for downstream diversity analyses consisted of 37,064 sequences distributed across 211 ASVs.

Rarefaction analysis revealed that samples of Site 2 and Site 3 harbored the highest eukaryotic richness, requiring greater sequencing effort to capture their diversity. In contrast, Site 1 and Site 4 showed lower richness and reached saturation with fewer sequences. The plateauing of all curves indicates that sequencing depth was generally sufficient to capture species richness across sites (Fig. S1).

No significant differences in alpha diversity were detected among sites based on Kruskal-Wallis tests for Observed richness (χ² = 3.65, p = 0.302), Shannon index (χ² = 6.08, p = 0.108), or Simpson index (χ² = 6.05, p = 0.109). These findings were further supported by generalized linear models: the quasi-Poisson GLM for observed richness did not reveal significant effect of site, and the Gaussian GLMs for Shannon and Simpson indices also indicated no statistically significant differences. A consistent, though non-significant, trend toward higher diversity was observed in Site 3 across both Shannon and Simpson indices (Fig. 2). These results suggest a relatively homogeneous distribution of gut eukaryotic alpha diversity across the studied sites.

Native vegetation percentage is the main factor explaining the eukaryotic gut community structure between sites

To assess differences in community composition across sampling sites, we conducted a permutational multivariate analysis of variance (PERMANOVA) using the Bray–Curtis dissimilarity matrix. Results indicated a significant effect of site on community structure (R² = 0.52, F = 2.90, p < 0.001), indicating that 52% of the variance in community composition was explained by location. Using presence–absence data and Jaccard dissimilarities, a second PERMANOVA also revealed significant differences in community composition among sites (R² = 0.41, F = 1.83, p < 0.001). This indicates that 41% of the variation in species composition was attributable to site. The ‘envfit’ function revealed that environmental variables significantly associated with community composition differed slightly depending on the dissimilarity index used. In NMDS based on Bray–Curtis distances, variables such as natural vegetation (NaV), anthropized cover (AnC), water bodies (Wat), and riparian vegetation (Rip) showed strong associations with the ordination (p < 0.001, p < 0.001, p < 0.001, and p = 0.011, respectively) (Fig. 3). In the Jaccard-based NMDS (Fig. S2), similar variables were significant, but orchards (Orch, p = 0.001), prairies (Pra, p = 0.003), urban areas (Urb, p = 0.002), and water bodies (Wat, p < 0.001) exhibited higher r² values, suggesting that these factors may be more influential in shaping species presence/absence patterns rather than abundance. In contrast, vineyards (Vin), which were significant in the Bray–Curtis ordination (p = 0.019), were only marginally significant under Jaccard (p = 0.057).

Principal component analysis (PCA) was used to summarize variation in site-level environmental variables (Fig. S3). The first two principal components (PC1 and PC2) explained major gradients in landscape composition and were significantly associated with gut eukaryotic community structure. Vector fitting analysis on the PCA ordination revealed that both PC1 (r² = 0.90, p = 0.001) and PC2 (r² = 0.59, p = 0.019) were significantly correlated with community composition. Similarly, PERMANOVA models using PC scores also revealed significant associations with beta diversity, both for Bray-Curtis (PC1 p = 0.001, PC2 p = 0.011) and Jaccard distances (PC1 p = 0.001, PC2 p = 0.001).

In contrast, generalized linear models using PC1 and PC2 as predictors did not reveal significant associations with alpha diversity indices, including Observed Richness (PC1 p = 0.797, PC2 p = 0.577) and Shannon diversity (PC1 p = 0.127, PC2 p = 0.758). Sites clustered in accordance with these environmental gradients, with Site 4 strongly associated with AnC and Vin, and Site 3 with higher NaV values (Fig. S3).

High representation of Apicomplexa and other ecologically relevant eukaryotes in bat fecal samples

For visualization purposes, only ASVs assigned at the phylum level were included. At the phylum level, only Apicomplexa, Ascomycota, and Basidiomycota were present at all the sampling sites. Among the total reads, Apicomplexa was the most abundant, although its relative abundance varied by Site: Site 1 had 86.1% of the total reads assigned to this phylum, followed by Site 2 (75%), Site 3 (37.6%), and Site 4 (16.2%) (Fig. 4).

The phylum Ascomycota was particularly abundant at Site 3 (41.9%), and Basidiomycota was more abundant at site 4 (33.5%). Platyhelminthes and Zoopagomycota were primarily detected at Site 3, whereas Euglenozoa were present only at Site 4 (9.9%), and Dinoflagellata were present only at Site 1 (0.7%). On the other hand, Mucoromycota was present at Site 2 (1.1%) and Site 3 (0.4%) (Fig. 4).

At the genus level, Eimeria from the Apicomplexa phylum was the most abundant genus found in our samples (Fig. S4).

Distribution of Apicomplexa across sampling sites

A total of 14 amplicon sequence variants (ASVs) belonging to the Apicomplexa phylum were detected across sites (Fig. 5). These ASVs included members of the genera Eimeria, Cryptosporidium, Cystoisospora, and Adelina, along with several taxonomically unclassified lineages.

The most abundant ASVs corresponded to the genus Eimeria (ASV17 and ASV57), were present in multiple sites, suggesting a relatively wide distribution and possibly a high prevalence among host organisms. In contrast, ASV111 (Cryptosporidium) and ASV109 (Cystoisospora) showed more restricted distributions, appearing only in a subset of sites. This could reflect ecological specialization or lower prevalence within the sampled communities.

A substantial proportion of the detected ASVs were unclassified, indicating either limited representation of these taxa in current reference databases or the presence of potentially novel Apicomplexan lineages. Notably, several unclassified ASVs were found across all sites, underscoring the need for further taxonomic and functional characterization. Of particular interest is ASV129, assigned to the genus Adelina, a group rarely reported in metagenomic studies.

Special attention is drawn to Site 3, this site harbored a broader range of ASVs, including both classified and unclassified taxa, suggesting a more complex parasitic community or a richer diversity of hosts. Conversely, Site 4 showed the lowest diversity, with fewer ASVs detected overall, potentially reflecting a simpler community structure or more constrained ecological conditions.

These results highlight the spatial heterogeneity of Apicomplexan communities across sites, shaped potentially by differences in host availability, environmental factors, or anthropogenic impacts.