Next generation sequencing of serum and saliva RNA

RNA was isolated from serum and saliva using the miRNeasy Serum/Plasma Kit (Qiagen, Inc.) according to the manufacturer's instructions. Serum: frozen serum samples were thawed on ice, and 200μL of serum was added to 1mL of QIAzol lysis reagent. Following vigorous vortexing, 200μL of chloroform was added and the samples were incubated for 5 minutes at room temperature (RT), then centrifuged at 12,000 x g for 15 minutes at RT. The resultant aqueous phase was removed, mixed with 1.5 volumes of 100% ethanol, transferred to an RNeasy MinElute spin column, and centrifuged for 15 seconds. The column was washed with Buffers RWT and RPE at the manufacturer's indicated volumes, and the RNA was eluted with 30μL of RNase-free water. Saliva: refrigerated saliva samples originally collected in an Oragene vial or swab collection kit were incubated at 50°C for 1 hour. A 250μL aliquot was then removed, transferred to a microcentrifuge tube, incubated at 90°C for 15 minutes, and cooled to RT. 750μL of QIAzol lysis reagent was added, and the sample was vortexed vigorously for 1 minute, and incubated for 5 minutes at RT. Chloroform (200μL) was added, and the sample was vortexed for 1 minute, then centrifuged at maximum speed (>13,000 x g) for 10 minutes. 450μL of the resultant aqueous phase was transferred to a new tube, mixed with 675μL of 100% ethanol, transferred to an RNeasy MinElute spin column, and centrifuged for 15 seconds. The column was sequentially washed with Buffers RWT and RPE at the manufacturer's indicated volumes, and the RNA was eluted with 30μL of RNase-free water. RNA quality was assessed using the Agilent Technologies Bioanalyzer on the RNA Nanochip.

Stranded RNA-sequencing libraries were prepared using the TruSeq Stranded Small RNA Kit (Illumina) according to manufacturer instructions. Samples were indexed in batches of 48, with a targeted sequencing depth of 10 million reads per sample. Sequencing was performed using 36 bp single end reads on an Illumina NextSeq 500 instrument at the SUNY Molecular Analysis Core (SUNYMAC) at Upstate Medical University. FastQ files were trimmed to remove adapter sequences, and alignment performed to the mature miRbase21 database using the Shrimp2 algorithm in Partek Flow (Partek, Inc., St. Louis, MO).

The aligned reads were quantified and normalized to an internal invariant reference miRNA (miR-24-3p) and converted to log2 scale. As with the functional data, each subject’s normalized miRNA post-fight data was contrasted with their respective pre-fight/baseline values (obtained at either 1 week pre-fight or immediately prior to the fight), yielding a total of 141 sample difference values (n = 62 saliva, 79 serum). These normalized miRNA difference values were screened for sphericity using principal component analysis (PCA) prior to statistical analysis and filtered to eliminate those with more than 60% missingness.

We initially used a two-way analysis of variance (ANOVA) to examine the main effects of Sample Type and mTBI Risk (as defined by HTH) as well as their interaction to screen for miRNAs with a significant effect of the mTBI Risk on their expression level. This was performed in all of the samples from both biofluids with a False Discovery Rate (FDR) correction < 0.15. The miRNAs which passed this filter were then used in a stepwise linear regression to establish the miRNAs that best predicted the actual HTH values. A logistic regression classification analysis was then completed to assess the ability to distinguish all of the Very Likely and Low probability TBI samples from each other (holding out the Moderate group). 100-fold Monte-Carlo Cross-Validation (MCCV) was performed to estimate empirical accuracy across biofluids. miRNAs that showed the strongest predictive utility were then subjected to functional analysis using Diana Tools miRpathv3. The difference levels of miRNAs showing strong discriminatory and predictive power were also assessed in relation to various functional measures using Pearson correlation analysis.

Because the first miRNA analysis combined all the initial samples from each subject post-fight into the same TBI probability class, it was possible some miRNAs may have eluded detection if they only had acute or delayed effects at particular time points. Such temporal-dependent responses, if present, could be as meaningful as any derived from the subject binning. To reveal potential acute or delayed effects we used a General Linear Model to examine the effects of Time and Sample Type, and their interaction, on relative miRNA expression based on four different temporal bins. A total of 122 of the samples were used in this analysis. Time 1 contained samples from subjects who showed up to the MMA match but did not participate in a fight, and still provided a biofluid sample (these serve as controls for non-specific effects of the event) as well as subjects that participated in a match but experienced no hits to the head (these serve as exercise controls). Collectively, these are referred to as Time 1 HTH negative (HTH -) Controls. The remaining temporal bins were from fighters who participated in a match and received at least 2 hits to the head (HTH+). These HTH+ samples were grouped by collection time point into Time 1 HTH+ (within 1 hour post-fight), Time 2 HTH+ (2–3 days post-fight), and Time 3 HTH+ (7 days post-fight). The temporal profiles of all miRNAs with significant Time effects were visualized using line plots and subjected to supervised classification analysis to identify the most salient patterns. miRNAs with expression profiles of interest were then subjected to functional analysis using Diana Tools miRpathv3 and compared with the miRNAs from the Subject Binning analysis.

After identifying miRNAs with expression profiles of interest, we examined the balance and cognitive score data along with the molecular data using principal component analysis (PCA) to detect the molecular and functional features that show the most similarity across time. For this analysis, only miRNAs with defined patterns (Acute Saliva Response miRNAs or Delayed Serum Response miRNAs) were used with the functional data from all of the post-fight samples (n = 39 saliva, n = 31 serum). Iterative principal axis PCA was performed using a quartimax root curve extraction. Factor weights were examined to identify functional variables most similar to the miRNA variables, with line plots created for visualization purposes.

To validate the expression patterns in the NGS-based miRNA data, we elected to use a custom FirePlex assay, which has been shown to perform comparable to qRT-PCR in the validation of RNA sequencing data. A subset of the miRNAs were chosen, including 11 changed and 3 unchanged (5 increased and 3 decreased), along with the reference miRNA used for normalization of the NGS data, and a total of 65 samples (from 63 of the NGS samples with sufficient RNA remaining) were sent to the FirePlex Assay Service lab at AbCam to perform the validation in a blind fashion. Due to low yields in serum, only saliva samples were sent. Included in the samples were two pairs of stabilized saliva for comparison of the signals with purified salivary RNA.