The effect of the stress hormone cortisol on the metatranscriptome of the oral microbiome

Imbalances of the microbiome, also referred to as microbial dysbiosis, could lead to a series of different diseases. One factor that has been shown to lead to dysbiosis of the microbiome is exposure to psychological stressors. Throughout evolution microorganisms of the human microbiome have developed systems for sensing host-associated signals such as hormones associated with those stressors, enabling them to recognize essential changes in their environment, thus changing their expression gene profile to fit the needs of the new environment. The most widely accepted theory explaining the ability of hormones to affect the outcome of an infection involves the suppression of the immune system. Commensal microbiota is involved in stressor-induced immunomodulation, but other biological effects are not yet known. Here we present the impact that cortisol had on the community-wide transcriptome of the oral community. We used a metatranscriptomic approach to obtain first insights into the metabolic changes induced by this stress hormone as well as which members of the oral microbiome respond to the presence of cortisol in the environment. Our findings show that the stress hormone cortisol directly induces shifts in the gene expression profiles of the oral microbiome that reproduce results found in the profiles of expression of periodontal disease and its progression.

vortexing and aliquoted in 1ml volume per well in a 24-well Corning™ Costar™ Flat Bottom Cell Culture Plate (Thermo Fisher).
To three of those wells, we added hydrocortisone-water soluble (Sigma-Aldrich) to a final cortisol concentration of 3.89g/ml (109l per well of a dilution 1/100 from a stock solution of 38mg/ml), while we did not add anything to the other three wells that were used as controls. The plate was incubated at 37 o C for 2 hours under anaerobic conditions. Cells were collected by centrifugation at 10,000 x g for 8 minutes, and RNA was extracted immediately for further analysis as described below.

Effect of cortisol on Leptotrichia goodfellowii and Fusobacterium nucleatum
We also challenged L. goodfellowii and F. nucleatum with cortisol to assess the effect that this hormone has on their expression profiles. Both organisms were grown on 10mL of BD BBL Schaedler broth (Thermo Fisher) at 37 o C under anaerobic conditions until reaching an OD 600 = 0.8. Cells were collected by centrifugation at 10,000 x g for 8 minutes, washed and resuspended in 10 ml of modified saliva medium (MSM) as described by Pratten et al. 2 . We used artificial saliva medium to avoid variability among actual saliva samples from volunteers that could influence the results. MSM composition was: yeast extract 2 g/l (BD Biosciences), proteose peptone 8 g/l (Sigma-Aldrich), porcine gastric mucin 2.8 g/l (Sigma-Aldrich), sodium chloride 0.2 g/l, potassium chloride 0.2 g/l (Sigma-Aldrich), calcium chloride 0.3 g/l; 1.28 ml/l of a 0.2 9m filter-sterilized solution of 40% urea was added after autoclaving. MSM medium was allowed to stay at least 60 h under anaerobic conditions before inoculation.
1ml of the suspension per well was added to a 24-well Corning™ Costar™ Flat Bottom Cell Culture Plate (Thermo Fisher). To three of those wells, we added hydrocortisone/cortisol (Sigma-Aldrich) to a final cortisol concentration of 3.89g/ml while three wells with the suspension of the microorganisms but without cortisol were used as controls. The plate was incubated at 37 o C for 2 hours under anaerobic conditions. Cells were collected by centrifugation at 10,000 x g for 8 minutes, and RNA was extracted immediately for further analysis as described below. In the case of these pure cultures, RNA was not amplified for analysis. The rest of the procedure is identical to the one used for metatranscriptomic analysis of the whole community.

Metatranscriptomic analysis.
Detailed protocols for community RNA extraction, RNA amplification, and Illumina Sequencing are described in Yost et al. 3 . Briefly, 6009L of mirVana kit lysis/binding buffer and 300 9l of 0.1-mm zirconia-silica beads (BioSpec Products) were added to the samples. Samples were bead beaten for 1 min at maximum speed. RNA was extracted following the protocol of the mirVana TM Isolation kit for RNA. MICROBExpress (Life Technologies) was used to remove prokaryotic rRNA. All kits were used following the manufacturer's instructions. RNA amplification was performed on total bacterial RNA using MessageAmp TM II-Bacteria RNA amplification kit (Life Technologies) following the manufacturer's instructions. Sequencing was performed at the Forsyth Institute.
For the bioinformatic analysis, we first calculated the biological variation (BCV) after estimating the common dispersion using the R package edgeR. These values were used as a cv (coefficient of variation) cutoff in NOISeq. For the whole community, analysis cv was 2 (200 cutoff NOISeq), for the Leptotrichia goodfellowii libraries cv was 1.2 (120 cutoff NOISeq) and for the Fusobacterium nucleatum libraries cv was 0.28 (28 cutoff NOISeq).
Low-quality sequences were removed from the query files. Fast clipper and fastq quality filter from the Fastx-toolkit (http://hannonlab.cshl.edu/fastx toolkit/) were used to remove short sequences with a quality score >20 in >80% of the sequence. Cleaned files were then aligned against the bacterial/archaeal database using bowtie2. We generated a .gff file to map hits to different regions in the genomes of our database. Read counts from the SAM files were obtained using bedtools multicov from bedtools 4 .
To identify deferentially expressed (DE) genes from the RNA libraries, we applied non-parametric tests to the normalized counts using the NOISeqBio function of the R package 'NOISeq' with 'tmm' normalization, with batch and length correction and removing genes whose sum of hits across samples was lower than 10. We used a significance threshold value of q=0.98, which is equivalent to an FDR adjusted p-value of 0.08 8 .
To evaluate functional activities deferentially represented we mapped the DE genes to Gene Ontology (GO) terms (http://www.geneontology.org/). GO terms for the different ORFs were obtained from the PATRIC database (http://patricbrc.org/portal/portal/patric/Home). GO terms not present in the PATRIC database and whose annotation was obtained from the HOMD database or the J. Craig Venter Institute were acquired using the program blast2GO under the default settings 6 . Enrichment analysis on these sets was performed using the R package 'GOseq,' which accounts for biases due to overdetection of long and highly expressed transcripts 6 . Gene sets with ≤ ten genes were excluded from analysis. We used the REVIGO web page 7 to summarize and remove redundant GO terms. Only GO terms with FDR adjusted p-value < 0.08 in the 'GOseq' analysis were used.

Phylogenetic assignment of transcripts (LEfSe)
Counts from the mRNA libraries were used to determine their phylogenetic composition for bacteria and archaea. Phylogenetic profiles of the metatranscriptomes were obtained using Kraken 8 .
We generated a custom Kraken library with the oral microbiome genomes indicated in the above section with a filtering threshold of 0.08. Phylogenetic profiles were used to identify significant differences between active communities under the different conditions studied by performing linear discriminant analysis (LDA) effect size (LEfSe) as proposed by Segata et al. 9 with an alpha value for the Wilcoxon test to 0.01. Significant P-values associated with microbial clades and functions identified by LEfSe were corrected for multiple hypothesis testing using the Benjamini and Hochberg false discovery rate correction 10 using the p.adjust function in R with a cutoff of FDR < 0.08.