Examining potential confounding factors in gene expression analysis of human saliva and identifying potential housekeeping genes

Isolation of RNA from whole saliva, a non-invasive and easily accessible biofluid that is an attractive alternative to blood for high-throughput biodosimetry of radiological/nuclear victims might be of clinical significance for prediction and diagnosis of disease. In a previous analysis of 12 human samples we identified two challenges to measuring gene expression from total RNA: (1) the fraction of human RNA in whole saliva was low and (2) the bacterial contamination was overwhelming. To overcome these challenges, we performed selective cDNA synthesis for human RNA species only by employing poly(A)+-tail primers followed by qRT-PCR. In the current study, this approach was independently validated on 91 samples from 61 healthy donors. Additionally, we used the ratio of human to bacterial RNA to adjust the input RNA to include equal amounts of human RNA across all samples before cDNA synthesis, which then ensured comparable analysis using the same base human input material. Furthermore, we examined relative levels of ten known housekeeping genes, and assessed inter- and intra-individual differences in 61 salivary RNA isolates, while considering effects of demographical factors (e.g. sex, age), epidemiological factors comprising social habits (e.g. alcohol, cigarette consumption), oral hygiene (e.g. flossing, mouthwash), previous radiological diagnostic procedures (e.g. number of CT-scans) and saliva collection time (circadian periodic). Total human RNA amounts appeared significantly associated with age only (P ≤ 0.02). None of the chosen housekeeping genes showed significant circadian periodicity and either did not associate or were weakly associated with the 24 confounders examined, with one exception, 60% of genes were altered by mouthwash. ATP6, ACTB and B2M represented genes with the highest mean baseline expression (Ct-values ≤ 30) and were detected in all samples. Combining these housekeeping genes for normalization purposes did not decrease inter-individual variance, but increased the robustness. In summary, our work addresses critical confounders and provides important information for the successful examination of gene expression in human whole saliva.


Materials and methods
Sample collection. Whole saliva samples were collected using ORAgene®RNA (RE-100, catalog number: RE-100) vial collection kits from DNA Genotek according to the manufacturer's instructions (DNA Genotek Inc., Kanata, Ontario, Canada). The kit is an all-in-one system for the collection, stabilization and transportation of RNA from saliva. Unstimulated whole saliva was collected from 61 healthy donors (27 females, 34 males, average age 38.5 ± 16.4 years, Fig. 1). The following exclusion criteria were applied: age below 18 years old, donors with history of immunodeficiency, autoimmune disorders, viral hepatitis, HIV infection, current or previous cancer, current oral problems (infection). From all donors, whole saliva was collected from 9 to 11 am and then preserved after having been shaken vigorously in the vial. Additionally, in 10 out of 61 donors, whole saliva samples were collected at 9 am, 3 pm, 9 pm and 9 am again the next day. After completing normal oral hygiene, donors were not allowed to eat or smoke 2 h prior to collection or to drink at least 1 h prior to collection. Samples were stored at room temperature overnight and placed in a freezer (− 20 °C) for storage. All samples were anonymized and obtained with informed consent from the donors. Sampling was carried out in accordance with the institutional guidelines and regulations. At the time of sampling, donors were asked to fill in a questionnaire (supplemental Figure 1) about demographics (e.g. sex, age), social habits (e.g. alcohol, smoking), oral hygiene (e.g. flossing, mouth wash), and previous radiological diagnostic procedures (e.g. number of CT-scans).
RNA extraction. Total RNA comprising a mixture of human and bacterial RNA, was isolated from whole saliva samples following a combination of the ORAgene® RNA purification protocol 16 and the mirVana™ kit protocol (Invitrogen™, ThermoFisher Scientific, Carlsbad, CA 92008; USA/Life Technologies, Darmstadt, Germany) as described in detail elsewhere 15 . In brief, the samples were heated at 50 °C (1 h), three aliquots (of 1000 µl) were generated, incubated at 90 °C (15 min), cooled to room temperature, 40 µl ORAgene® neutralizer solution (1/25 of total volume) was added, incubated on ice, centrifuged at 13,000g (3 min) and the cell-free clear supernatant was collected for further processing. At this step, we switched to the mirVana™ kit protocol 17 by adding the Lysis/ Binding Solution. With the mirVana™ kit, total RNA, including human and bacterial RNA species, was isolated by combining a Phenol-Chloroform RNA precipitation with further processing using silica membranes. After several washing procedures to purify RNA from other residual debris, DNA residuals were digested on the membrane (RNAse-free DNAse Set, Qiagen, Hilden, Germany). RNA was eluted with 100 µl RNAse free water in a collection tube and the aliquots were pooled for each sample. In order to increase the input RNA amount for downstream gene expression analysis, samples were steamed at 45 °C for 90 min followed by elution with 30 µl of RNase free water before freezing at − 20 °C.
Quality and quantity of isolated total RNA were measured spectrophotometrically using NanoDrop™ One Microvolume UV-Vis spectrophotometer (NanoDrop, PeqLab Biotechnology, Erlangen, Germany). RNA integrity was assessed by the 2100 Agilent Bioanalyzer (Life Science Group, Penzberg, Germany) and DNA contamination was controlled by conventional PCR using actin primers.
Conventional cDNA synthesis-high-capacity cDNA reverse transcription kit. For analyzing gene expression of human rRNA (18S) and pan-bacterial rRNA (16S, see below), total salivary RNA was converted into complementary DNA (cDNA) via reverse transcription using the High-capacity cDNA reverse tran- Adjustment of human RNA input for cDNA synthesis. A main modification of the previously described workflow 15 was the adaption of human RNA input for cDNA synthesis with the SuperScript® III First-Strand Synthesis System (Fig. 2). As NanoDrop™ spectrophotometer only provides total RNA values comprising an unknown mixture of human and bacterial RNA species, we calculated the ratio of the detected raw Ct-values (threshold cycles) of human rRNA (18S) to pan-bacterial rRNA (16S, high-capacity cDNA reverse transcription kit cDNA synthesis followed by 1st qRT-PCR), regarding rRNA as representative for all human and bacterial RNA specimens. The ratio was determined by calculating the relative ratio of 18S rRNA to 16S rRNA for each sample individually as follows: www.nature.com/scientificreports/ Using the generated ratio, a RNA input of 4 ng total human RNA was determined individually for each sample for a second cDNA synthesis followed by 2nd qRT-PCR (confirmation for methodological reasons in this context, Fig. 2

boxes in darker grey).
Modified cDNA synthesis-SuperScript III First-Strand Synthesis System. In order to subsequently perform gene expression analysis of human origin, only eukaryotic messenger RNA (mRNA) was reverse transcribed to cDNA by using the SuperScript® III First-Strand Synthesis System with Oligo (dT) 20 primers. As a result, pan-bacterial RNA and other human RNA species were excluded from further processing, minimizing the inhibition effects of those RNA species on the following reactions. According to the kit description, the amount of starting material can vary from 1 pg to 5 μg of RNA and the maximum input volume of RNA is 8 µl 19 . The amount of human RNA input was set to 4 ng per reaction. Using the concentration from repeated NanoDrop™ measurements and the calculated 18S/16S ratio in each sample, the corresponding amount of measured total RNA input was calculated, conforming to the maximum input volume of 8 µl. The RT was performed according to the manufacturer's instructions 19 .

Pre-amplification-TaqMan® PreAmp Master Mix.
To detect low abundance mRNA species, preamplification was required. We used TaqMan® PreAmp Master Mix (Thermo Fisher Scientific Baltics UAB, Vilnius, Lithuania) to increase the amount of specific cDNA targets, synthesized with the SuperScript® III First-Strand Synthesis System. According to the manufacturer, pre-amplification with this kit is linear when a minimum amount of cDNA molecules is present (minimum of 1-250 ng and Ct-values without pre-amplification should be < 35) and multiplex amplification can be performed by pooling up to 100 TaqMan® Gene Expression Assays. Pre-amplification was performed according to the TaqMan® PreAmp master mix kit protocol 20 . In the present work, 10 different TaqMan® Gene Expression Assays (ACTB, Hs01060665_g1; B2M, Hs00187842_m1; GUSB, Hs00939627_m1; MT-ATP6, Hs02596862_g1; PGK1, Hs00943178_g1; PP1A, Hs99999904_m1; RPL13A, Hs04194366_g1; RPLP0, Hs02992885_s1; TBP, Hs00427620_m1; YWHAZ, Hs01122445_g1) were utilized and pooled to enable the multiplex amplification of specific cDNA targets. Those are commonly used house-keeping genes already employed in other experimental set ups (www. genom ics-online. com).

Real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR).
Using different sets of primers, two kinds of cDNAs were utilized for qRT-PCR: for human (18S rRNA, Hs99999901_g1) and pan-bacterial (16S rRNA, Ba04230899_s1) primer probe designs, cDNA from High-capacity cDNA reverse transcription kit was used, whereas for the primer probe designs representing the potential house-keeping genes (ACTB, B2M, GUSB, MT-ATP6, PGK1, PP1A, RPL13A, RPLP0, TBP, YWHAZ), SuperScript™ III First-Strand Synthesis SuperMix synthesized, i.e. human cDNA with and without 14× pre-amplification was used for detection of each of these genes in each sample. The qRT-PCR reaction contained TaqMan® Universal PCR Master Mix and one of the inventoried TaqMan® Gene Expression Assays for separate detection of transcripts. All meas-   were calculated for continuous variables such as gene expression data and age. Frequency tables of categorical data were examined for statistical differences using the chi-square test for equal proportion. Comparisons of categorical variables with gene expression values were performed using the non-parametric Kruskal-Wallis (KW) test. We assessed the assumptions of normality (Kolmogorov-Smirnov) and boxcox transformed the continuous variable "age". All calculations were performed using SAS (release 9.4, Cary NC, USA). Graphical presentations were performed using Sigma Plot 14 (Jandel Scientific, Erkrath, Germany).
Data availability and approval statement. The merged set of raw data is provided within supplemental  .7) implied on average about 2,702 (2 28.5-17.1 ) times more bacterial RNA copy numbers relative to human RNA in whole saliva (Fig. 3A). This means that for each copy of a human gene, on average 2,702-times more copies of bacterial genes can be found in the samples. Human 18S rRNA raw Ct values showed a broad variance with almost 10 Ct values in the 50% interquartile range, indicating about 1000-fold differences in RNA copy numbers (Fig. 3A).
Adjustment of human RNA input for cDNA-synthesis and its confirmation. Using the generated ratio (ratio = 2^( Ct18S rRNA-Ct16S rRNA) ), a defined amount of human RNA (4 ng) was reverse transcribed in a second cDNA synthesis, again using High-capacity cDNA reverse transcription kit, followed by qRT-PCR analysis for detection of human 18S rRNA and bacterial 16S rRNA Ct-values (Fig. 3B). The 50% interquartile range of human 18S rRNA quantification in 51 saliva samples changed significantly (P < 0.001) from 7 Ct values of unadjusted input RNA (comprising an unknown mixture of human and bacterial RNA) to 1.5 Ct values after considering the degree of bacterial contamination and adjusting for it.

Identifying inter-and intra-individual differences in saliva RNA isolates. Sociodemographic and
epidemiological characteristic of donors. The 61 healthy donors comprised almost equal proportions of females (44.3%) and males (55.7%, Table 1). Donors were mostly Caucasians (77.1%) and about half of them were aged 19-30 years at time of saliva collection. Only 18% smoked over at least 5 years and another 16% of former smoker smoked at least 2 years and stopped smoking about a year ago. Alcohol consumption and diet was reported by 41% and 28%, respectively. Most (83.6%) of the donors brushed their teeth at least twice a day and about half of them used flossing and about one-third mouthwash. Braces (3.3%) and dentures (11.5%) were reported less frequently and oral problems like periodontitis in 23%. Acute and chronic diseases such as rheumatism or disease of the thyroid gland were reported in 6.6% and 16.4%, respectively. Radiological examinations during the last six months (mainly X-rays and CT-scans) were reported in 18% and none of our donors received radiotherapy. www.nature.com/scientificreports/ group showed significantly more alcohol consumption per week (P = 0.012), higher frequency of former smoker (P = 0.0058), radiological examinations (P = 0.012) and dentures (P = 0.032). In donor 1, we observed a minimum human 18S rRNA raw Ct value of 21.7 (9 am the next day) and a maximum 18S rRNA raw Ct value of 35.1 (3 pm), resulting in a delta Ct value of more than 13 or in other words more than 8,000fold difference in human RNA amount (data not shown). Relative to the earliest sampling time at 9 am, slightly increased median RIN values and about two-fold decreased interquartile ranges were observed at 3 and 9 pm (Fig. 4B, P = 0.022), but distribution of RIN values at 9 am of the following day were similar to the 9 am values of the previous day (Fig. 4B).  Table 2).

Identifying inter-and intra-individual differences in housekeeping genes. Baseline gene expres-
Correlation of housekeeping gene expression with sociodemographic and epidemiological characteristics. Eight of ten housekeeping genes showed a weak (P = 0.02-0.05) but significant association of altered gene expression and oral hygiene like mouth wash and flossing (Table 2). Further weak associations were found for alcohol consumption with PGK1 (P = 0.03), age with RPLPO (P = 0.034) and radiological examination during the last 6 months with ATP6 gene expression (P = 0.046). Only TBP showed no significant association with sociodemographic and epidemiological characteristics ( Table 2).

Correlation of housekeeping gene expression with a circadian periodic of saliva collection.
None of the housekeeping genes revealed significant gene expression changes associated with the time of saliva sampling, but different patterns in gene expression changes over time of sampling were observed for individuals (Fig. 6). For example, we found a decrease of normalized Ct values among all genes in donor 9 at 9 am on the next day,  . The box plots in (B) represent the human 18S rRNA raw Ct-values before and after adjustments accounting for input differences from left to right. The left part shows Ct values from 1st qRT-PCR performed using cDNA with an input of 500 ng total RNA (bacterial and human RNA) and the right boxplot the corresponding results when taking 4 ng of human RNA (calculated via the 18S/16S-ratio together with total RNA concentration values measured). Asterisks (**) refer to a P value < 0.001 using 500 ng total RNA measurements as reference.  (Fig. 6).

Discussion
The search for simple less-invasive sampling methods plays an important role for high-throughput diagnostic tests like for victims of radio/nuclear incidents. Besides whole blood, saliva, a non-invasive easily-accessible biofluid, has been shown to contain mRNA biomarkers for prediction and diagnosis of several diseases 21,22 . Examinations in this regard are challenged, because our previous studies indicate that most of the isolated RNA originates from the oral microbiome, thus, reducing the amount of isolated human RNA considerably. We previously modified the methodology to better analyse the low abundance of the human RNA fraction from whole saliva 15 .  www.nature.com/scientificreports/ In the current study, we confirmed previous results increasing the sample size from 12 to 91 whole saliva samples. Furthermore, we modified the previously described workflow to ensure equal human RNA input for cDNA-synthesis as a prerequisite for comparability among samples when performing quantitative RT-PCR. In addition to normalization using a housekeeping gene, this normalization step via RNA quantification proved to be a robust method when measuring RNA expression between samples 23 . Spectrophotometrically, only total RNA (non-specific human and bacterial) can be measured. Measurement of total RNA quantity is relatively uninformative considering the high and inhomogeneous bacterial contamination of saliva samples. To resolve this issue, we relatively quantified human RNA using 18S rRNA as surrogate and bacterial contamination using 16S rRNA as surrogate, and introduced a correction factor for the same starting material (human RNA amount) for downstream cDNA-synthesis. By performing a second cDNA-synthesis with the High-capacity Kit (followed by qRT-PCR with detection of human 18S rRNA and bacterial 16S rRNA) we confirmed that we were able to adjust the amount of human RNA input for downstream cDNA synthesis with the SuperScript® III First-Strand Synthesis System (amount of starting material: 1 pg-5 µg 19 ) to equal amounts.
In the current study we also examined factors with potential impact on RNA quantity and quality. Those included different collection time points for saliva sampling or extraneous factors, demographic and epidemiological. Concerning circadian periodic rhythmicity, published protocols from other research groups recommend highest yields of RNA as well as best quality when sampling between 9 and 11 am 9,24-28 . Our study indicated no significant differences in RNA amounts between the collection time points, but intra-individual differences were high (> 1000-fold). Also, RIN values slightly improved at later (3 pm and 9 pm) relative to early collection time points (9 am), indicating that saliva samples can be collected during the whole day, thus, widening the applicability of this approach, e.g. for clinical use.
We did not observe any significant differences concerning sociodemographic and epidemiological characteristics that would explain the observed magnitude of inter-individual variance of human RNA yields. Furthermore, the differences in amounts of human RNA (raw Ct value for 18S rRNA was ≤ 30) between the samples could not be explained by sociodemographic or epidemiological characteristics. In this study, the addressed sociodemographic and epidemiological conditions seemed to be of minor relevance for interpretation of saliva gene expression results.
Human 18S rRNA, as a commonly known housekeeping gene, cannot be used as a normalizer in gene expression analysis in the current application due to the lack of a poly(A)+-tail 15 . We examined the baseline gene expression values of ten commonly used housekeeping genes. These genes appeared not or only weakly altered by sociodemographic or epidemiological factors which adds to their robustness. ATP6, ACTB and B2M appeared www.nature.com/scientificreports/ most suitable, because sufficiently high copy numbers ensured detection in all samples indicating methodological robustness. However, inter-individual differences in gene expression and certain time points were found in several donors. These patterns occurred among all genes, suggesting they may be caused by methodological reasons. Combining all three housekeeping genes in our study did not reduce the variance significantly, but normalization based on more than one reference gene has been increasingly suggested by others [29][30][31] . Certainly, combining ATP6, ACTB and B2M as housekeeping genes for expression studies using human saliva will increase the robustness and, therefore, would be suggested. Nevertheless, planned future studies on saliva samples from irradiated donors will finally show whether radiation impacts these three identified genes, which would render them unsuitable as housekeeping genes. The applicability of suggested housekeeping genes has to be proven in future independent studies. Finally, some limitations of this manuscript need to be considered. The major limitation lies in the fact that epidemiological data was gathered from the healthy donors by anamnesis. For example, the oral health status was self-reported and not confirmed by medical examination.
Nevertheless, a main strength of this study is that it represents a comprehensive examination of different facets when working with human whole saliva. The enhancement of the methodology and the proper examination of potential confounders like the influence of sociodemographic and epidemiologic characteristics that could potentially influence salivary isolates are completely novel findings, not described before in the literature. Considering saliva as an emerging source of body fluid for gene expression examinations underlines the importance of those findings. A key strength of the present study was the sample size: 91 samples from 61 donors in total are remarkable numbers considering molecular biological studies. These numbers and the numerous endpoints in this study are sufficient for creating reliable hypotheses.
In summary, we (I) improved the comparability of gene expression measurements among different saliva samples, (II) demonstrated that quality and quantity of RNA isolates is highly robust considering potential confounding factors such as demographics/epidemiologic and the saliva sampling time, making the approach of saliva collection even more attractive for further biomarker studies and (III) identified a set of potential housekeeping genes (ATP6, ACTB and B2M) and suggested their combination to increase robustness of saliva-based gene expression studies.