A clinically validated human saliva metatranscriptomic test for global systems biology studies

Chronic diseases are a leading cause of morbidity and mortality globally, yet treatment and prevention options have shown limited success. Importantly, chronic diseases have been shown to be weakly associated with genetics; instead, the human microbiome and alterations to human and microbial gene expression patterns have been identified as the underlying driver of many chronic diseases [1,2]. To identify the etiology of chronic diseases and develop more effective preventative measures, comprehensive gene expression analysis of the human body and associated microbiomes is needed. Saliva is noninvasive, is easy to obtain and contains a diverse microbiome involved in many fundamental aspects of human physiology [3–6]. Saliva is also a suitable clinical specimen for the identification of human pathogens, such as SARS-CoV-2 [7]. Human transcripts in saliva can also produce an informative snapshot of human gene expression and its possible roles in human health and disease [8–10]. Saliva presents an opportunity to investigate the oral microbiome and oral human transcriptome and their roles in disease pathogenesis [6,11].

There is a critical need for noninvasive, at-home sample collection that can be used to predict, diagnose and inform therapeutic options for chronic diseases. Critically, methods for investigating multiple aspects of human physiology need to be pioneered to adequately shed light on the diversity of factors impacting chronic diseases. For example, systems biology transcriptomic methods have already been employed with great success in cancer therapies [12–14]. Methods to analyze the transcriptome are becoming invaluable in the understanding of chronic diseases with unknown etiologies, but large-scale adoption of metatranscriptomic analysis has not yet been possible due to the lack of low-cost and scalable methods.

The saliva microbiome has been associated with numerous chronic diseases, including Alzheimer's disease [15], Parkinson's disease [16], autism spectrum disorder [8], cancers [17], cardiovascular disease [18], diabetes [19], obesity [20] and autoimmune disorders [21,22]. As an example of the direct effect of salivary microorganisms on disease, Fusobacterium nucleatum originating in the oral cavity has been shown to contribute to colorectal cancer development [23,24]. The oral microbiome has also been shown to have a strong influence on the gut microbiome and the immune system, which are involved in a variety of chronic diseases [25]. In addition to the clear role of the saliva microbiome in human health and disease, there is substantial information to be gained from human gene expression patterns in saliva. Alterations in salivary human gene expression and epigenetic markers have been observed in a range of disorders such as autism spectrum disorder [9], Parkinson's disease [26] and traumatic brain injury [27]. Saliva is clearly an informative sample type in human health and disease, with predictive biomarkers of several chronic diseases.

Despite the abundance of literature showing clear connections among the oral microbiome, human gene expression and chronic diseases, most of the evidence is based on 16S rRNA gene sequencing or metagenomic sequencing methods. Both 16S and metagenomics have intrinsic limitations that minimize the ability to discover biomarkers and to provide actionable therapeutic interventions. 16S is limited in its ability to reliably identify microorganisms at the species or strain level and, importantly, cannot provide insights into microbial or human gene expression [28–31]. Metagenomics is limited in its ability to detect RNA viruses and, again, critically cannot identify active microbial or human gene expression. The analysis of gene expression is an important component of understanding chronic diseases and has been shown to have a greater role than simply looking at the microbial composition [32,33]. These limitations of the standard methodologies result in an incomplete picture of the saliva microbiome ecosystem, which minimizes the ability to discover novel microbial or human-derived biomarkers [34]. Metatranscriptomic analysis of the saliva microbiome presents a more comprehensive and relevant snapshot of the microbial and human transcripts compared with traditional methodologies. Metatranscriptomic analyses address the limitations of both 16S rRNA gene sequencing and metagenomics by providing high-level taxonomic resolution (strain level) and by allowing for assessments of both human and microbial gene expression [35,36]. Metatranscriptomics is uniquely capable of providing much of the information needed for detailed biomarker discovery.

Saliva metatranscriptomic tests have already been used by the authors' laboratory to create comprehensive and accurate diagnostic indicators of oral cancer [37]. Importantly, these predictive models rely on microbial taxonomy, microbial gene expression and human gene expression signatures.

In addition to their utility in biomarker discovery, functional outputs are also critical in developing a comprehensive understanding of disease etiology and subsequent treatment options [38]. Gene expression analyses can also provide information about microbial pathogenicity, which is lost when simply looking at the microbiome composition. Porphyromonas gingivalis, for example, can be present in small amounts but has highly expressed functions with a strong negative impact on the microbial community and ultimately the host [32,39].

Scalable, affordable and at-home sample collection had limited metatranscriptomic clinical applicability in the past [40]. Effective RNA preservation has been challenging, traditionally requiring a cold chain that is expensive and complicated. In addition, the high abundance of noninformative RNA sequences (such as rRNA and tRNA) in clinical samples requires high sequencing depth and associated costs to overcome. By utilizing the selective depletion of noninformative RNAs, mRNA sequences can be enriched, which dramatically reduces sequencing costs while improving data resolution [35,41].

Here the authors present a comprehensive (sample collection-to-result) method for the quantitative metatranscriptomic analysis of saliva samples that can easily be applied to clinical studies and trials globally. The method is automated, high-throughput, inexpensive and clinically validated, and it includes a fully automated bioinformatic suite for strain-level taxonomic classification of all microorganisms and human genes and their quantitative gene expression levels.

Methods

Ethics statement

All procedures involving human subjects were performed in accordance with the ethical standards and approved by a federally accredited institutional review board committee. Informed consent was obtained from all participants, who were residents of the USA at the time.

Sample collection

Viome has developed a patented and convenient sample collection device that allows anyone, even children, to easily collect sufficient volumes of saliva and preserve it for metatranscriptomic analysis (Figure 1). For future studies, saliva can be collected in any device that enables mixing of the Viome RNA preservative buffer (RPB) with the sample as soon as it is collected.

Participants were instructed to fast for 8 h prior to collection, to gently rinse their mouth with water for 10 s before collecting the sample, to collect 1.2 ml of saliva within 1 h of waking up and then to dispense the RPB and shake the collection tube vigorously for 15 s.

Method validation

The performance of the saliva method was assessed by determining the method precision, sample stability and longitudinal changes in saliva. Precision was assessed in four technical replicates of saliva from ten donors. Sample stability was assessed by storing samples for 7 and 28 days at room temperature compared with the technical replicates analyzed on the day of collection (day 0). A subset of replicate samples was shipped and returned to Viome to assess the impact of shipping on sample stability. All sample stability analyses were performed with four technical replicates among three donors. Longitudinal changes and Hellinger distances in saliva were assessed by collecting saliva weekly for 5 weeks from eight donors.

To develop a comprehensive saliva cohort, saliva samples were collected from 1102 individuals. For the saliva cohort, the average age was 48 years (range: 12–88 years) with 436 males, 665 females and one participant identifying as another gender.

Laboratory analysis

The method for the analysis of saliva samples was similar to that which has been previously described by the authors' lab [42,43]. Briefly, saliva samples were chemically lysed using the RPB, followed by mechanical lysis using bead-beating beads. RNA extraction was performed with silica beads on 87.5 μl of saliva. DNA was digested using DNase. rRNAs, both prokaryotic and eukaryotic, were removed via a subtractive hybridization method. Directional, dual-barcoded libraries were generated and analyzed with Qubit dsDNA (Thermo Fisher Scientific) and Fragment Analyzer (Advanced Analytical) methods. Library pools were sequenced on Illumina NovaSeq instruments using 300-cycle kits.

Bioinformatic analysis

Viome's bioinformatic methods include quality control, strain-level taxonomic classification, microbial gene expression and human gene expression characterizations. The quality control includes per sample and per batch metrics, such as the level of barcode hopping, batch contamination, positive and negative process controls, DNase efficacy and number of reads obtained per sample. Following the quality control, the paired-end reads are aligned to a catalog containing rRNA, the human transcriptome and 53,660 genomes spanning archaea, bacteria, fungi, protozoa and viruses. Reads that map to rRNA are filtered out. Strain-level relative activities are computed from mapped reads via the expectation-maximization algorithm [44]. Relative activities at other levels of the taxonomic tree are then computed by aggregating according to taxonomic rank. Relative activities for the biological functions are computed by mapping paired-end reads to a catalog of 52,324,420 genes, quantifying gene-level relative activity with the expectation-maximization algorithm and then aggregating gene-level activity by Kyoto Encyclopedia of Genes and Genomes ortholog (KO) annotation [45]. The identified and quantified active microbial species and KOs for each sample are then used for downstream analysis.

For analysis of the human and microbial transcriptome, each sample was analyzed by two separate pipelines (for the microbiome and human transcriptome, respectively). It is possible to map reads to both pipelines.

To assess the ability of Viome's custom catalog to interrogate the oral microbiome, the authors compared the genera and species present in the Viome catalog with those present in the expanded Human Oral Microbiome Database V3, which contains 2123 genomes representing 539 species (www.homd.org/genome/genome_table). Of the 539 species, 109 are uncultured oral taxa with undefined taxonomic classification. Of the remaining 430 species, 307 (71%) are present in Viome's catalog. The presence rate is higher at the genus level, with 158 out of 191 genera (83%) present in Viome's catalog. Viome's custom catalog is sufficiently comprehensive to provide insight into the oral microbiome.

For the presentation of metrics from the saliva cohort (Supplementary Table 6), the percentage of total reads that align to archaea, bacteria, eukaryotes, viruses and humans are reported. These percentages include only reads that align to their target and meet Viome's quality-control criteria. Reads that align to rRNAs, internal control sequences, rRNA depletion probes and unaligned reads are removed.

Data analysis

Statistical parameters, including transformations and significance, are reported in the figures and figure legends. To compare pairs of samples, the authors report Spearman correlation coefficients (which are invariant to absolute expression levels of the genes and only consider the similarity of ranked expression), Pearson correlation coefficients (which compare the linear relationship between features and is sensitive to differences in expression levels) and Hellinger distance (an appropriate distance measure for compositional data). For correlations of species and KOs, the union of features between sample pairs was analyzed. For correlations of human genes, the overlap of features between sample pairs was analyzed. The multiplicative replacement method [46] was employed to deal with missing values, the data were transformed with centered log-ratio (CLR) [47] and Spearman correlation coefficients were computed. CLR transformation is commonly done to reduce false discoveries due to the compositional nature of sequencing data. CLR transformation breaks the dependence between features (i.e., transcripts) and makes data more normally distributed, reducing the impact of highly abundant transcripts that could be artificially driving high correlation values. Statistical analyses were performed in Python.

The relative abundance of transcripts was calculated by calculating the relative activity value of every molecular feature (gene, KO, taxa). The relative activity was calculated by taking the total number of mapped reads to that molecular feature, dividing it by effective sequencing depth (the total number of reads aligned to human transcripts, KOs or species that meet Viome's quality-control requirements) and then normalizing so that the sum across all features in a sample is equal to one. Multiplicative replacement was then performed to fill in the zero values.

To assess the relationship between expression bins (low, medium and high) within a sample, correlation coefficients were individually computed for each bin. The threshold between the low and medium bins was set to represent the 50th percentile, and the threshold between the medium and high bins was set to represent the 95th percentile of the CLR values of the samples. For the assessment of the lower limit of detection, the expression threshold that was needed to achieve an overall correlation of 0.6 was computed. All correlation values were above 0.6 for species and KOs regardless of thresholding, indicating that no specified lower limit of detection was needed for high-quality data. A detection threshold was needed for robust analysis of human genes, which is reported in the associated figures.

Results & discussion

Saliva cohort metrics

For the cohort of 1102 saliva samples, the average number of sequencing reads per sample was 24,353,309. On average, there were 487 species, 2539 microbial KOs and 5174 human transcripts detected per sample. See Supplementary Tables 1–5 for a rank-ordered list of all detected genera, species, strains, microbial KOs and human transcripts in the saliva cohort. See Supplementary Table 6 for a breakdown of total reads, species/KO/human gene richness and the percentage of reads aligned to microorganisms and human genes per donor. The metatranscriptomic data demonstrate that the species and KO richness are less variable between people (%CV of 14.7% for species richness and 14.2% for KO richness) than the human gene richness (%CV of 59.8% for human gene richness). Additionally, the data show that a higher proportion of reads align to microorganisms (on average, 5.9%) than human genes (on average, 2.4%).

Method precision

For longitudinal clinical studies in large populations that measure gene expression changes as a function of health and disease states, a method with high precision is extremely important, as it allows for the measurement of small changes in the levels of gene expression that are associated with or predictive of chronic disease. To determine the precision of the metatranscriptomic method, Spearman correlations were calculated for microbial taxonomy (species) and microbial functions (KOs) for four technical replicates from ten participants (Figure 2A & B). For one participant (PID-0286), one replicate was determined to be an outlier, with the log-transformed effective sequencing depth (the number of sequencing reads aligned to microbial species) being greater than three SDs away from the mean and was removed from the analysis, as it failed quality control. The method precision for microbial taxonomy and microbial functions is high, showing that the method is able to reliably reproduce the results (Figure 2A & B).

Human transcripts are present in low amounts in saliva and therefore require higher sequencing depth to analyze. Samples from four donors were randomly chosen for resequencing with higher sequencing depth to determine the method precision for human genes. The effective sequencing depth (the number of reads aligned to human transcripts) needs to be about 1 million reads to yield high-quality human gene expression data (Table 1 & Figure 2C). This demonstrates that the technology does allow for high-precision human gene expression analyses from saliva but requires higher sequencing depth.

In addition to the Spearman correlations, two random technical replicates were compared from each participant for microbial taxonomy, microbial KOs and human gene expression (Figure 3). To further categorize the correlations between technical replicates, correlations were computed for the low, medium and high expression bins (Supplementary Figure 1). The data indicated generally high correlations across expression bins. For a robust analysis of human genes, an average detection threshold of 10.25 p.p.m. was shown to be sufficient (Supplementary Figure 1C). These data show that Viome's metatranscriptomic saliva test is able to produce data with high precision.

Sample stability

RNA is prone to rapid degradation, making proper sample preservation a critical component of any metatranscriptomic test [48]. As previously described in the authors' publications, their laboratory uses a proprietary RPB to prevent RNA degradation [42,43]. The ability of RPB to preserve RNA in saliva was validated by collecting saliva from three participants and storing it at ambient temperatures (72°F, 22°C) for 0, 7 and 28 days. Subsets of samples were also packaged and were shipped via the US Postal Service to a designated address, then shipped back to the laboratory using the same shipping method to emulate an at-home testing process. Four technical replicates per condition were analyzed using the test, and microbial taxonomy and microbial functions were compared among all samples using Spearman correlations (Figure 4) and Pearson correlations (Supplementary Figure 2). The correlation between microbial species in saliva samples in all conditions was very high, with Spearman correlation coefficients above 0.928 for all conditions tested (Figure 4A). The correlation between microbial functions in saliva samples in all conditions was also very high, with Spearman correlation coefficients above 0.901 for all conditions tested (Figure 4B). Venn diagrams showing the number of overlapping and unique features between each storage condition and their expression levels can be found in Supplementary Figure 3. The data demonstrate that the majority of features overlap even across storage conditions and that the expression levels of the overlapping features are significantly higher than the unique features (p < 0.05). These data show that Viome's saliva metatranscriptomic analysis method can adequately preserve RNA for up to 28 days at ambient temperature, inclusive of sample shipping.

Longitudinal changes in the saliva transcriptome

It is important to understand the longitudinal stability of the salivary transcriptome so that large-scale studies can be reliably carried out. Toward this goal, the authors recruited eight subjects, who collected saliva samples weekly for 5 weeks. On average, 537.18 microbial species and 2693.54 microbial KOs were detected per sample. For all of the taxa and KOs that were detected, the correlation remained very high across 5 weeks for each participant (Figure 5 & Supplementary Figure 1). Correlations for one participant across expression bins (low, medium and high) are presented in Supplementary Figure 4. The correlations for each expression bin were generally high, indicating robust method performance. The correlations observed between the 1-week and 5-week samples for taxonomy ranged from 0.87 to 0.94 (Figure 5A & Supplementary Figure 1A–H) and for KOs ranged from 0.88 to 0.94 (Figure 5B & Supplementary Figure 1I–P). These data show that the saliva microbiome is longitudinally stable in terms of both composition (taxonomy) and activity (KOs).

Intrasample versus intersample precision

One of the more important parameters of the saliva test is to be sufficiently precise to distinguish very small changes among thousands of measured features. Such high precision would enable the identification of even minor transcriptome changes over time and changes related to health and disease. To assess the test precision in this context, the authors compared Hellinger distances among biological samples from the same person collected over time (intraperson distances in Figure 6) and biological samples from different people (interperson distances in Figure 6). Empirical cumulative distribution function plots of Hellinger distance show that when taking microbial species (Figure 6A), microbial functions (Figure 6B) and human genes (Figure 6C) into account, samples from one individual over 5 different weeks (intraperson distances) tend to be more similar than pairs of samples coming from two different participants (interperson distances).

Conclusion

This paper outlines a novel saliva analysis method that can produce high-quality human and microbial gene expression data. This paper also includes, to the best of our knowledge, the largest list of population-scale data from saliva samples analyzed through a metatranscriptomic method (taxonomic classifications, microbial gene expression and human gene expression). The data demonstrate clinical utility with high levels of method precision, adequate sample stability for global shipping (including to and from traditionally underrepresented communities in clinical research) and longitudinal stability suitable for population studies. This saliva method complements metatranscriptomic pipelines previously developed by our laboratory for the analysis of stool and blood [42,43]. With the integration of these tests, our lab has developed a robust and novel systems biology platform for the investigation of biomarkers critical to human health and disease.

Despite recent advancements in a variety of molecular biology techniques, limited insights have been generated into the cause of chronic and noncommunicable diseases. Systems biology approaches that take into account human and associated microbial features are rapidly becoming invaluable tools for the identification of biomarkers critical in human health and disease [37]. The saliva metatranscriptomic method described in this paper is easily deployable across the world, precise, cost effective, automated and clinically validated, and it overcomes many of the limitations previously encountered with saliva microbiome characterizations and metatranscriptomic methods (e.g., sample stability [49]). This method is suitable for use in population-scale studies to elucidate the role of the saliva microbiome and associated human gene expression in human health and disease.

Future perspective

The methods presented in this paper and other papers from our laboratory demonstrate the feasibility of cost-effective methods for interrogation of the holistic human system (microbial components and human components). We anticipate that over the course of the next 5–10 years these molecular biology techniques will be used in large-scale, population-based studies for the identification of chronic disease etiologies. This will pave the way for the development of accurate and effective diagnostic, therapeutic and preventative therapies to treat chronic diseases.

Since the methods presented by our laboratory allow for at-home sample collection, we also believe that these methods will allow for the democratization of healthcare. This will enable individuals to easily collect samples outside of primary care settings and to obtain critical health information on their own accord. We anticipate that this increased accessibility to health infrastructure will result in better and more widespread screening, which should lower the overall disease burden in the population. Additionally, we hope that shifting the healthcare setting from hospitals to the home will result in an overall reduction in healthcare costs and will enable the effective prioritization of healthcare resources.

Finally, we anticipate that the analysis of microbiomes and human gene expression will open the door for optimized personalized therapies. This personalized approach to healthcare will likely be more effective than traditional treatment strategies.