Mapping gastrointestinal gene expression patterns in wild primates and humans via fecal RNA-seq

Background Limited accessibility to intestinal epithelial tissue in wild animals and humans makes it challenging to study patterns of intestinal gene regulation, and hence to monitor physiological status and health in field conditions. To explore solutions to this limitation, we have used a noninvasive approach via fecal RNA-seq, for the quantification of gene expression markers in gastrointestinal cells of free-range primates and a forager human population. Thus, a combination of poly(A) mRNA enrichment and rRNA depletion methods was used in tandem with RNA-seq to quantify and compare gastrointestinal gene expression patterns in fecal samples of wild Gorilla gorilla gorilla (n = 9) and BaAka hunter-gatherers (n = 10) from The Dzanga Sangha Protected Areas, Central African Republic. Results Although only a small fraction (< 4.9%) of intestinal mRNA signals was recovered, the data was sufficient to detect significant functional differences between gorillas and humans, at the gene and pathway levels. These intestinal gene expression differences were specifically associated with metabolic and immune functions. Additionally, non-host RNA-seq reads were used to gain preliminary insights on the subjects’ dietary habits, intestinal microbiomes, and infection prevalence, via identification of fungi, nematode, arthropod and plant RNA. Conclusions Overall, the results suggest that fecal RNA-seq, targeting gastrointestinal epithelial cells can be used to evaluate primate intestinal physiology and gut gene regulation, in samples obtained in challenging conditions in situ. The approach used herein may be useful to obtain information on primate intestinal health, while revealing preliminary insights into foraging ecology, microbiome, and diet. Electronic supplementary material The online version of this article (10.1186/s12864-019-5813-z) contains supplementary material, which is available to authorized users.


Background
Difficulties in obtaining nucleic acids from representative samples make it necessary to implement non-invasive approaches to explore genetic traits in wild animals and humans [1,2]. In wild apes, non-invasive genotyping methods from feces have been invaluable in determining long term variation in social group dynamics, estimating population abundance, and performing individual tracking [3][4][5][6][7]. Other applications of feces-based, noninvasive genetic methods in wild non-human primates include inferences on evolutionary history [8], seed dispersal patterns [9], and assessment of diet [10][11][12]; all facilitated by advances in high throughput genomic methods [13].
The application of other noninvasive genetic methods from fecal samples, such as transcriptome analyses, more directly reflective of intestinal immune or metabolic phenotype at the moment of sample collection, are more challenging [14]. Apart from difficulties in obtaining representative tissue samples, significant challenges associated with these analyses include preservation and integrity of samples and nucleic acids in field conditions, contamination with other biological materials in feces (e.g. microbes, dietary materials) and hence, inability to detect nucleic acid signals reflecting animal physiological status. Nonetheless, transcriptomic (mRNA) analyses from fecal samples in humans and animals, targeting intestinal epithelial cells (IECs) via RT-qPCR or RNA-seq, have been previously used to monitor inflammatory biomarkers during diarrhea [15], to profile intestinal metabolism and immunity after birth [16,17], to assess colon cancer risk [18], and to monitor mucosal transcriptomics in response to non-steroidal anti-inflammatory drug enteropathy [19]. In addition, other techniques, such as droplet digital PCR, have been useful in detecting transcriptional markers of inflammation from human fecal samples [20].
Thus, transcriptomic analyses from fecal samples can be a promising, non-invasive approach to assess intestinal function, when obtaining tissue samples is a challenging task. However, this approach is much less common in wild or captive animals, or human populations living in challenging field conditions, where this noninvasive analysis tool could be invaluable in monitoring for physiological and health status in the context of host health and conservation. Therefore, this study focuses on a preliminary assessment of this noninvasive approach to profile intestinal function in wild primates and human populations in situ.
To that end, we validate RNA-seq analyses from fecal samples, with the goal of quantifying differential intestinal gene expression patterns between wild western lowland gorillas (Gorilla gorilla gorilla) and a human forager population; the BaAka hunter-gatherers, from the Dzanga Sangha Protected Areas, Central African Republic. We chose these two closely related, but different primate species, expecting to detect wide distinctions in intestinal gene expression, while validating the potential usefulness of the technique for monitoring health and physiological status in habituated or captive wild primates, and isolated human populations when repeated sampling is possible.
The RNA-seq pipeline followed herein sought to preserve the integrity of intestinal epithelial cells in field conditions, while maximizing the amount of transcripts (mRNA reads) obtained by comparing different eukaryotic mRNA enrichment methods (poly(A) mRNA enrichment and/or rRNA depletion). Our data showed that although the proportion of gorilla and human mRNA signals recovered from fecal samples via RNAseq is significantly small, differential immune and metabolic gene expression patterns denoting intestinal functional distinctions between the two closely related primate species can still be detected.

Results
Finding recovery of human and gorilla mRNA reads from fecal samples using poly(a) mRNA enrichment plus rRNA depletion vs rRNA depletion alone To evaluate two methods for eukaryotic mRNA enrichment from fecal samples, we conducted RNA-seq in 4 sample sets (1 from gorilla and 3 from humans) using poly(A) mRNA enrichment plus rRNA depletion and rRNA depletion alone, generating between 4.5-16.7 M paired end reads (Table 1). After quality control and alignment of pair-end reads to the latest genome drafts for G. g. gorilla (gorGor4.1/gorGor4) and H. sapiens (GRCh38/hg38), the percentage of reads aligning to the respective host genomes ranged from 0. 44 to 1.41% for samples undergoing poly(A) mRNA enrichment plus rRNA depletion. In contrast, samples in which only rRNA depletion alone was used showed from 0.04 to 0.27% alignment rate to the target genomes. Thus, the combined method showed from 11 to 15 times more mappability to the host genomes (paired t-test, P = 0.02), while effectively removing rRNA (0.0 to 4.4%) (Additional file 1: Table S1). Non-host RNA-seq reads in the fecal samples primarily map to the microeukaryotes As poly(A) enrichment plus rRNA depletion yielded the most mappability to human and gorilla genomes, we used this approach for RNA-seq analyses on n = 10 fecal samples of BaAka foragers and n = 9 samples of wild gorillas, with sequencing depth varying from 2.2 million to 21. 52 million in all samples (Additional file 1: Table S1). However, before proceeding with differential gene expression analyses, we first determined the identity of reads not mapping to the target genomes. To that end, 10,000 randomly selected unmapped reads from each sample were aligned against the NCBI NR database. A strict identity cut-off was used to achieve maximum specificity and hence a large proportion of non-host reads remained unclassified (78.08-98%) (Additional file 2: Table S2). Out of all classified reads, most were of eukaryotic origin (~74.77 ± 26.16%) followed by reads mapping to bacteria (~15.13 ± 15.06%) (Fig. 1a, Additional file 1: Table S1). Upon closer inspection to the eukaryotic fraction, the majority of reads mapped to microeukaryotes (~53.4 ± 23.9% of all eukaryotic reads, excluding gorillas and humans), followed by reads mapping to metazoans (~11.6 ± 9.33%), plants (~6.4 ± 6.5%) and fungi (~3.2 ± 3.5%) (Fig. 1b, Additional file 1: Table S1). Of note, gorillas harbored significantly higher number of reads mapping to metazoans (q = 2.17e-05, Wilcoxon signed-rank test) and plants (q = 0.001, Wilcoxon signed-rank test), and tended to have higher percentage of reads mapping to fungi (q = 0.09, Wilcoxon test). Humans tended to show higher percentage of reads mapping to bacteria compared to gorilla (q = 0.07, Wilcoxon signed-rank test) (Fig. 1c). Further inspection of the metazoan-derived reads indicated similar proportion of reads mapping to arthropods (q = 0.68, Wilcoxon signedrank test) and nematodes (q = 0. 51, Wilcoxon signed-rank test) in both gorillas and humans (Fig. 1d).

Differential expression profiles in IECs of gorillas and humans
The percentage of reads mapping to human genomes in BaAka fecal samples ranged from 0.40 to 4.99% (1.37 ± 1.36%); while in gorilla fecal samples, the mapping rate to gorilla genomes ranged from 0.26 to 3.1% (1.37 ± 1.00%). Thus, no differences in fecal mappability between gorillas and human were detected (t-test, p = 0.99). From the alignments, transcripts per million (TPM) values for 45, 717 gorilla transcripts and 179,974 human transcripts were obtained and then a non-redundant (NR) set of 15, 602 common genes, based on gene symbol, and known function in both species, was identified. This NR gene set was used for all downstream gene expression comparisons between humans and gorillas using the DESeq2 R package [21,22]. A principal coordinate analysis based on Bray-Curtis distances showed significant differential gene expression patterns in the gastrointestinal tract between humans and gorillas (PERMANOVA, R 2 = 0.33, p < 0.001) ( Fig. 2a), which was further supported by a hierarchical clustering analysis (Fig. 2b). Furthermore, to validate these differential gene expression patterns, we downloaded RNA-seq data targeting IECs shed in fecal samples of human infants. This analysis showed that intestinal gene expression profiles in fecal samples of human infants are closer to those seen in the BaAka hunter-gatherers, based on Bray-Cutis distances (Principal coordinate analysis, PERMANOVA, R 2 = 0.43, p < 0.001, Wilcoxon rank sum test, p < 0.001, Fig. 2c-d).
Fold change analyses on the total NR gene set supported by false discovery rate-adjusted q-values (Benjamini & Hochberg correction) identified 212 genes showing significant differential expression between gorillas and the BaAka foragers. One hundred eighty four genes were up-regulated in gorillas while only 28 genes were upregulated in the BaAka foragers (> 5 fold, q < 0.05, Fig. 3a). The top 50 differentially expressed genes were used to construct a heatmap showing a pool of markers involved in diverse regulatory functions ( Fig. 3b and Table 2). Additionally, 26 gene markers specific to gastrointestinal tract metabolism and immune responses were identified (Additional file 3: Table S3). The biological relevance of the differentially expressed genes detected in gorillas and humans was determined using a functional gene set enrichment analysis, as implemented in Ingenuity Pathway Analysis [23]. For this analysis, differentially expressed genes with an FDR q-value cut-off < 0.05 between the two groups were considered. This procedure identified 23 pathways (based on the z-scores) significantly regulated between gorillas and humans (Top 5 upregulated and downregulated pathways in each species can be seen in Fig. 4, all pathways in Table 3). The most upregulated pathway in gorillas was oxidative phosphorylation (ATP5F1E, ATP5MC3, ATP5MF, ATP5PD, COX5B, COX6B1, SDHB, UQCRB and UQCRQ), followed by phagocytosis in macrophages and monocytes (e.g., ACTA1, ACTB, ACTC1, ARPC3 and FCGR2A), dendritic cell maturation (e.g., FCGR2A, IL32, LEP, NFKBIA and PLCD3), interferon signaling (e.g., IFI6, IFITM1, IFITM2 and IFNGR2) and PI3K signaling in B lymphocytes (e.g., ATF3, CHP1, NFKBIA and PCLD3) ( Fig. 4 and Table 3). In contrast, humans exhibited highly up-regulated expression of genes associated with IL-8 signaling (e.g., GNB1, LASP1 and RHOB), G-protein (GNB1 and RHOB), RhoGDI (e.g., ARHGEF12, GNB1 and RHOB) and phospholipase C (e.g., ARHGEF12, GNB1 and RHOB) signaling and production of nitric oxide and reactive oxygen species in macrophages (e.g., RHOB) ( Fig. 4 and Table 3).

Discussion
Here, we attempted to maximize functional genomic information from the intestinal tract of wild lowland gorillas and the BaAka foragers from the Dzanga Sangha Protected Areas in the Central African Republic. Our main goal was to determine if the non-invasive RNA-seq approach would detect regulatory signals under selection in the gastrointestinal tract of two closely related primate species, potentially allowing us to answer critical questions on gut physiology of free-ranging gorillas and syntopically living humans. The main obstacle faced by this method was the recovery of sufficient mammalian transcripts from human and gorilla fecal samples. However, as healthy humans can shed approximately 10 billion intestinal cells per day in feces, it has been suggested that detecting mRNA signals from intestinal cells in feces is possible [24,25].
Despite success in recovering information from IECs and immune cells, percent mappability obtained in this dataset is still lower in comparison to the mappability obtained in human studies (range 0.7-37.9%, mean = 9.21 ± 13.6%) [17]. These differences in mappability may be due to distinctions in sample processing and storage; for instance, we did not perform any bowel fluid preparation to filter out fecal debris nor use a direct poly(A) RNA enrichment extraction kit as described previously [26]. Also, we could not collect fecal samples immediately after voiding, which promotes cellular death. In addition, although we only selected RNA samples of high quality (2 to 5μg of RNA, with a RIN from 7 to 8), RNA purity may have also affected the number of host-derived transcripts obtained. To address these issues, we attempted to test the combined effect of rRNA depletion plus poly(A) enrichment, which improved mammalian genome mappability. In this regard, it has been shown that poly(A) enrichment and rRNA depletion provide similar rRNA removal efficiency, coverage, mappability and transcript quantification, besides providing a less biased coverage of 3′ ends in genes [27]; however, increasing sequencing effort (~4 times sequencing depth) may be necessary when using rRNA depletion alone at the same exonic coverage [28]. This is consistent with the increased mappability reached herein, using the combined rRNA depletion and poly(A) enrichment method, as opposed to rRNA depletion alone, at the same sequencing depth.
Nonetheless, this combined approach can be considerably more costly compared to using RNA depletion alone. As such, approaches that are gene-centric and that rely on either micro-arrays or RT-qPCR to target gene markers of gastrointestinal disease [29], nutrition [30], and immune and metabolic health [31] may be more cost effective.
Comparison of human and gorilla transcriptomic data from fecal samples using expression profile of total 15,602 common or orthologous genes: a) Principal coordinate analysis based on Bray-Curtis distances, b) Hierarchical cluster analysis using hclust function based on the distances calculated using binary distance measure, c) Comparison of human-BaAka, gorilla and human-infants (downloaded) transcriptomic data using principal coordinate analysis based on Bray-Curtis distances, d) Bray-Curtis distances between these groups However, targeted approaches are unable to detect additional information on extrinsic host factors (non-host derived) such as diet and microbiome.. For instance, our results showed that from all the mRNA reads unmapped to gorilla and human genomes, those of plants, metazoa and fungi were more prevalent in gorillas than in the BaAka fecal samples. Thus, nonhost RNA-seq reads, as detected here, could potentially offer preliminary insights on dietary behaviors and microbiome. For example, it is a fact that, compared with humans, wild gorillas rely significantly more on plant-based diets [32]. Also, gorillas have been reported to have specialized gut fungal populations to process and ferment highly fibrous and complex foods, an adaptation believed to be of less importance for humans [32]. The observation that both gorillas and the BaAka humans showed similar proportions of arthropod-derived reads is more intriguing; but it may be connected with similar patterns of consumption of insects and insect-derived products by both gorillas and human foragers [33][34][35][36]. Likewise, similar prevalence of gastrointestinal nematode infections have been reported by our group in these same gorilla groups and human populations previously [37].
Nonetheless, we were unable to determine why the proportion of all nematode reads was higher in gorillas. This issue is likely a result of the limited taxonomic information that can be recovered from such a small transcript fragment, influence of food processing by the host, and issues related to the technical challenges of this endeavour. Also, due to the liable nature of RNA, DNA-based metabarcoding may be more effective to infer dietary behaviors from feces surveys, in the context of plant and insect consumption [10]. As such, all results related with non-host factors (diet and microbiomes) presented herein warrant further investigation, and comparison with targeted approaches.
Another important consideration of the proposed noninvasive approach is the kind of functional information that can be obtained from exfoliated intestinal epithelial and immune cells in fecal samples. Exfoliated cells in fecal samples may originate from any site along the gastrointestinal tract, including stomach, villus tips in the small intestine and from crypt surfaces at colon level. This information can potentially offer valuable insights on gastrointestinal homeostasis in the context of nutrition and health [25,38]. However the exact  origin of these cells is unknown, and once detached from the extracellular matrix, cells enter a stage of programmed apoptosis and anoikis. Additionally, an important concern is the degradation of RNA coming from small intestinal cells in the lower GI tract and in fecal samples; these limitations may have caused us to miss important expression patterns. Regardless, the results presented here, along with validation of our data by contrasting fecal mRNA reads of with a previously published human dataset, demonstrate the biological value of fecal RNA-seq to mine for differential gene expression patterns in gut tissues [17,19]. For instance, we report a higher number of upregulated pathways in gorillas (18 pathways) than in humans (5 pathways), with oxidative phosphorylation in gorillas being the most upregulated pathway in the dataset. Colonic expression of genes involved in energy metabolism, such as oxidative phosphorylation, may be enhanced via the colonic microbiome and its metabolites, specifically butyrate [39]. Thus, it is expected that diets that promote increased fermentation selectively favor this regulatory trait [40]. This observation explains the dramatic downregulation of oxidative phosphorylation in the BaAka humans, who incorporate significantly less plant material in their diets. Other gorilla-specific pathways mainly denote maintenance of intestinal barrier integrity, cell architecture and immune functions (Table 3). [A.G2] For example, the upregulation of adaptive immune pathways in gorillas such as phagocytosis in macrophages, dendritic cell maturation and interferon signaling may reflect capacity to recognize and eliminate potential external insults and maintain a balance between tolerance and inflammation in the intestinal environment. Regulatory adaptations to counteract potential pathobionts in the gorilla gut, while tolerating commensals, may be reinforced by reliance on diets that could also inhibit pathogens and stimulate beneficial microbes, such as those rich in phenolics and fermentable dietary substrates (fiber) [41,42]. These observations may explain why gorillas, compared with the BaAka humans, exhibit down regulation of pathways involved in proinflammation.
For instance, in the BaAka, most upregulated pathways (4 out of the 5 detected) denote increased inflammatory responses. The more upregulated pathway in the BaAka was IL-8 signaling, which mediates recruitment of proinflammatory cytokines (CXCL8) to sites of infection [43,44]. In this regard, it has been reported that African foragers exhibit evolutionary adaptations to counteract increased susceptibility to intestinal infection, by enhancing immune responses to pathogens [37,[45][46][47]. Other upregulated pathways seen in the BaAka that may support increased responses to pathogens and inflammation are Gɑq signaling and production of nitric oxide (NO) and reactive oxygen (RO) by macrophages. Gɑq signaling plays a role in Paneth cell maturation, intestinal barrier integrity and protection against luminal pathogens by secreting alpha-defensins [48], and NO and RO production in macrophages is correlated with levels of pathogenic infections in monocytes and inflammation [49,50]. Moreover, phospholipase C signaling is also involved in protection from a hostile luminal environment and tissue restitution after intestinal epithelial cell injury, in conjunction with Gɑq [51,52]. Based on these results, and the differential expression of pro-inflammatory pathways detected, it is tempting to speculate that humans, in comparison with nonhuman primates are evolutionarily primed for inflammatory responses in response to specific diets, environments or the microbiome; nonetheless, these hypotheses warrants further investigation.

Conclusions
In summary, we advance a noninvasive approach that evaluates primate intestinal physiology, and that can be used in situ, by measuring gene expression patterns in fecal samples. Although the proportions of IEC mRNA patterns obtained in these conditions are not fully representative, compared to standard RNA-seq approaches, we show that this method has the potential to reveal critical gene expression information when comparing primates exposed to different physiological and environmental conditions. In addition to providing insights into differences among two primate species in intestinal cell function, this approach could reveal fine-grained distinctions among conspecific populations, including differences by age or sex. Broadly, these results have important implications for understanding the immune and metabolic gut environment of humans and animals in filed conditions, when invasive samples are unattainable. In the future, applying this procedure in tandem with methods that provide a more reliable and targeted picture of diet and microbiome (e.g. DNA barcoding, metabolomics and metagenomics) may offer a powerful tool to characterize diet-host-microbe interactions in the context of primate nutrition, health and conservation.

Study site and sampled objects
The study was conducted in Dzanga-Sangha Protected Areas (DSPA) in the Central African Republic. Activities in DSPA are directed by the DSPA administration under the collaborative management of the CAR Government and World Wildlife Fund. The climate is characterized by marked seasonal variation [53], with dry (November-March) and wet (April-October) months [54]. Human population density is low, estimated at one person per square kilometer, with the greatest concentration (60% of people) located in the village of Bayanga [55,56]. In May 2016, we sampled three habituated gorilla groups (Makumba, Mayeli, Mata) around two permanent

Sample collection and preparation
Fecal samples (1 g) of lowland gorillas and humans were collected within 5 min after defecation and divided in three subsamples (to prevent repeated thawing-freezing if samples need to be analyzed more times) which were immediately in the field preserved in RNAlater (Qiagen, Hilden, Germany) (to prevent repeated thawing-freezing if samples need to be analyzed more times) and left overnight at room temperature. Next day the samples were stored in the mobile freezer powered by solar battery. The necessity to collect the sample immediately after defecation restricts this method only to habituated apes as samples from the unhabituated ones are usually several hours old. In the field, the samples were stored in a mobile freezer powered by solar battery. During the transport to the laboratory in the Czech Republic samples were kept at − 20°C in a mobile freezer powered by a car battery, while during the international air transport, samples were transported in dry ice.
After the arrival to the laboratory in the Czech Republic, samples were stored at − 80°C and RNA isolated within a month.

RNA isolation
RNAlater-preserved and frozen feces were thawed on ice to allow collection of approximately 250 mg of feces. The RNA was isolated using the RNeasy Midi Kit (Qiagen, Hilden, Germany) following manufacturer's instructions. For the cell lysis, β-mercaptoethanol was added to the RLT buffer and RNA was eluted in two steps to two separate tubes using 150 μL of RNAse-free water each time and the isolated RNA was immediately stored at − 80°C.

Next-generation sequencing
Total RNA was treated with DNAse and followed by RiboZero rRNA removal kit (bacteria, MRZB12424 Illumina). Samples that yielded from 2 to 5μg of RNA, with a RIN from 7 to 8 were used for downstream applications. RNA was QC'd on the Bioanalyzer using Agilent RNA 6000 Nano Kit (5067-1511), with accepted values of remaining rRNA at 20%. The RNA samples were subjected to mRNA isolation in an effort to capture the host mRNA. The entire sample was used for mRNA isolation using NEBNext poly(A) mRNA Magnetic Isolation Module (NEB E7490L) on the Apollo 324 system. The resulting mRNA was used for generating Illumina libraries using the PrepX RNA-seq for Illumina Library Preparation Kit, 8 sample (Wafergen cat#400039) on the Apollo 324 system with the scriptSeq Index PCR primers (epicenter SSIP1234). Final libraries were side selected at 0.9X. The libraries were QC's for proper sizing and the absence of adapter dimers on the Bioanalyzer using the Agilent DNA high sensitivity kit (5067-4626) followed by qPCR using the library quantification kit (Kapa Biosystems KK4835) on a 7900HT Fast real-time PCR system (Thermo-Fisher). The libraries were sequenced on the NextSEq 500 at 2*75 base pair read length.

Gene expression quantification
Quality reports for raw reads were generated using FastQC toolkit. Trimmomatic was used for the trimming of raw reads using "ILLUMINACLIP:$TRIMMOMATIC/ adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:36". Trimmed reads obtained from the gorilla samples were mapped against the gorilla genome and from the human samples were mapped against the Human genome using Kallisto [57]. Aligned counts (TPM values) for each gene from each sample were used for the subsequent analyses. Variance for each was calculated using DESeq2 package after performing the variance stabilizing transformation. Principal coordinate and hierarchical clustering analysis were performed in R. Fold changes and adjusted p-values (q-value) for Gorilla vs Human comparison were calculated using DESeq2 at alpha = 0.05 and pAdjustMethod = "BH". Top 50 differentially expressed genes (q-value < 0.05 & log2 fold change > 10) and having maximum variance were selected for the association analysis using cluster heatmap (distfun = "euclidean", hclustfun = "complete").

Ingenuity pathway analysis
Ingenuity Pathway Analysis (IPA) version 2.0 was used to perform the functional enrichment analysis (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/). All differentially expressed genes (q-value < 0.05) were uploaded along with their log2 fold change information. Out of total 212 genes, 184 genes significantly enriched in gorilla (log2 fold change > 5) and 28 genes significantly enriched in humans (log2 fold change > 5) were used for this analysis. The z-score was used to infer the activation state which was calculated using the log2 fold change values of each gene.

Identification of different taxonomic groups
To identify the taxonomic fractions other than the host genome, 10,000 randomly selected unmapped reads from