2. Materials and methods
2.1 Study area and sample collection
FNNR is located in the transition from Yunnan Guizhou Plateau to the hills of western Hunan and includes the main peak of Wuling Mountains (27°49′50′′ to 28°1′30′′N, 108°45′55′′ to 108°48′30′′E). It is located at an altitude of 500–2,572 m, with an area of 419 km2. In terms of vegetation, it is composed of evergreen broad-leaved, deciduous broad-leaved, evergreen deciduous broad-leaved mixed, and needle- and broad-leaved mixed forests (Xiang et al., 2009). The FNNR is one of the areas with the highest forest biodiversity protection priority in the upper reaches of the Yangtze River. The natural vegetation in the protected area is relatively well preserved, providing a good habitat for various rare and endangered animals and plants, and this area is the only habitat for the wild population of R. brelichi (Wu et al., 2006).
From December 2021 to November 2022, fresh fecal samples were collected from all suspected R. brelichi in the northern part of FNNR. During this sampling, the surface soil was removed manually using disposable medical gloves. Furthermore, plant branches, leaves, and other parts attached to the samples were removed with disposable sterile tweezers and scalpels, and the middle uncontaminated part of the sample was selected and placed in 50-mL sterile centrifuge tubes, which were then marked and packed into a sealed bag. The geographical information about the site was recorded on the sealed bags (Figure 1). To prevent DNA degradation and microbial reproduction in the fecal samples, they were immediately placed in liquid nitrogen tanks after collection and were brought back to the laboratory for storage at −80℃ until use. Fecal samples from captive animals were collected simultaneously using the same method. These samples were collected from five R. brelichi in captivity at the rescue station of FNNR Administration. A total of 42 and 15 fecal samples were collected from animals in the wild and those in captivity, respectively.
2.2 DNA extraction, polymerase chain reaction (PCR) amplification, and sequencing
Mitochondrial CO Ⅰ-based DNA barcoding was used to identify 42 fecal samples collected in the wild, of which 31 belonged to R. brelichi and 11 belonged to other monkey species. Total genomic DNA samples were extracted using OMEGA DNA Kit (M5635-02) (Omega Bio-Tek, Norcross, GA, USA), following the manufacturer’s instructions, and stored at −20°C prior to further analysis. The quantity and quality of extracted DNA were measured using a NanoDrop NC2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and 0.8% agarose gel electrophoresis, respectively.
The V3–V4 highly variable region of the bacterial 16S rRNA gene, with a length of approximately 468 bp, was selected for PCR amplification. The primer sequences were as follows: forward primer 338F (5′-barcode+ACTCCTACGGGAGGCAGCA-3′, the barcode is a seven-base oligonucleotide sequence used to distinguish different samples in the same library) and reverse primer 806R (5′-GGACTACHVGGGTWTCTAAT-3′). PCR reactions were performed in triplicate 25 μL mixtures containing 0.25 μL of Q5 high-fidelity DNA polymerase, 5 μL of 5× Reaction Buffer, 5 μL of 5× High GC Buffer, 2 μL of 10 mmol L−1 dNTP, 1 μL of forward primer (10 μmol L−1), 1 μL of reverse primer (10 μmol L−1), and 2 μL of template DNA, with ddH2O added up to a final volume of 25 μL. Thermal cycling consisted of initial denaturation at 98°C for 5 min, followed by 25 cycles consisting of denaturation at 98°C for 30 s, annealing at 53°C for 30 s, and extension at 72°C for 45 s, with final extension for 5 min at 72°C (PCR instrument: ABIGeneAmp® 2720, USA). PCR amplicons were purified with Vazyme VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China) and quantified using the Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA). After the individual quantification step, amplicons were pooled in equal amounts, and paired-end 2×250 bp sequencing was performed using the Illumina NovaSeq platform with NovaSeq 6000 SP Reagent Kit (500 cycles) at Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China).
2.3 Data analysis
Microbiome bioinformatics was performed with QIIME2 (https://docs.qiime2.org/2019.4/tutorials/). Specific analysis was conducted as follows: (i) Raw sequence data were demultiplexed using the demux plugin followed by primer cutting with the cutadapt plugin (Martin, 2011). (ii) Sequences were then quality-filtered, denoised, and merged, and chimeras were removed using the DADA2 plugin (Callahan et al., 2016). The obtained sequences with 100% similarity were merged, and amplicon sequence variants (ASVs) and abundance data tables were generated. (iii) The QIIME2 classify-sklearn algorithm was used (https://github.com/QIIME2/q2-feature-classifier). Specifically, for the feature sequence of each ASV, QIIME2 software was used with the default parameters and species annotation was performed using a naive Bayes classifier that had been pre-trained to obtain taxonomic information corresponding to each ASV (Bokulich et al., 2018).
We evaluated the alpha diversity of the microbial communities using Chao1, ACE, Shannon, and Simpson indexes, as calculated with QIIME2 and visualized as box plots. Additionally, we performed Kruskal–Wallis tests to examine differences in alpha diversity between the two groups (Xi et al., 2023). Principal coordinate analysis (PCoA) and non-metric multidimensional scaling (NMDS) were used to analyze beta diversity between the different groups to compare the intestinal microbial composition between wild and captive R. brelichi . Based on the Bray–Curtis similarity distance algorithm, analysis of similarities was used to test the difference between the two groups (Anderson et al., 2006).
Linear discriminant analysis effect size (LEfSe) analysis was performed to identify the taxa that differed between the wild and captive groups (Segata et al., 2011). For LEfSe, Kruskal–Wallis tests were first performed among all groups of samples, followed by comparison of the selected species that differed between the two groups via Wilcoxon rank sum tests. Linear discriminant analysis (LDA) was used to sort the selected differences to obtain an LDA histogram (LDA score of > 3, p < 0.05), after which the evolutionary branching diagram was obtained by mapping the differences to the classification tree with a known hierarchical structure. Clustering ASV information was compared with the sequenced microbial genome database using PICRUSt2 software to obtain the functional types and abundance of the corresponding species in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Douglas et al., 2019) (https://www.kegg.jp/). Differences in KEGG pathways between wild and captive R. brelichi were analyzed using Wilcoxon rank sum tests (p < 0.05).
The rarefaction ASV curves, Venn diagrams, stacked column charts, LEfSe difference analysis diagrams, and KEGG function annotation diagrams were drawn using Parsenno Genomics Cloud Platform (https://www.genescloud.cn/home). The clustering heatmaps at the genus level, and PCoA and NMDS analysis charts were drawn using OmicStudio Cloud Platform (https://www.omicstudio.cn/tool). Box plots were drawn using Origin software.