This study was approved by the Beijing Municipal Committee of Animal Management before sample collection. All experiments were performed in accordance with the approved guidelines and regulations under approval number #SYSU-IACUC-MED-2021-B0123.

This study was conducted from July to August in 2014 at Shennongjia Forestry District natural reserves (SR), a zoo located in Beijing (BR) and a wildlife zoo located in inner Monglia (ER). All fecal specimens of rhesus monkeys (Macaca mulatta) were collected following defection at three sampling locations: one wild (SR), one semi-captive (ER), and a zoo-captive population (BR). Details of the sample collection sites, sample groups and food usage are presented in Table 1. DNA extraction of the fecal samples was performed using the TruSeq™ DNA Sample Prep Kit (Illumina) following the manufacturer's instructions.

Fecal samples from each monkey were subject to 16S rRNA amplicon sequencing. The V3–V4 hypervariable regions of the bacterial 16S ribosomal RNA (rRNA) gene were amplified using barcoded primers, 341F- 5′-CCTACACGACGCTCTTCCGATCTN (barcode) CCTACGGGNGGCWGCAG-3′ and 805R-5′-GACTGGAGTTCCTTGGCACCCGAGAATTCCA (barcode) GACTACHVGGGTATCTAATCC-3′, according to the Illumina 16S Metagenomic Sequencing Library Guide. The amplicons generated were sequenced on an Illumina HiSeq platform in a 2 × 250 paired-end mode. All sequencing and library preparation procedures were performed by Sangon Biotech (Beijing, China).

The raw amplicons generated were screened, trimmed, filtered, denoised, and chimera-depleted using QIIME2 version 2018.2 (http://qiime.sourceforge.net). Short, ambiguous sequences and chimeras deriving from the PCR process were removed using DADA2 plugins. Sequences were clustered into Operational Taxonomic Units (OTUs) and then assigned to bacterial sequences with at least 99% similarity to representative sequences from the SILVA 132 database (http://www.arb-silva.de/). For statistical analysis, all the sequences were rarefied to 1,112 reads for the downstream analysis. For each sample, the relative abundance of each bacterium identified was expressed as the percentage of total reads. QIIME2 was applied to profile the taxonomy of microbial composition in each group and to calculate alpha diversity matrices (including ACE, Shannon diversity index and Simpson index) (13, 14). To evaluate the variation between different groups, beta-diversity distance matrices (including Bray-Curtis distances, weighted and unweighted UniFrac values) were performed using rarefied data sets, and subsequently principal coordinate analysis (PCoA) was conducted to visualize the dissimilarities in the fecal bacterial communities among different groups of rhesus monkeys.

All reads from the high-throughput DNA sequencing data were mapped to reference genomes of Macaca mulatta (NCBI txid:9544) using Bowtie2 (15) to remove genetic material of host origin after quality-trimming by Trimmomatic (16). To profile the bacterial results from microbial composition, CCMetagen (17) was used for taxonomic annotation against nt database. To screen for viruses, host-filtered reads from the metagenomic sequencing data sets were assembled using MEGAHIT (18) then compared against the entire nr database in GenBank using Diamond BlastX e value <10⁵) (19). Any viral reads and contigs identified by Blast were then extracted and reassembled using the assembler implemented in Geneious v.11. This process identified abundant adenovirus sequences that were then reassembled into an entire adenovirus genome. This genome was then translated into amino acid sequences for gene annotation and functional prediction using Conserved domain databases (CDD).

The assembled sequences were then aligned using the MAFFT version 7 with implemented E-INS-I algorithm (20). Conserved domains within the E1A and 100k protein of adenoviruses were used for subsequent phylogenetic analyses. After removing all ambiguously aligned regions using TrimAl (21), the final lengths of E1A and 100K protein alignments were 832 and 1,379 amino acid residues, respectively. Phylogenetic trees of these data were inferred using the maximum likelihood method (ML) implemented in PhyML version 3.0, employing a Subtree Pruning and Regrafting topology searching algorithm. Statistical support for specific groupings in the tree was assessed using the approximate likelihood-ratio test (aLRT) with a Shimodaira-Hasegawa like procedure with 1,000 replicate bootstrap. The phylogenetic trees were visualized using the FigTree program (http://tree.bio.ed.ac.uk/software/figtree).

To determine the presence of putative antimicrobial resistance (AMR) genes in the data, we analyzed the shotgun sequencing data using the KMA program (22) combined with the ResFinder reference database (23). To reduce false-positive results, genes were only considered in downstream analyses when p-values for the conclave score were lower than 0.05, only two genes were excluded due to their p > 0.05 (22). We also excluded the blaTEM116 gene which has been previously identified as a common laboratory contaminant (24). AMR diversity and abundance was visualized in R with the package ggplot2.

The 16S rRNA amplicon sequencing generated a total of 2,572,794 reads and 2,680 OTUs. The total number of raw reads across all groups from the high-throughput shotgun DNA sequencing data was 1,425,675,194. Rarefaction curves showed a similar trend in all three populations. Observed numbers of OTUs (Observed_OTUs), an indicator of alpha diversity, is a qualitative measure of community richness. By this metric, the population of SR (wild monkeys) harbored the highest numbers of OTUs among all three groups under the same sequencing depth (orange: SR; blue: BR; cyan: ER) (Supplementary Figure 1A). The Shannon-Wiener curves showed that the samples from all groups had plateaued (Supplementary Figure 1B). The rarefaction curves indicate that sequencing depth was sufficient to capture the bacterial diversity in all samples (Supplementary Figure 1).

Based on the OTU data, we examined the bacterial richness and diversity of captive (BR), semi-captive (ER), and wild (SR) groups using ACE, the Shannon index and the Simpson index (Figure 1). The number of OTUs identified in the samples depicted species richness, as estimated by ACE. A non-parametric Kruskal-Wallis test was performed in all groups. The richness indices (ACE) revealed no significant difference (p > 0.05) between the wild and zoo-captive groups of macaques (Table 2). However, bacterial diversity was significantly different (p < 0.05) among all groups, as evaluated with the Shannon and Simpson indexes. Furthermore, Shannon indexes revealed significant differences between the captive (BR) and semi-captive (ER) groups, whereas no differences between wild (SR) and semi-captive groups (ER) were found by any of the methods.

Principal coordinate analysis (PCoA) was performed based on unweighted UniFrac (Figure 2A) and Bray-Curtis distances (Figure 2B) to visualize the dissimilarities in the bacterial communities among different groups of monkeys. The unweighted UniFrac analysis provided a much stronger clustering by population than either the weighted UniFrac or Bray Curtis distances, indicating that the clustering is likely driven by presence or absence of key taxa in different populations, rather than by shifts in the ratios of dominant members of the microbiota. In addition, PCoA plots based on Bray–Curtis distance matrices revealed that the samples from different locations formed distinct clusters, indicating that bacterial community composition conforms with the groups they were in, and hence that there were clear differences among wild, captive, and semi-captive monkeys. Analyses of distances based on relative abundance showed semi-captive groups overlapped more with captive group than with the wild group. We additionally performed Permutational Multivariate Analysis of Variance Using Distance Matrices (PERMANOVA) based on unweighted-uniFrac dissimilarity matrices (Supplementary Figure 2). Accordingly, the PERMANOVA results, indicated that (p = 0.001, number of permutations is 999) higher pseudo-F value in comparison of SR and ER groups with others (Supplementary Table 1).

Based on 16S rRNA sequencing, the clustered operational taxonomic units identified in fecal samples were assigned to 32 bacterial phyla. Both 16S and WGS approaches identified Firmicutes, Bacteroides, Proteobacteria, and Actinobacteria as the most abundant phyla in all samples, although the proportion of Bacteroides and Proteobacteria differed substantially (Figures 3A,C). At the class level, the two approaches revealed different bacterial compositions (Figures 3B,D). The main differences lie in the Epsilonproteobacteria that only appeared at high abundance in the wild (SR) group from metagenomic sequencing, but not in the corresponding group from 16S sequencing. Furthermore, the proportion of the class Bacilli also varied greatly between the two approaches. In general, the 16S sequencing resulted in relatively consistent results across three groups, whereas WGS sequencing revealed relatively high levels of variation. Further comparisons were performed at the family level using 16S sequencing results for microbial composition between the three groups (i.e., wild, semi-captive and zoo-captive) (Figure 4).

To assess the adenoviral reads and contigs identified from group ER, a near complete genome was derived from reassembled reads sequences that were mapped to a reference adenovirus genome (a double-strand DNA virus). The total length was 34,291 nucleotides with a GC content of 56.9%. To further characterize the adenovirus, phylogenetic trees were estimated based on a sequence alignment of the conserved region of the E1A and the 100 K protein, and utilizing reference adenovirus sequences downloaded from NCBI/GenBank. The novel virus shared 70.5% (E1A) and 88.7% (100 K) sequence identity with the closest relative—Simian adenovirus 3—within the Simian adenovirus clade (Figure 5). Based on its level of sequence divergence, the newly discovered virus likely represents a new virus species that we tentatively termed “simian adenovirus ER” (GenBank accession number: MZ062897).

A total of 67 acquired AMR genes were detected in the DNA-seq data sets, representing resistance against nine classes of antibiotics (Figure 6). Genes providing resistance to aminoglycosides, beta-lactams, MLS (including macrolides, lincosamides, streptogramin) and tetracyclin were found across all locations tested (Supplementary Table 2). The semi-captive group (ER) showed highest variety and abundance of antibiotic genes, followed by captive group (BR). Diversity measures indicate the number of AMR genes detected against the ResFinder database in each class. Abundance was calculated as the sum of Reads Per Kilobase of each class of AMR maker per Metagenome (RPKM) in each library. Accordingly, the wild group(SR) had the lowest variety and abundance, while genes conferring resistance to sulphonamide, rifampicin, trimethoprim and phenicol were only detected in ER group (Figure 6). In both zoo-captive groups, genes conferring resistance against Tetracycline had the highest relative abundance [i.e., AMR genes abundance (RPKM)/Total AMR genes abundance (RPKM)] at 68% in ER and 88% in BR, while genes conferring resistance against Vancomycin were the most abundant in the SR group (97.5%).