Five white-faced flying squirrels (Petaurista alborufus lena) were captured from the mountains of Taiwan. A wild animal collecting permit (No. 0990007029) was approved by Yushan National Park Headquarters, Taiwan. Animals were dissected immediately after death. All experiments were performed in accordance with the Wildlife Conservation Act (http://law.moj.gov.tw/Eng/LawClass/LawAll.aspx?PCode=M0120001). All five individuals were used for measurements of gut structure and feed contents, whereas three were used for the metabolomes and the remaining two were used for metagenomes/metatranscriptomes.

The length and weight (with feed contents included) of the stomach, small intestine, cecum, and large intestine (i.e., four main gut compartments) were measured and reported (Lu et al., 2012). Here, we additionally examined the pH of feed contents. Moreover, to understand how plant-based diets were processed along the digestive tract, feed contents (~10 g from each gut compartment) were passed through the graded sieves (pore sizes: 4, 2, 1, 0.5, 0.25, and 0.125 mm) by wet-sieving (Clauss et al., 2002). Thereafter, feed particles of each size fraction (seven groups: > 4.0, 2–4, 1–2, 0.5–1, 0.25–0.5, 0.25–0.125, and <0.125 mm) were transferred onto a Petri dish, dried at 60°C for 24 h, and weighed after cooling to room temperature. Mean particle size of feed contents from each gut compartment was calculated by fitting a normal distribution (Fritz et al., 2009). Based on these measurements, we sketched the shape of the flying squirrel's digestive tract and marked physiologic characteristics of each gut compartment (Figure 1). In addition, for a visual image of gut microbes, feed contents from the cecum were fixed and observed with a Hitachi S-4700 field emission scanning electron microscope (FE-SEM). Specifically, specimens were pre-fixed with 2.5% glutaraldehyde in 0.1 M cacodylate buffer at 4°C, post-fixed with 1% osmium tetroxide in 0.1 M cacodylate buffer at room temperature, then dehydrated through a graded series of ethanol until they were in 100% ethanol, and then they were put into 100% acetone. Specimens were critical point dried, sputter coated with gold, and examined with a FE-SEM at 15 kV.

For constructing metabolomes, feed contents from three individuals were frozen (−20°C), immediately after death. For each individual, a total of 15 samples from the stomach, small intestine, cecum, and large intestine (2, 5, 4, and 4, respectively; each ~0.5 g) were collected (Supplementary Figure 1). Each sample was homogenized with 3-fold distilled water for 5 min. The supernatant containing metabolites was retrieved by centrifugation (10,000 × g for 10 min) and mixed with 100% methanol (Sigma-Aldrich, CHROMASOLV®, for HPLC) at a ratio of 1:3, and repeated once. The supernatant was lyophilized and dissolved in 100 μl distilled water for metabolomic detection by liquid chromatography (LC)—electrospray ionization (ESI)—mass spectrometry (MS), composed of an ultra-performance LC (Ultimate 3000 RSLC, Dionex) coupled with an ESI source of quadrupled time-of-flight MS (maXis HUR-QToF system, Bruker Daltonics). Metabolites were separated by reversed-phase liquid chromatography on an HSS T3 C18 column (2.1 × 100 mm; Walters). The LC parameters were autosampler temperature at 4°C and injection volume of 10 μl with flow rate at 0.4 ml/min. Elution started from 99% mobile phase A (0.1% formic acid in pure water) and 1% mobile phase B (0.1% formic acid in ACN). Mobile phase B was held at 1% for 0.5 min, raised to 60% in 6 min, raised to 90% in 0.5 min and held for 1.5 min, and lowered to 1% in 5 min. Finally, the column was equilibrated by pumping 99% mobile phase A for 4 min. The LC–ESI–MS chromatograms were acquired by the following parameters: 190°C dry gas at 8 L/min flow rate, and 1.4 bar nebulizer gas with 4,500/3,500 capillary voltage for positive/negative ion modes. The m/z values in mass spectra were recorded for further data processing.

Metabolomic data acquired from the LC–ESI–MS chromatograms were processed by TargetAnalysis (Version 1.1, Bruker Daltonics) and XCMS (Smith et al., 2006) with optimized parameters for Bruker Q-TOF mass spectrometer (Tautenhahn et al., 2008). Metabolites were identified by MetaboSearch (Zhou et al., 2012) and matched with theoretical m/z values against the Madison Metabolomics Consortium Database (Cui et al., 2008). Identified metabolites were further annotated according to the reference library of KEGG compounds (Hattori et al., 2003), with the tolerance of LC peaks within 0.3 min and signal intensities >1,000 counts. These compound annotations were used for profiling chemical environments of each gut compartment.

The overall similarity/dissimilarity among metabolomic profiles was assessed with Jaccard distance and displayed in an ordination diagram using principal coordinates analysis (PCoA), generated by the QIIME (Version 1.9) pipeline (Caporaso et al., 2010). Detected KEGG compounds were further grouped according to their chemical structures and cellular functions (Hattori et al., 2003), with an emphasis on the quantity and signal intensity of phytochemicals (leaf-derived compounds) across gut compartments. Multiple comparisons were done by ANOVA coupling with Scheffé's test, using “stats” and “agricolae” packages (Mendiburu, 2016) in R (R Development Core Team, 2016).

Finally, metabolomic data were compiled with metagenomic and metatranscriptomic data (described below) to detect potentially crucial functional reactions of cecal microbiota under the leaf-based diet.

For constructing metagenomes and metatranscriptomes, feed contents from two individuals were collected immediately after death and placed in RNAlater solution (Ambion, Life Technologies). For each individual, samples from the cecum (Supplementary Figure 1, sites 8–11, each ~5 g) were pooled to represent the entire cecal microbiota. Total DNA and RNA were isolated using the AllPrep DNA/RNA Mini Kit (QIAGEN), according to manufacturer's instructions. Briefly, cecal contents were centrifuged (10,000 × g for 10 min) to remove RNAlater and re-suspended in lysis buffer (TE buffer with 5 mg/mL lysozyme) at room temperature for 5 min. After addition of RLT buffer (with beta-mercaptoethanol), lysate was homogenized by passing it 10 times through a 20-G needle with a 1-ml syringe. The DNA molecules were purified through an AllPrep DNA spin column, and RNA molecules (>200 bp, without small rRNAs and tRNAs) were isolated with an RNeasy spin column with DNase treatment (QIAGEN). Furthermore, 16S rRNAs and 23S rRNAs (typically accounting for >80% of total RNA molecules) were removed using a MICROBExpress™ Bacterial mRNA Enrichment Kit (Ambion, Life Technologies), according to the manufacturer's instructions. Quantity and quality of DNA and RNA samples were estimated using a NanoDrop 2000 Spectrophotometer (Thermo Scientific) and an Agilent 2100 bioanalyzer (Agilent Technologies).

Coupled DNA and RNA samples were subjected to shotgun sequencing on a GS-FLX Titanium system (Roche Life Science) for characterizing metagenomes and metatranscriptomes. For RNA samples, double-stranded cDNA libraries were constructed according to the cDNA Synthesis System Kit (Roche Life Science). The DNA and cDNA libraries were converted into single-stranded DNA fragments for sequencing using a GS-FLX Titanium Rapid Library Preparation Kit (Roche Life Science). Two full runs of shotgun sequencing were conducted. Sequence data were submitted to the NCBI Sequence Read Archive (SRA; BioProject accession PRJNA267179).

Raw reads of metagenome (DNA-based) and metatranscriptome (RNA-based) sequences were filtered to remove exact duplicates (possible artifacts from emulsion PCR) using CD-hits (Li and Godzik, 2006), and low-quality parts (Phred quality < 25, N content > 3%, sequences < 100 bp) were trimmed using SeqClean (http://seqclean.sourceforge.net). After quality control, a total of 569,349 metagenome reads and 483,241 metatranscriptome reads (both with average length > 300 bp) were used for bioinformatic analyses. Potential rRNA sequences were identified using BLASTn (cut-off e-value < 1e-5, with > 90% alignment) against the SILVA (containing eukaryotic and prokaryotic ribosomal RNA sequences) database (Pruesse et al., 2007). Remaining non-rRNA sequences were considered as mRNA reads for protein-coding gene annotation. Detailed sequence statistical information is shown in Supplementary Table 1.

To provide taxonomic profiles of metagenomes and metatranscriptomes, putative mRNA reads were mapped to NCBI-nr (non-redundant protein database of National Center for Biotechnology Information) database (Sayers et al., 2012) using BLASTX (cut-off e-value < 1e-5) to identify potential taxonomic origins. Taxonomic classification of sequences was retrieved and summarized according to NCBI Taxonomy. For domain—and phylum-level taxonomic profiles, all significant hits (cut-off e-value < 1e-5) were considered. For fine-level taxonomic identification, minimum identity thresholds of amino acid sequences were further restricted at 45 and 65% for family- and genus-level taxonomic groups, respectively (Konstantinidis et al., 2017).

To provide functional profiles of metagenomes and metatranscriptomes, FragGeneScan (Rho et al., 2010) based on the hidden Markov model was used to predict open reading frames (ORFs) on putative mRNA reads. The ORFs were queried to identify conserved protein families and domains with amino acid sequences against: (1) COG (Clusters of Orthologous Groups of proteins) database (Galperin et al., 2015); (2) Pfam (Protein families) database (Finn et al., 2016); and (3) KEGG (Kyoto Encyclopedia of Genes and Genomes) database (Kanehisa and Goto, 2000). Specifically, COG annotation was performed using a reverse position-specific BLAST algorithm (RPS-BLAST, cut-off e-value < 1e-5) against NCBI COG database (latest update 2017/3/28). Pfam annotation was performed using a HMMER3 (cut-off e-value < 1e-5) against Pfam database (version 31.0). In addition, to focus on carbohydrate-degrading enzymes, glycoside hydrolases (GHs) were determined according to correspondence between Pfam and CAZy (Carbohydrate-Active enZYmes; Cantarel et al., 2009). KEGG Orthology (KO) annotation was performed using KEGG GhostKOALA (Kanehisa et al., 2016), which is designed for metagenome sequences. Moreover, to reconstruct potentially crucial functional reactions conducted by the cecal microbiota, compounds, and genes detected in the metabolome, metagenome, or metatranscriptome were all mapped to KEGG modules (Takami et al., 2012) and pathways (Hattori et al., 2003).

To evaluate individual variation, Pearson correlation coefficients (r) were used to assess degree of similarity between results derived from two flying squirrel individuals, using the “stats” package of R (R Development Core Team, 2016).

To provide an overview of functional characteristics of mammalian gut microbiota, in addition to 2 metagenome datasets from the flying squirrel's cecum, 4 metagenome datasets from cow's rumen (Brulc et al., 2009), and 39 metagenome datasets from fecal samples of zoo mammals (Muegge et al., 2011) were downloaded and analyzed. All datasets were processed with the same procedures for the ORF prediction and functional annotation, as mentioned above. Overall similarity/dissimilarity of mammalian gut metagenomes was assessed with Bray-Curtis dissimilarity and displayed in an ordination diagram using nonmetric multidimensional scaling (NMDS). An ANOSIM was done to determine whether there were significant differences in functional profiles (based on COG or Pfam) among five metagenome groups (i.e., flying squirrel's cecum, cow's rumen, zoo carnivores' feces, zoo omnivores' feces, and zoo herbivores' feces). Moreover, we performed a Kruskal-Wallis test with a threshold of adjusted p-value (FDR) <0.05 to assess significant gene families and detect different relative abundances among the five metagenome groups. Then, Dunn's test was used as a post-hoc procedure for pairwise multiple comparisons to detect the uniqueness of the flying squirrel's cecal metagenome, compared to all other gut metagenomes. Those statistical assessments were conducted using “vegan” (Oksanen et al., 2017), “RVAideMemoire” (Hervé, 2017), and “stats” packages in R (R Development Core Team, 2016).

The leaf-eating flying squirrel (Figure 1A) had a skull and molars (Figures 1B,C) that were specialized, facilitating mastication of tree leaves. In addition, shape and size both varied substantially across four main gut compartments, with the cecum being the largest chamber (Figure 1D). The pH of feed contents also differed among gut chambers (Figure 1D), with stomach contents having a relatively low pH (4-5), small intestine contents being slightly alkaline (pH 7-8), and contents in the cecum and large intestine being mildly acidic (both pH 6-7). Regarding feed particle sizes (Figure 1D and Supplementary Figure 2), feed particles in the stomach were already fine (on average, ~1.5 mm). Nevertheless, compared to dry and coarse fibers present in the stomach, feed contents in the cecum were even finer (on average, ~0.6 mm) with a sludge-like texture (Supplementary Figure 3). Moreover, in microscopic evaluation of cecal contents, various microbial cells were closely attached to plant debris (Figures 1E–G).

Based on LC–ESI–MS metabolomic profiles, each gut compartment contained its own unique compound composition in feed contents, with clear differences between samples from the upper (stomach and small intestine) vs. lower (cecum and large intestine) portions of the digestive tract (Figure 2). Overall, there were ~600 phytochemicals (mainly leaf-derived metabolites) detected among all gut samples. Several phytochemicals (~20%) were present only in the stomach and small intestine, with relatively fewer phytochemicals detected in the cecum and large intestine (Figure 3), especially flavonoids, phenylpropanoids, and polyketides (Figures 3D–F). Moreover, considering the relative intensity of those detected phytochemicals, their signal intensities were also much lower in the hindgut than in the foregut (Supplementary Figure 4).

Based on the annotation of protein-coding sequences, at both DNA and RNA levels, cecal microbiota were dominated by bacteria (average, 97.94% of DNA reads and 88.24% of RNA reads; Supplementary Table 2), followed by eukaryotes (0.77 and 11.37%), archaea (0.46 and 0.11%), and viruses (0.09 and 0.03%). Targeting bacterial composition, the phylum Firmicutes was extremely abundant in both metagenomes (~90%) and metatranscriptomes (~85%), with Actinobacteria, Proteobacteria, Bacteroidetes constituting only a small fraction (Supplementary Table 3). At a genus level, 577 bacterial taxa were detected, whereas only 30 taxa contained >0.5% reads in any library (Figure 4). Among these abundant genera, Clostridium, Ruminococcus, and Eubacterium (all belonging to Firmicutes) were the top three in all libraries, contributing to the majority (up to 50%) of annotated protein-coding genes in the cecal metagenome and metatranscriptome.

According to the COG annotation (Figure 5 and Supplementary Figure 5), functional gene families had distinct abundance patterns at DNA and RNA levels. Specifically, the most abundant COG in the cecal metagenome was assigned to the ABC-type multidrug exporter component (COG1132; involved in defense mechanisms; representing ~2% of hits; Figure 5). By contrast, the most abundant COG in the cecal metatranscriptome was assigned to the flagellin protein (COG1344; involved in cell motility), followed by several ABC-type sugar importer components (e.g., COG3839; involved in carbohydrate transport and metabolism). In addition, ABC-type importer components for short peptides (e.g., COG0747, involved in amino acid transport and metabolism) also had relatively high abundances at the RNA level.

Contrasting functional gene profiles in the metatranscriptome vs. metagenome were also revealed by KEGG annotation (Figure 6 and Supplementary Figure 6), with distinct types of membrane transport components dominant at RNA—vs. DNA-levels, similar to COG results (Figures 5, 6). In combination with metabolomic data, we identified KEGG pathways with high coverage of KOs and relevant compounds, including pathways involved in substrate-induced cell motility (i.e., bacterial chemotaxis and flagellar assembly; Supplementary Figures 7, 8) and pathways for biosynthesis of cellular components from plant-sourced nutrients (e.g., cellobiose, xylose, and arabinose were transported into bacterial cells and converted into other macromolecules, such as amino acids, nucleotides, and peptidoglycans; Supplementary Figures 9–16).

Among those abundantly detected COGs and KOs (Figures 5, 6), it was notable that apart from ABC-type transport systems and genes associated with “housekeeping functions” (such as DNA gyrases, RNA polymerases, and translation elongation factors), one particular carbohydrate-degrading enzyme, beta-glucosidase (COG1472 or K05349), had high abundance in both metagenomes (~0.5%) and metatranscriptomes (~1.0%). To better characterize diversity of carbohydrate-degrading enzymes, we focused on glycoside hydrolase (GH) groups, based on the Pfam annotation. A total of 60 GH groups were detected at the DNA level, of which, 39 GH groups were also detected at the RNA level (Supplementary Table 4). According to enzymatic activities of those GH groups, diverse carbohydrate-degrading genes specific for beta-glycosidic linkages in plant polysaccharides / oligosaccharides were detected in the cecal metagenome and metatranscriptome (Figure 7). Among them, GH3 was the most abundant group at both DNA and RNA levels. Moreover, compared to the DNA-level background, GH9, GH48, GH43, and GH53 were enriched at the RNA level (Figure 7).

Basically, taxonomic and functional profiles from individuals (FS1 and FS2) were comparable (Figures 4–7), with a high similarity (r > 0.9 and P < 0.01) for either metagenomes or metatranscriptomes.

To provide an overview of functional characteristics of mammalian gut microbiota, in addition to 2 metagenomes from the flying squirrel's cecum, we conducted a meta-analysis on mammalian gut metagenomes with published datasets (including four datasets from the cow's rumen, and 39 datasets from fecal samples of zoo carnivores/omnivores/herbivores). Functional gene compositions based on the COG annotation (Figure 8) demonstrated that: (1) functional profiles of cecal metagenomes and ruminal metagenomes were distinct from those based on fecal samples; and (2) functional profiles of cecal metagenomes were also distinct from those of ruminal metagenomes. Similar functional relationships were obtained based on the Pfam annotation (Supplementary Figure 17).

To reveal unique functional characteristics of the flying squirrel's cecal microbiota, we conducted statistical tests among 5 metagenome groups (i.e., flying squirrel's cecum, cow's rumen, zoo carnivores' feces, zoo omnivores' feces, and zoo herbivores' feces). In comparisons based on relative abundances of COG functional categories, many categories had significant differences among 5 metagenome groups (Supplementary Table 5). Moreover, post-hoc pairwise comparisons focusing on specific differences between the flying squirrel's cecum and the other four groups (Supplementary Table 6) indicated that, despite several significant pairs, the flying squirrel's cecum had clear differences from all other groups in one particular category, namely defense mechanisms (V). Specifically, compared to other mammalian gut metagenomes, the flying squirrel's cecal metagenome contained a significantly higher percentage of gene function involved in defense mechanisms (Supplementary Table 5), especially for those COG gene families assigned to multidrug efflux pumps (Supplementary Table 7).

In addition, in a comparison of glycoside hydrolase (GH) groups among all 5 metagenome groups, the flying squirrel's cecal metagenome contained relatively high levels of GH3, GH43, and GH53 (Supplementary Table 8); those GH groups were specific for degradation of oligosaccharides and were actually enriched in cecal metatranscriptions (Figure 7).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.