Salivary miRNAs as non-invasive biomarkers of hepatocellular carcinoma: a pilot study

Patient recruitment and sample collection

Saliva samples were collected from 20 individuals with HCC and 19 from individuals with cirrhosis seen at the Cleveland Clinic (Cleveland, OH). Participants were adult patients (>18 years of age) who underwent liver transplantation for HCC, surgical resection for liver tumors or liver biopsy. Some patients with HCC had previously received treatment but all had active HCC at the time of collection. Prior treatments and clinical staging for patients with HCC can be found in Table S1. A description of the cohort is provided in Table 1, including the frequencies of chronic liver disease (CLD) between the two groups. Twelve out of the twenty patients with HCC had CLD, and analyses were performed with these samples combined and stratified by liver disease, as described below. Initial disease diagnoses were made from a combination of clinical presentation, imaging and laboratory techniques. Subsequently, these diagnoses underwent a secondary confirmatory pathological diagnosis. All participants provided written informed consent and the study was approved by the Cleveland Clinic IRB (IRB #10-347).

Small RNA-seq extraction and sequencing

Small RNA library preps were prepared using the QIAseq miRNA Library Kit (QIAGEN, Hilden, Germany). Adapters are first ligated sequentially to the 3′ and 5′ end of the miRNAs followed by cDNA synthesis with UMI assignment, cDNA cleanup, library amplification and final library cleanup. All protocol steps were followed based on the use of the miRNeasy Serum/Plasma kit used upstream for purification of RNA, which has been shown to be effective for saliva samples (Zahran et al., 2015)⁠. The recommended starting amount of total RNA is 5 µl of the RNA eluate when 200 µl of sample has been processed using the miRNAeasy Serum/Plasma Kit. Adapter dilutions throughout the protocol followed the serum/plasma recommendations. Cycles of library amplification followed that of a 10 ng input sample with 19 total cycles, consistent with manufacturer recommendations. Two RNA control samples (Human XpressRef RNA, QIAGEN, Hilden, Germany) with an input of 10 ng RNA were processed alongside the saliva samples. Final libraries were validated by Qubit Fluorometer (Invitrogen, Waltham, MA, USA) and Fragment Analyzer (Agilent Technologies, Inc, Santa Clara, CA, USA), and quantified via qPCR using NEBNext Library Quant Kit for Illumina (New England BioLabs, Inc, Ipswich, MA, USA). Pooled libraries were diluted, denatured and loaded onto the Illumina NextSeq 550 System, following the NextSeq User Guide. All 42 libraries were sequenced on one NextSeq High Output flow cell, single read 75 cycle run. FastQ files were developed for downstream analysis.

Data processing

FastQ files were first evaluated for read quality, using FastQC v0.11 (Andrews, 2010)⁠. Individual reads were then trimmed for phred quality scores (Q > 20) and presence of adapters, using fastp v0.19 software (Chen et al., 2018)⁠. Processed FastQ files were then aligned to the human genome (hg38) (Schneider et al., 2017)⁠ and annotated using known miRNA sequences from miRBase v22 (Griffiths-Jones, 2010)⁠. Alignments and miRNA counts were performed using the Rsubread package v2.2⁠, following the authors’ guidelines for miRNA. Mature miRNA and miRNA hairpin precursors were both used for alignment, and are designated “hsa-miR” and “hsa-mir”, respectively. Data files, including miRNA counts are located in Files S1 and S2, and R code for the analysis described below can be found at: https://github.com/rotroff-lab/salivary-miRNA-HCC. Raw sequences can be found at: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA755300.

Differential expression analysis

miRNA expression was compared between HCC and cirrhosis samples using the R package DESeq2 v1.28 (Love, Huber & Anders, 2014). First, raw read counts were normalized to account for sequencing depth, gene length, and RNA composition. Wald’s test was then performed on the normalized counts to identify any differentially expressed miRNA, using the cirrhosis samples as reference. All association tests included age, sex, race, body mass index and smoking status as model covariates to prevent confounding. In addition, patients were stratified to perform comparisons between those with CLD and the cirrhosis reference samples to distinguish miRNA associations specific to HCC vs those that may be due to impaired liver function. Log₂ fold changes (logFC) were estimated using empirical Bayes procedure and implemented with the R package apeglm v1.10 (Zhu, Ibrahim & Love, 2019)⁠. All P values were false discovery rate (FDR) corrected (Benjamini & Hochberg, 1995)⁠. Any miRNA with an FDR P < 0.05 was considered significantly differentially expressed. The level of confidence for significant miRNAs were categorized based on read counts per million (CPM) into low, medium or high categories. miRNAs with <1 CPM in >20% of the samples were considered to be low confidence. miRNAs were classified as medium confidence if at least 20% of the samples had CPM ≥ 1 and high confidence if at least 20% of the samples had CPM ≥ 10. Only miRNAs meeting the criteria for medium or high confidence were considered for inclusion in the predictive models, as described below. Pathway analysis was performed using QIAGEN Ingenuity Pathway Analysis (IPA) software (Krämer et al., 2014). Here, differentially expressed miRNA (FDR P < 0.05) were mapped to mRNA and tested for overrepresentation in cancer and disease-specific pathways using the IPA knowledgebase. The analysis was performed (1) restricted to only cancer-related pathways and (2) cancer-related and disease-specific pathways combined.

Comparison to tissue-based miRNA profiles

Differentially abundant miRNAs were compared to microarray expression data from Martinez-Quetglas et al. (2016), GEO Accession Number: GSE74618 (Martinez-Quetglas et al., 2016)⁠. These data included miRNA expression data from 218 samples from human HCC tumors and 10 samples from cirrhotic non-tumoral tissue. Differential expression was conducted on the normalized expression values between HCC and cirrhosis samples, using the GEO2R tool available from the GEO database (Barrett et al., 2012)⁠. Overlapping miRNAs detected in both GSE74618 and the salivary miRNA data were compared and were adjusted using a FDR approach, as described above.

Predictive modeling for biomarker development

Model development was performed using the statistical software, R (Team, 2020). Significant miRNAs (FDR P < 0.05) were assessed for their ability to differentiate HCC and cirrhosis samples. Here, we used support vector machine (SVM) models with a radial basis function using the R caret package (Kuhn, 2008). Because some patients with HCC had no evidence of CLD, we also performed modeling on the subset of patients with CLD to compare with the cirrhosis group. Hyperparameters for both models with and without covariates are listed in Table S2. miRNA expression was mean centered and scaled by the standard deviation prior to modeling. Leave-one-out cross-validation was used to limit model over-fitting. Recursive feature selection with leave-one-out cross-validation was used to characterize feature importance and 10 miRNAs were selected for inclusion in the final models (Kuhn, 2008). Model performance was assessed based on the AUC, sensitivity, specificity, positive predictive value (PPV), negative predictive value (NPV), and balanced accuracy. PPV is the ratio of true positives out of all identified positives, whereas NPV is the ratio of true negatives out of all identified negatives.

Differential expression of miRNAs in the saliva of HCC patients

A total of 4,565 miRNAs were detected and 365 miRNAs were significantly differentially expressed between HCC and cirrhosis samples and met the threshold for medium or high confidence (FDR P < 0.05), demonstrating a broad shift in miRNA profiles (Table S3).

The majority of the differentially expressed miRNAs (N = 283) were downregulated (logFC < 0) in the HCC samples compared to cirrhosis samples, and 50 of those had a logFC < −2 (Fig. 1A). The five most significantly differentially expressed miRNAs that also met the threshold for high confidence were: has-mir-3198-2 (logFC = −5.00, FDR P = 6.81 × 10⁻⁶), hsa-miR-3198 (logFC = −5.01, logFC = 6.80 × 10⁻⁶), hsa-mir-1246 (logFC = −6.03, FDR P = 7.98 × 10⁻⁶), hsa-miR-1246 (logFC = −6.86, FDR P = 7.98 × 10⁻⁶) and hsa-mir-3648-2 (logFC = −7.10, FDR P = 1.24 × 10⁻⁵). The 20 most significantly differentially expressed miRNA are shown in Table 2. A heatmap showing the expression of the significant miRNAs (FDR P < 0.05 and absolute |logFC| > 2) is shown in Fig. 2. Interestingly, a cluster of majority HCC patients was identified (N = 14 out of 16, 87.5%). The other four clusters had cirrhosis patients in majority. The number of cirrhosis patients in the four clusters was three (100%), three (75%), four (80%) and seven (58%), respectively. Notably, a cluster of three individuals from the cirrhosis group clustered together with elevated miRNA expression (Fig. 2). The Child-Pugh scores for these individuals were between five and seven, consistent with the rest of the control group, and no demographic or physical characteristics were identified that could explain the observed pattern (Fig. 2).

Pathway analysis revealed 96 (27%) differentially expressed miRNAs that were represented in at least one significant pathway in the IPA metabolic and signaling pathways knowledgebase (Table S4). Approximately a third of these miRNAs were associated with cancer-related diseases and pathways (N = 31 of 96) (Table S4). When the analysis was limited to cancer signaling pathways, the top three diseases were “early-stage invasive cervical squamous cell carcinoma” (miRNA N = 12, FDR P = 8.04 × 10⁻¹³), “hypopharyngeal squamous cell carcinoma” (miRNA N = 12, FDR P = 1.51 × 10⁻¹¹), and “early stage solid tumor” (miRNA N = 13, FDR P = 4.59 × 10⁻⁸). The cancer pathway with the greatest number of overlapping miRNAs differentially expressed in this cohort was “mammary tumor” (miRNA N = 35, FDR P = 2.26 × 10⁻²).

Comparison between patients with and without chronic liver disease (CLD)

Twelve out of 20 patient samples with HCC also had CLD, and these samples were compared with cirrhosis control samples. This resulted in 281 differentially expressed miRNAs, of which 248 (88%) were also significant in the analysis between all HCC samples and those with cirrhosis only (Table S5), suggesting that these differentially expressed miRNAs are likely to be due to the presence of HCC rather than the presence or lack of CLD. Out of 281 differentially expressed miRNA, 240 (85%) were downregulated (logFC < 0) in the HCC samples with CLD compared to only those with cirrhosis, and 49 of these had a logFC < −2 (Fig. 1B). The five most significantly differentially expressed miRNAs were: hsa-mir-122 (logFC = −6.65, FDR P = 7.54 × 10⁻⁵), hsa-mir-122b (logFC = −6.64, FDR P = 7.54 × 10⁻⁵), hsa-miR-122-5p (logFC = −6.64, FDR P = 7.54 × 10⁻⁵), hsa-mir-122b-3p (logFC = −6.64, FDR P = 7.54 × 10⁻⁵) and hsa-mir-603 (logFC = −4.86, FDR P = 1.14 × 10⁻⁴). These top five miRNAs were also considered to be of high confidence based on the CPM criterion. Pathway analysis after stratifying patients based on CLD status resulted in the same top rankings of cancer-related and disease-specific pathways (Table S6).

Comparison to tissue-based miRNA

Three significant miRNAs detected in saliva were also found to be differentially expressed in HCC tissue samples compared to cirrhotic liver tissue samples (FDR P < 0.05) (GSE74618) (Fig. 3 and Table S7). Forty-seven miRNAs were significantly different between HCC and cirrhosis tissue samples (FDR P < 0.05), 27 of which were evaluated in the saliva; whereas out of the 365 significant salivary miRNAs, 100 were also evaluated in tissue. However, only three miRNAs were significantly different in both cohorts (FDR P < 0.05) (Fig. 3). Out of the three significant miRNA found in each dataset, all were downregulated in saliva, but two were downregulated and one was upregulated in tissue (Fig. 3, Table 3). The miRNA comparisons are provided in Table S7. The significant miRNAs common to both datasets were hsa-mir-92b, hsa-mir-548i-2 and hsa-mir-548l.

Saliva miRNA predictive modeling

The predictive model trained on HCC +/− CLD vs cirrhosis optimized with the following 10 miRNAs: hsa-mir-576, hsa-miR-576-5p, hsa-mir-6727, hsa-mir-27b, hsa-miR-27b-3p, hsa-mir-4664, hsa-mir-125a, hsa-miR-6727-5p, hsa-mir-190b and hsa-miR-125a-5p. All of the above miRNAs met the threshold for high confidence except for hsa-mir-4664 and hsa-miR-6727-5p, which met the threshold for medium confidence. The model classified HCC with salivary miRNA alone (AUC = 0.74), and was improved with the inclusion of demographic and lifestyle variables (AUC = 0.87) (Fig. 4, Table 4). The model trained on the subset of patients with both HCC and CLD vs cirrhosis optimized with following miRNAs: hsa-miR-1262, hsa-mir-1262, hsa-mir-216a, hsa-mir-484, hsa-mir-30d, hsa-miR-216a-5p, hsa-miR-30d-5p, hsa-miR-484, hsa-mir-10401 and hsa-miR-454-3p. Five miRNAs above reached the threshold for high confidence: hsa-mir-30d, hsa-miR-30d-5p, hsa-miR-454-3p, hsa-mir-484 and hsa-miR-484. The model performed well, discriminating patients with HCC using salivary miRNA (AUC = 0.78) and with the inclusion of covariates (AUC = 0.87) (Fig. 4). Table 4 shows the AUC, sensitivity, specificity, and balanced accuracy for the predictive models, demonstrating that salivary miRNAs show potential for predicting which patients with CLD are likely to have HCC.