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.