Detection of Circulating mRNA Variants in Hepatocellular Carcinoma Patients Using Targeted RNAseq Open Access

Certain high risk-groups, such as patients with liver cirrhosis (LC), chronic hepatitis B, and Metabolic Dysfunction-Associated Steatotic Liver Disease (MASLD) with advanced fibrosis, are recommended to undergo routine surveillance for hepatocellular carcinoma (HCC) to improve early detection and clinical outcomes [1, 2]. Current standard of care consists of liver ultrasound (US) serum alpha-fetoprotein testing every 6 months [3‒5]. Both US and alpha-fetoprotein, however, have limited sensitivity for early cancers; US is also limited by operator-dependence and access challenges for some patients, and utilization remains low. A more sensitive blood test-based surveillance approach would be expected to improve uptake and increase the overall impact of HCC screening.

Liquid biopsy, utilizing blood or body fluids to search for cancer cells or other molecules such as DNA or RNA that are cancer-specific, is a promising platform to be utilized for HCC screening [6‒8]. However, technical challenges have limited broad application. For example, detection of low-frequency somatic mutations from circulating tumor cells and circulating tumor DNA has proven challenging [9‒11].

Circulating tumor RNA is an underexplored alternative approach to liquid biopsy that may offer better sensitivity for early cancer detection. Using RNAseq on RNA isolated from plasma from patients with HCC, we previously reported the detection of dysregulated mRNA transcripts in circulating exosomes also known as small extracellular vesicles (sEVs) [12]. In our previous report, NGS- and PCR-based tools were used to identify thousands of circulating mutated RNA (ctmutRNA)/variants specifically detected in the plasma from HCC patients (N = 15) compared with samples from those with LC without HCC or healthy subjects (N = 15) [13]. A set of 288 variant transcripts were identified as those with the greatest promise as early detection markers of HCC because they were (i) exclusively detected in multiple HCC samples but not in any LC samples; (ii) detected in tumors and matching plasma; and finally (iii) were “high impact” mutation lesions that would be predicted to specify aberrant or dysfunctional polypeptides. However, although provocative, the single cohort and small sample sizes and difficulty in sEV preparation limited the conclusions that could be made.

In this communication, we optimized sample preparation and the approach that was used to determine the presence or absence of the 288 ctmutRNA variants in plasma from an entirely new cohort of age- and gender-matched early stage HCC (n = 25), late-stage HCC (n = 25), and non-cancer LC (n = 31) patients. HCC tumor tissues (n = 11), normal liver tissue (n = 1), and liver cell lines (n = 3) were also included. In the present study, deep sequence analysis was also conducted for regions flanking each of the 288 candidate variant transcripts isolated from the EVs from each patient sample. Our goal was to explore a panel of circulating RNA (ctRNA) markers in LC patients not diagnosed of HCC which can be used for surveillance and screening to detect HCC at an earlier stage.

A QIAseq RNA Fusion XP system was used, and the sequence files were analyzed using the GeneGlobe (QIAseq Qiagen), a cloud-based data analysis pipeline. Variants with annotations such as SNPs, indels, and splice variants originally seen in the first samples sets were detected in many of the new patient samples. Several novel fusion transcripts associated with HCC were identified. In many cases, novel variants in the same transcript, flanking the original lesion (corresponding to the 288-test set) were detected. A panel with 36 variants was identified that could identify all HCC samples. The specific identities of the transcripts most commonly mutated in the cancer sample, their possible pathogenic significance and use as early detection and management of HCC is discussed.

Plasma samples from diagnosed HCC patients were obtained from our collaborators at University of Pennsylvania under approved IRB protocols. All plasma samples were obtained from subjects prior to receiving any treatments. HCC diagnosis was made by either (i) biopsy or (ii) typical enhancement patterns on dynamic contrast-enhanced CT or MRI (later codified as meeting Li-RADS 5 criteria). Staging was performed per the Barcelona Clinic Liver Cancer system (BCLC). We investigated separately 25 BCLC 0 or A samples (very early or early stage HCC with single tumor or up to 3 lesions <3 cm) and 25 BCLC B or C (intermediate or late-stage HCC – exceeding stage A size criteria or with vascular invasion or with metastatic disease) in order to investigate a possible discrimination between the earlier stages and the later stages of HCC. There were no stage D patients in the cohorts. Patient information is shown in Table 1. Plasma samples from LC patients without HCC (n = 31), HCC tumors (n = 11), cell lines (n = 3), and 1 normal liver were used as controls. Demographic and clinical details of patients whose tumors were investigated are shown in online supplementary Table 1 (for all online suppl. material, see https://doi.org/10.1159/000545366). The older cohort of HCC patients reported earlier [13] comprised of ∼36% stage A, 36% stage B, ∼21% stage C, and ∼7% stage D corresponding to BCLC criteria. Overall, the older and the new cohorts of HCC patients had similar stage distribution and similar age. All clinical specimens were acquired through our collaborators in University of Pennsylvania, Philadelphia, PA, Mayo Clinic, Rochester, MN and Capital Health Cancer Center, Pennington, NJ using IRB approved protocols. Plasma samples of 0.5 mL were processed using our optimized EV/RNA isolation protocol, and RNA from tumors and cell lines was isolated as reported earlier [12, 13].

Whole blood was drawn, and serum and plasma were prepared. Serum samples were prepared by incubating at room temperature for 30 min and centrifuging the samples at 2,000 g at 277 K for 10 min. Plasma was prepared using 4 different anticoagulants (EDTA, ACD, sodium-citrate, and sodium heparin) and samples were centrifuged at 2,000 g for 15 min at 277 K. Two 0.5 mL aliquots of each preparation were used to isolate RNA using the miRNeasy Serum/Plasma kit (Qiagen). A Tapestation 4200 instrument (Agilent) with HSRNA Screentape was used to analyze and quantify the RNA content in each sample.

A combination of methods was used to isolate high quality RNA from plasma samples. Plasma samples of 0.5 mL were thawed and both sEV and sEV-free fractions corresponding to each sample were processed. sEVs were precipitated using ExoQuick (System Biosciences) precipitation followed by purification using ultracentrifugation (UC) as reported earlier. The sEV-free plasma was further subjected to RNA isolation using serum plasma miRNA isolation kit (Qiagen). sEV fraction from UC was also lysed and subjected to serum/plasma miRNA Qiagen column RNA isolation. RNA from sEV and sEV-free plasma fractions was combined, and RNA samples were concentrated using a SpeedVac. QC analysis was conducted using Bioanalyzer.

EDTA plasma was collected from one individual. Increasing volume aliquots from 0.5 mL to 2.0 mL were used to isolate sEVs utilizing four different methods. Three methods use precipitation reagents, and one method is based on UC. All plasma samples were precleared of any residual cellular debris by centrifugation at 3,000 g for 15 min.

UC was performed as previously reported [13], and the sEV pellet is resuspended in 0.5 mL of PBS. sEV isolation using ExoQuick (Systems Biosciences) was performed according to the manufacturer protocol. Polyethylene glycol (PEG) has been reported to enhance EV-derived nucleic acids [14]. Test matrix 1 and 2 are made of two different average molecular weight PEGs, 6,000 and 8,000, respectively, and 0.5 m NaCl. A 2x stock solution is prepared and added in equal volume to the plasma samples and left to incubate overnight at 277 K. Samples were centrifuged at 3,250 g for 1 h at 277 K. The supernatant was aspirated, and the exosome pellet resuspended in 0.5 mL of PBS. Nanoparticle Tracking Analysis (NTA) was performed using a ZetaView (Particle Metrix) to determine the size distribution and particle count in the EV samples. RNA samples were run on the Agilent 4200 TapeStation using RNA ScreenTape. DNA contamination in RNA samples is a potential risk and to ensure that sEV/RNA isolation methods yield pure RNA, RNA isolated from plasma-derived EVs using Exoquick precipitation, exoRNeasy columns (Qiagen), and UC was either treated or left untreated with DNase before conversion to cDNA. The effectiveness of these methods was tested to measure levels of some liver specific transcripts in human circulation using qPCR.

Plasma isolated from 4 normal healthy individuals each in duplicate sets was subjected to Exoquick precipitation followed by purification of EVs by UC. The yield of EVs from 1 mL plasma samples was characterized by NTA using a ZetaView (Particle Metrix). RNA was isolated from all plasma-derived sEVs and EV-free plasma from 4 individuals using serum plasma miRNA isolation kit (Qiagen) and RNA characterized using Tapestation 4200 (Agilent).

Blood was taken from 3 healthy individuals three times over an 8-day period and EDTA-prepped plasma was prepared from each sample. Plasma samples of 1 mL with duplicate set for each time point were compared. The plasma was stored at 193 K for at least 1 day to ensure freeze and thaw consistency in results. Exoquick precipitation was performed to isolate sEVs followed by purification using UC and NTA analysis. Total RNA was isolated from sEV using the miRNeasy Serum/Plasma Kit (Qiagen) following the manufacturer’s protocols. Quantification and quality analysis were performed for each time point’s RNA using the Tapestation 4200 instrument (Agilent) using HSRNA Screentape according to the manufacturer’s published protocols.

RNA was extracted from isolated plasma-derived sEVs using the miRNeasy Serum/Plasma kit (Qiagen) and quantified on a NanoDrop One (ThermoFisher). Here, 25 ng or 50 ng of cDNA was used for each reaction. DNase I digestion was followed with cDNA synthesis. qPCR was performed in triplicate sets, in a 24 μL reaction for several gene targets on a LightCycler 480 (BioRad) for 45 cycles using PowerUP SYBR Green Master Mix (ThermoFisher). Primer information about these targets was reported earlier [12]. A synthetic DNA ultramer was used as an exogenous control to create a standard curve to calculate copy numbers [12].

To explore the 288 genomic targets of interest (online suppl. Table 2), 894 target specific primers were generated to ensure capturing all SNPs, splice variants and fusions. Primer design was based on GRCh38 genome build and designed using a proprietary algorithm to capture transcriptional variants. Briefly, QIAseq RNA Fusion XP custom panel primers were designed to flank the targeted SNP/indel site, at either the 5ʹ end, 3ʹ end, or both, depending on the target, so that the 230 bp read sequence contains sufficient endogenous sequence (≥25 bp) from both directions. High primer specificity within the Gencode Basic transcriptome model set was chosen. For SNPs and insertions/deletions (indels), chromosomal coordinates with 50 bp padding were targeted. Detailed information about target regions and primers is shown in online supplementary Table 3.

QIAseq RNA Fusion XP targeted panels were used. These require just one gene-specific primer per target allowing for considerable flexibility for targeting splicing variants, fusions, and SNP regions of interest. The QIAseq RNA Fusion XP Panels use single-primer extension (SPE) and unique molecular index (UMI) technologies in NGS to help identify and characterize fusion gene events, gene expression, and SNP/Indel at the RNA level with high efficiency, sensitivity, and flexibility. The QIAseq RNA Fusion XP Panels rely on highly efficient RNA conversion, gene-specific single-primer enrichment, and molecular barcoding for sensitive fusion, gene expression, and RNA SNP/Indel detection. Total RNA was isolated from plasma samples using Qiagen miRNA serum/plasma kit and QC performed using Tapestation. Isolated total RNA (10 ng) was reverse transcribed with high efficiency into first-strand and second-strand cDNA synthesis adapter complexes containing unique molecular indices. Sample indices are incorporated into the dsDNA. Enrichment of targets was carried out using a single gene-specific primer and a universal forward primer and the libraries were amplified via fast universal PCR, in which a second index was added. Sequencing files were fed into the QIAseq pipeline (Qiagen), a cloud-based data analysis pipeline, which enables filter, map, and align reads, as well as count unique molecular barcodes associated with targeted genomic regions, and call variants with a barcode-aware algorithm.

UMIs also called “molecular barcoding,” are applied prior to any amplification such each original target molecule is “tagged by” a unique sequence “barcode.” This is accomplished by the ligation of double-strand cDNA with a sample index adapter containing a 12-base random sequence. Statistically, this provides 4¹² = 16,777,216 unique molecular tags for each adapter and each converted double-strand cDNA molecule in the sample receives a unique UMI sequence. The barcoded cDNA molecules are then amplified by SPE for target enrichment and library amplification. Due to intrinsic noise and sequence-dependent bias, barcoded cDNA molecules may be amplified unevenly between different enriched targets. Therefore, target transcripts can be better evaluated by counting the number of UMIs in the reads rather than counting the number of total reads for each transcript. Sequence reads having distinct UMIs represent different original molecules, while sequence reads having the same UMI are the results of PCR duplication from 1 original molecule and are counted together as 1 molecule (QIAseq RNAscan handbook).

The QIAseq RNA Fusion XP Panels are provided as a single tube of primer mix, with up to 20,000 primers per tube (custom panel). RNA samples are initially converted to first-strand cDNA. A separate, second-strand synthesis is used to generate double-stranded cDNA (ds-cDNA). This ds-cDNA is then end-repaired and A-tailed in a single-tube protocol. The prepared ds-cDNAs are then ligated at their 5′ ends to a sequencing platform-specific adapter containing UMI and sample index. Adapter-ligated cDNA molecules are subject to limited target-barcode enrichment with SPE. This reaction ensures that intended targets are enriched sufficiently to be represented in the final library. Universal PCR is then carried out with highly efficient, low error rate, fast processing Taq enzyme to amplify the library.

Libraries were prepared using QIAseq RNA Fusion XP Panel (QIAGEN) according to manufacturer’s instructions from 20 ng RNA for the tissue samples and 5 μL of the concentrated plasma samples (500–1,000 pg). A total of 31 cycles were performed for the universal PCR amplification step. The libraries size distribution was validated, and quality inspected on a Bioanalyzer 2100, DNA7500 Chip (Agilent Technologies). Libraries were pooled in equimolar concentrations based on the bioanalyzer automated electrophoresis system (Agilent Technologies). The library pool was quantified using qPCR and optimal concentration of the library pool used to generate the clusters on the surface of a flow cell before sequencing on a NextSeq 550 instrument (Read 1: 229, Read 2: 69, Index 2 × 10) according to the manufacturer instructions (Illumina Inc.). Raw data were demultiplexed and FASTQ files for each sample were generated using the bcl2fastq software (Illumina inc.).

The QIAseq RNA Fusion XP Analysis workflow of GeneGlobe on https://geneglobe.qiagen.com/gb/ was used for the analysis of the samples. Panel CJHS- 14875Z-894 was chosen with default parameters. Data were captured in terms of detected variants, gene expression, and fusions. The filters are designed to catch false positive calls that have incorrectly high mutation likelihood for various reasons. A non-reference allele needs to pass the quality score threshold and all filters to be reported as a variant. smCounter2 uses a logistic regression classifier to determine if an indel in homopolymer is real. The filters in smCounter2 are: STAR Parameters Setting align SJoverhangMin:12; alignSJDBoverhangMin:12; chimSegmentMin:12; outFilterMultimapNmax:10; Filters Setting endogenous: 15; longest_must_be: 25; minmts <dynamic>; minreads <dynamic> as summarized in QIAseq Fusion XP SNV/Indel calling table online supplementary Table 4.

Variant calling was performed using GATK tool kit as reported earlier (Block T. Frontiers). Variants with 20 or higher quality scores as well as depth of 5 or more were considered as true positive. Filtration and annotation of high-risk variants were carried out using SnpEff analysis (version 4.3t) to identify high risk variants in all sample subsets [15, 16]. SnpEff is a software tool that predicts the functional effects of genetic variants, such as single nucleotide polymorphisms (SNPs), indels, and structural variants. It is used in genomic studies to determine how variants may impact genes, proteins, and phenotypes. It can predict coding effects such as synonymous or non-synonymous amino acid replacement, missense, splice site, start/stop codon gains or losses, or frameshift changes.

We first optimized methods to isolate EVs and RNA from human blood. Blood from the same healthy individual was drawn and processed with various anticoagulants or without any anticoagulant to extract serum. Total RNA was isolated from 0.5 mL plasma and serum samples using Qiagen serum/plasma miRNA isolation kit and characterized using Tapestation. Comparable amounts of RNA from 10 to 20 ng were obtained (Fig. 1a). Serum and EDTA, ACD, and sodium-citrate prepared plasma have relatively the same amount of RNA, about 12.5 ng in 0.5 mL of sample. The serum appears to have more RNAs in the 175–500 nucleotide (nt) range while most of the RNA in the plasma samples is in the 150–200 nt range. Regardless of sample type, the greatest amount of RNA is always in the small range, about 25 nt with the 75–150 nt species also prominent in plasma. While the sodium-heparin plasma appears to have much more RNA per sample, the extra steps that must be taken to make the RNA viable for downstream applications do make it a less enticing candidate.

RNA isolated from plasma and serum specimens was also treated with polynucleotide kinase, an enzyme known to stabilize circulating fragmented RNA [17]. However, this did not result in any significant change in RNA yield (online suppl. Fig. 1A).

Since RNA in the circulation is mostly concentrated in sEVs, we evaluated various sEV isolation methods. In addition to an ExoQuick precipitation and UC for sEV isolation, we also tested two matrices produced in our lab, we called matrix 1 and 2. sEVs were characterized by NTA (Fig. 1b). UC performed poorly in terms of particle count with a range of 4.1 × 10¹⁰ EVs in the 0.5 mL plasma sample to 2.0 × 10¹¹ particles in 2.0 mL of plasma. Test matrix 1 and 2 (TM1 and TM2) did result in a slightly wider range of sEV particle size, 50–300 nm, than ExoQuick which had a narrower distribution, but both did capture the most particles in the 100 nm range. To evaluate the effect of one cycle of freezing and thawing on sEV-derived RNA yield, fresh and frozen plasma samples were compared. As expected, freezing and thawing led to a 25–50% loss of RNA depending on the method of sEV isolation (online suppl. Fig. 1B). All four sEV-isolation methods, however, show comparable levels of liver specific transcripts when tested by qPCR. Beta actin transcript was used as a control (Fig. 1c). To test the purity of ctRNA isolated by EQ, UC, and ExoRNeasy approaches, ctRNA samples were either treated with DNase or left untreated before conversion to cDNA followed by analysis of liver specific transcripts using qPCR. Comparable copy numbers of both Alb and FTL transcripts were observed confirming the absence of circulating tumor DNA contamination (Fig. 1d).

Using Exoquick precipitation method of capturing sEVs from plasma, we tested the temporal stability of sEV recovery and ctRNA levels. We demonstrate comparable levels of sEVs (Fig. 2a) and ctRNA levels (Fig. 2b) derived from 1 mL plasma in 4 different healthy individuals. To test the reproducibility and precision of sEV and ctRNA isolation, we sought to investigate the EV counts and RNA yield from increasing volumes of plasma samples. All four isolation methods show a consistent increase in both particle count and RNA amount directly proportional to the increasing volume of plasma used for isolation (online suppl. Fig. 2A, B). Compared with the commercial precipitation reagent ExoQuick (Systems Biosciences), our in-lab PEG-based solutions were as effective at capturing an equivalent number of particles and in total RNA recovered. The extent to which liver-associated circulating transcripts vary from day to day was examined. Plasma acquired at 3 variable time points over a period of 8 days from 4 different normal healthy individuals was tested. Linear longitudinal RNA samples corresponding to day 1, day 3, and day 8 were characterized by Tapestation, followed by investigating the levels of four liver specific transcripts using qPCR with the results shown in Figure 2c, d. Consistent number of sEV particles (Fig. 2c) and ctRNA levels (Fig. 2d) were recovered from linear time points across 3 normal healthy individuals. No significant changes were detectable in the levels of three liver specific transcripts with respect to time in the same individual (Fig. 2e). The amplified products are presented as relative copy numbers as a function of the day on which the blood was drawn. Importantly, the relative number of copies for each transcript did not vary by more than 3-fold for any given transcript. This suggests that the amount of these transcripts is generally consistent and does not vary significantly from day to day in human plasma.

From our previous total RNAseq work [13], we identified 288 ctmutRNA targets that correlated with a diagnosis of HCC that were concordant in plasma and tumors from HCC patients but were not present in LC or “Normal Healthy Controls” (NHC) plasma. These were classified as high-risk based on snpEff prediction of a high likelihood of having a detrimental effect and showed high prevalence among HCC samples (online suppl. Table 2). Since these variants were detected by RNAseq and from a relatively small number of patients, we sought to determine their detectability in a new cohort of liver disease patients (Table 1) using a “targeted” RNAseq approach, which would characterize these variant lesions and provide more sequencing depth than did RNAseq.

Primers and probes specific for the 288 targets were produced (online suppl. Table 3). RNA was isolated from plasma-derived EVs in 81 patient samples. Additionally, RNA from 11 tumors, 3 cell lines, and 1 liver tissue from a donor source who did not have liver disease (“normal liver”) was included. The investigation of 81 plasma samples using targeted seq analysis (QIAseq RNA-FusionXP, Qiagen) of the 288 specific lesions resulted in a total of 12,877 variants. Bar graphs representing the total number of variant counts corresponding to different clinical categories are shown in Figure 3. As shown in the figure, SNPs are the largest class of circulating variants, followed by indels, a profile mirrored in tumors. Missense variants are the largest class by function, followed by synonymous and frameshift variants. On the basis of impact, the representation of high-risk variants in circulation is the lowest compared to low and moderate-risk variants.

Variant calling using a depth ≥5, quality score Q ≥20 meaning that variants of a quality score <20 and depth <5 reads were removed. A total of 3,347 variants (165 per sample) were found in early stage HCC, 3,362 (134 per sample) in late-stage HCC, and 386 variants were shared between the two HCC groups. A total of 5,112 variants in total (165 per sample) were associated with LC samples and 1,194 total (149 per sample) were found in tumor tissues. The distribution of variants and range of recurrence corresponding to various sample subsets is shown in Figure 4a.

To further identify highly pathogenic or high-risk variants, SNPeff mediated annotation was carried out. SNPeff employs sequence and structure-based bioinformatics tools to predict the effect of coding variants on the structural phenotype of proteins and classifies RNA variants as high-, moderate-, low-risk variants. A “high-risk variant” in the context of SNPeff refers to a genetic variant that, based on its predicted functional impact, is likely to have a high likelihood of causing a detrimental effect. Venn diagrams (Fig. 4b) show the distribution of high-risk variants in various clinical subgroups. A total of 438 high-risk variants were exclusively detected in HCC plasma and 390 high risk variants exclusively with LC plasma samples. Among the high-risk variants in HCC plasma, 236 were associated with early stage, 216 with late-stage HCC and 59 were shared between the two groups. Some 110 high-risk variants were associated with tumor tissues but not detectable in the normal liver tissue. In summary, a higher number of high-risk variants were identified compared to the original test set of 288 ctmutRNA targets, presumably because of the very high sequencing depth achieved.

It was first of interest to know if any of the original exact variants from the list of 288 HCC-specific test candidates recur in the new cohort of HCC samples. A set of ∼74 variants annotated to the same exact genomic location as in the 288-target panel and thus represent the same original mutation. That is, 74 variants in the current cohort were identical to those found in the original cohort. Table 2 is a list of variants from the original test panel of 288 targets which are identical and validated in the new cohort of patients and thus represent HCC selective lesions in both studies. Of particular interest because of their biological significance are variants associated with the transcripts of EPB41, M6PR, ARHGAP5, PRDX6, FASN, ARCN1, SENP7, ECI1, IST1, ACOX1, EIF3G, which had lesions identical in the original and current cohort samples of those with HCC and were not present in plasma from those without HCC.

A significant number of new variants close to and flanking the original test lesions were identified in all clinical subsets. Targeted RNAseq allowed for much greater sequencing depth than was possible with conventional RNAseq we employed earlier [13]. Therefore, the possibility that recurrent high-risk/high-impact variants exclusively detected in HCC plasma in the vicinity of the initial 288 test lesions was also explored. High-risk/high impact variants exclusively detected in HCC plasma in at least 2 samples and not detected in any of the LC patients (n = 31) were identified. These variants are represented in a heat map shown in Figure 5. Each row shows a different variant, and each column shows a different sample. A total of 236 recurrent (≥2 patients) high-risk variants were found to be associated with early stage HCC patients and 215 high-risk variants associated with late-stage HCC patients. Fifty-nine variants were shared in both early- and late-stage patients. Variant annotation, effect, minor allele frequency and recurrence of early stage and late-stage HCC patients are tabulated in online supplementary Tables 5 and 6, respectively. Among the early stage variants, 137 variants represent indels with frameshift mutations and 99 are SNPs representing stop-gained mutation. A total of 180 variants represent noncoding transcript exon variants containing regulatory elements. Among the late-stage variants, 111 represented indels and 104 variants belonged to SNPs. A total of 159 late-stage variants represent noncoding transcript exon variants, containing regulatory elements. As in the case of early stage high-risk variants, indels associated with late-stage variants represent frameshift mutations, while the SNPs predominantly represent stop-gained mutation. As the data demonstrates, several highly recurrent variant transcripts correspond to genes with established onco-pathological associations, such as EPB41, ARCN1, ECI1, CCAR1, HDAC5, SENP7, and ACOX1. The possible significance of this is considered further in the Discussion.

To identify pure HCC-specific variants, our strategy of filtration and enrichment has excluded all variants which are detectable in non-cancer LC samples. However, since HCC often develops in the setting of LC it would be expected that some early markers of HCC would be present during LC. We, therefore, looked for ctmutRNA variants which are common in the circulation of both LC and HCC patients and HCC tumors but not detectable in normal liver tissue.

A set of 790 variants identified as shared between LC patients and HCC patients are represented in a heat map (Fig. 6a). The frequency of each of these variants ranged from 3 to 77% in LC patients, 2–80% in HCC patients and 0–80% in the tumor tissues. The most highly recurrent variant is the EIF3G SNP (Chr19:101155880 A>G), noncoding transcript exon variant, with a recurrence of 80% in tumors, 77% in LC plasma, and 80% HCC-plasma samples. The variant minor allelic frequency associated with this variant is very high in tumor and plasma samples, suggesting that it is a common variant. Similarly, an indel variant IST1 (Chr 16:71922608 AATGCCC>A, disruptive in frame deletion was recurrent in 67% LC plasma, 62% HCC plasma, and 75% of tumors. The annotation and proportions of LC-associated high-risk variants are shown in online supplementary Table 7. Among these 790 variants prevalent in non-cancer LC and HCC patients, 335 variants were identified as high-risk with varying recurrences in LC, HCC, and tumor samples. In summary, we demonstrate the detectability of highly recurrent circulating variants shared between LC and HCC patients. The possibility that these recurrent variants offer potential in identifying high risk chronic liver disease patients is considered in the Discussion.

Apart from variants detected in circulation, we wanted to find out the most recurrent variants associated with tumor tissues. Eleven tumor tissues from HCC patients (online suppl. Table 1) were investigated by targeted seq analysis. All 288 ctmutRNA test candidates were observed to be concordant between HCC-plasma samples and HCC tumor tissues [13]. Therefore, in addition to plasma samples, 11 HCC tumors, 1 normal liver tissue, the hepatoblastoma HepG2, hepatoma Huh7, and “normal immortalized” liver PH5CH cell lines were also investigated using the targeted RNAseq panel. A set of 204 variants not detected in normal liver tissue and with ≥0.1% recurrence in tumors are represented in a heat map in Figure 6b. As we have previously demonstrated, tumor tissues reflect a higher density of ctmutRNA variants than plasma samples from HCC patients. There are several variants that are enriched in tumors and HCC cell lines, yet not detectable in normal liver tissue. Some highly recurrent tumor tissue specific variants, also detectable in HCC and LC plasma samples, are summarized in Table 3. Characteristics, annotation, and recurrence of tumor tissue associated variants with ≥0.1% recurrence and not detected in normal liver tissue are shown in online supplementary Table 8.

In the previous analysis, we focused on identifying only high-risk variants or variants of highly pathogenic impact. Since the definition of “risk” (using snpEff annotation) is a judgment based on bioinformatic prediction, it is possible that variants not identified or “filtered” as “high risk” could also be significant and disease-correlative. We also observed that there are variants which are predominantly detectable in HCC plasma but also present in a few non-cancer LC plasma samples. We relaxed filtration frequency in order to investigate variants that are detectable in most of the HCC patients but rarely detected in LC samples. The ctmutRNA variants that were most common in and selective for the HCC samples, regardless of predetermined “risk” or “impact” are represented, with their prevalence, in Table 4. The two SNPs corresponding to ASGR1 and MEX3C transcripts, for example, were detected in 24% (or 12) of the HCC samples. The ASGR1 SNP was not detected in any of the LC samples but the SNP corresponding to MEX3C was present in 1 of the non-HCC (LC) samples. At the other extreme, a SNP in the CHMP2A transcript was present in only 3 HCC samples, but in none of the non-HCC (LC) samples.

In an effort to develop a panel with the highest diagnostic precision, we used various combinations of the ctmutRNA shown in Table 4 to generate 3 panels of variants that varied in their sensitivity and specificity for distinguishing the HCC from LC patient samples. The results are shown in Figure 7. In summary, all 12,877 variants are reduced to a panel of 23 in the following way. First, we enforce a specificity criterion that ≤0/31 LC patients and a sensitivity criterion that at least one early- or late-stage patient should be identified. Second, we remove the variants which identify the same set of early- and late-stage patients as any other variants. Third, we sort out the remaining variants in the decreasing order of the number of early- and late-stage patients they identify. Lastly, we go down this list of sorted variants and keep only those variants that identify additional patients than the variants higher up on this list. The panel with all 23 ctmutRNA candidates from Table 4 had 100% sensitivity (detected all 50 HCC patient samples) and 100% specificity (not detected in any of the 31 non-HCC samples). The number of ctmutRNA candidates in the panel could be reduced to, for example, to 18 or 17, and though 100% of the HCC samples were still detected, specificity is compromised as 1 or 2 LC (non-HCC) samples, respectively, would also reflect the variants. Online supplementary Table 9 provides more details about the panel construction and testing. In summary, we demonstrate a panel of circulating mRNA variants exhibiting high diagnostic specificity to identify HCC patients.

The main objective of this study was to validate the occurrence and detection of 288 circulating HCC-specific RNA variants [13] for consideration of use in HCC surveillance. In this work, we demonstrate detection and recurrence of high-risk circulating RNA variants in the circulation of a large cohort of HCC patients. Most of the ctmutRNA appear to reside within sEVs confirming earlier studies in which the observation of hundreds of mRNA variants in blood and tumors of HCC patients was first reported [13]. It extends those findings by profiling particular variants in a new cohort of early- and late-stage HCC patients with the creation of a panel of variants that can be used to detect 100% of the cancer samples and distinguish them from the non-HCC samples.

We also describe a few simple and practical methods of isolation of the circulating sEVs containing ctmutRNA without differentiating EV subclasses or their tissues of origin. Importantly, RNA transcripts are generally stable in the circulation upon phlebotomy, carried out at different time points from the same individual. Freezing plasma samples after blood draw led to a marginal reduction in the RNA yield compared to fresh plasma. These data suggest that, surprisingly, RNA can be a robust blood derived analyte for analyzing both expression and mutation profiles.

In the current study, at least 75 of the original 288 variants distinguished in the first cohort were also detected in samples from the current independent cohort of 50 HCC patient samples. These ctmutRNA variants were again not present in any of the original cohort of non-cancer patient samples (online suppl. Table 1). However, there are some variants with high prevalence in our new cohort of HCC samples that are also detected with very low prevalence in new cohort of non-HCC subjects (Table 2). The most commonly occurring HCC selective variants that are exact matches with the original 288 variants with their prevalence in the HCC samples correspond to GOLGA2 ∼20%; CHMP2A, ARCN1, and SLC25A39 ∼16%; SIRT2 and PTBP1 ∼15%; HDAC5, GET3, PPP1R12C ∼14%; EPB41, CXCL3, GPS2, and TRMT1 ∼12%; FASN, BIRC2, EIF3G, UBR4, GET3, TRIP12, NDUFS8 all at ∼11%, etc. Detection of these specific variants in the circulation of CLD patients could be useful in identifying high-risk patients. Future work will probe the functional consequences of these specific mutations on the dysregulation of the targets and molecular pathways.

As a point of extra validation, in another unrelated study which is underway in our laboratory, we are employing both exome sequencing and RNAseq to compare matching specimens of tumor, morphologically normal adjacent tissue and plasma from 5 diagnosed HCC patients. Preliminary RNAseq results from this small cohort confirm high impact somatic mutations in transcripts like ARCN1, CHD4, GDI1, MKNK2, EPB41, INF2, GPS2, HNRNPL, FASN, MRPL4, LRP1, GRN, PTPA, UBR4, NUMA1, and LRRFIP1 which are shown in Table 2 portraying 75 recurring variants in this study.

One might wonder why only 75 out of the original 288 high-risk candidates were detected in a new HCC cohort using targeted RNAseq approach which achieves higher sequencing depth and high confidence in calling a variant. This could be due to well-known genetic heterogeneity of HCC [18, 19]. Similar arguments could explain the low overlap between early- and late-stage HCC patients. In addition, all variants detected in LC patients (recurrence rate 3–77%) and tumor tissues (recurrence rate 9–100%) are all accounted for without any filtration. However, for HCC-specificity any variant detected in HCC samples which was also detected in LC patients was removed. Also, to ensure HCC specificity, if a variant in the original list of 288 HCC-specific variants was detected in any new LC controls it would be filtered out to ensure HCC specificity. Also, although the original cohort of HCC patients was comparable to this new validation cohort, differences in clinical characteristics between two cohorts might also have contributed to low overlap suggesting that the original list of 288 variants may still hold relevance. Technical differences between total RNAseq and targeted RNAseq could also be responsible. Finally, we cannot entirely exclude the possibility that some of the variants among the 288 list are bioinformatics tool artefacts as these pipelines are constantly being improved and updated. Since the objective of the current study is to validate the findings from the discovery cohort, a 26% validation rate is still significant.

For our purpose, targeted RNAseq allowed for the detection of a large number of additional variants within the vicinity of intended target lesions. Many of these are considered “high-risk” mutations, implying they cause a significant pathological impact on the expression or alter the structure of the corresponding protein.

Understanding basic aberrations in the early stages of cancer can reveal progressive changes from a healthy condition to a malignant state. Cumulative accumulation of these dysregulations marks the initiation of cancer. Aberrations associated with RNA processing and alternative splicing have been reported to contribute to cancer development and progression [20, 21]. Cataloging these alterations in RNA may provide early detection tools as well as highlight aberrant pathways and therapeutic vulnerabilities. Both SNPs and indels have been implicated in Mendelian and complex diseases. SNPs involve a single nucleotide change whereas an indel incorporates or removes one or more nucleotides. It is not clear whether indels in RNA are more likely to influence complex and pathogenic traits than the more abundant SNPs. Since reading frames should be maintained to preserve protein function, coding indels are subject to stronger purifying selection than SNPs [22, 23]. The cumulative contribution of indels compared with SNPs to oncological risk is not known.

Though HCC development can occur without signs of cirrhosis, in most cases, however, it develops in the context of LC. During the transition of the disease from LC to HCC-specific somatic, alterations likely occur at or near the advent of HCC. In such patients, we expect both LC specific as well as HCC-specific signatures. The strategy we have employed in our work is to use profiles of LC without HCC patients as controls to identify highly specific signatures associated with HCC. We demonstrate variants that are specific to each clinical category and others that appear to overlap (Fig. 4b). In Figure 6 heat map, we demonstrate variants which are common in tumors, HCC-plasma and non-cancer LC plasma samples. However, if a variant occurs in samples from those with a diagnosis of HCC as well as in those with only LC and no HCC, it may complicate early HCC detection because it might occur well before imaging can detect cancer. That said, it occurs to us that by ignoring or disqualifying variants detected in both HCC and LC, we may be missing the most early markers. Similarly, variants detected in the tumor tissues would be cancer specific and could be used to increase specificity of the assay. However, as we demonstrate some variants derived from tumor tissues are also detected in LC patients without HCC. We have to be careful enough since this type of very early HCC-derived variant might proceed all other methods of detection, including advanced imaging. Therefore, we understand that caution is needed in its use.

Some of the HCC-specific mRNA variants in the circulation of HCC patients we identified have previously been demonstrated to be associated with cancer based on cancer databases like catalog of somatic mutations in cancer , dbSNP, Clinvar and portal GDC from NCBI and NCI, IntOgen and others which document the mutations detected in human tumor tissues. Among the top six high-risk HCC associated variants we found in early stage HCC patients (with each detectable in at least 2 patients) include three indels associated with MYO1C, CCAR1, and ARHGAP45 transcripts, all resulting in frameshift mutations. Three additional indels in ARHGAP45 flanking the original test lesion were exclusively detected in early stage HCC samples leading to high-risk frameshift mutations. The CCAR1 transcript is a prognostic marker in liver, ovarian, and renal cancer [24] and exhibited an additional high-risk variant known to cause frameshift duplication in early stage HCC. Early stage HCC patients showed indels in SENP7 and TRIP12 resulting in high-risk frameshift mutations. Consistent with this, indels corresponding to SENP7 and TRIP12 were also detected in 75% and 12% tumor tissues, respectively. We detect high impact mutations in ATF2, RNF44, NUP54, LRRFIP1, SLC25A39, HDAC5, CDK4, and HERC1 in the circulation of early stage HCC patients [24‒26].

Interestingly, differential expression analysis of HERC1 showed significant upregulation in early stage HCC plasma compared to late-stage HCC patients (data not shown). We also identify high-risk SNP variants corresponding to POLR2E, ABTB1, DBNL, SH2B3, METTL26, and TRMT1 transcripts and annotated as stop gained and splice variants, associated with early stage HCC. Several other high-risk indels and frameshift variants corresponding to POLR2E, MKNK2, ARCN1, ABCD3, TPP1, SRP72, CUEDC2, PRDX6, PPP1R12C, ZNF691, FAH, and KHSRP transcripts, known prognostic markers of liver and other cancers [24, 26‒28] were observed to be associated with early stage HCC patients. Yet another high-risk frameshift deletion detected in ODC1, an FDA approved drug target and involved in polyamine metabolism [24, 29] is associated with early stage HCC. With further understanding of the role and physiological impacts of these high-risk variants on early cancer development, the findings are likely to have profound impact on current efforts toward development of early detection biomarkers for HCC.

Although more than ∼300 high impact variants were observed both in early- and late-stage HCC patients, after filtering out variants also detected in non-cancer LC patients, only a small fraction of HCC-specific variants was shared between early- and late-stage HCC patients. One can speculate the role of tumor heterogeneity which also suggests that HCC tumors progress through independent pathways. Though tumor staging reflects the size of the tumor and extent of spread, tumor grade essentially refers to its aggressiveness and degree of differentiation which is dictated by the dysregulated signaling pathways. High-risk variants associated with SENP7, TRIP12, SH2B3, METTL26, and TRMT1 are shared between early stage and late-stage patients. Late-stage HCC patients exhibit some unique alterations compared to early stage patients. Among the highly recurring high-risk variants in late-stage HCC patients are indels corresponding to EPB41, PNRC1, HNRNPL, MKNK2, EIF2S2, RBM42, PTTG1IP, SRM, PRRC2C, ZNF691, RPF1, GDI1, MAP4K2, SLC9A3R2, C16ORF72, and SLC9A3R2 known to be prognostic cancer markers [24]. ctRNA transcripts M6PR and GPAA1 harbor high-risk SNPs in late-stage HCC patients annotated as missense splice and stop-gained variant, respectively. In addition to the utility as diagnostic tools, late-stage biomarkers can be used in the assessment of tumor response to therapy as it permits a prospective end-point evaluation and provides a guide for clinicians to make future treatment decisions. Both early detection and the ability to assess the tumor response to treatment are critical aspects in the field of cancer.

LC is the dominant risk factor for HCC and was associated with 2.4% of global deaths in 2019. The burden of cirrhosis remains substantial owing to underdiagnosis and undertreatment of chronic liver disease, and the number of deaths and cases of cirrhosis are projected to rise in the next decade [30]. We reiterate that our focus is to identify circulating mRNA variants which facilitate identifying HCC patients among CLD patients with and without cirrhosis. Since we have employed LC patients without HCC as controls for ensuring HCC biomarker specificity, any variant which is detected in both HCC patients and LC patients without HCC is filtered out from HCC panels. This study, however, also highlights a repertoire of variants which are not only shared between the plasma samples from HCC and non-cancer LC samples but were also detected in 11 tumors studied here. Some of these variants may possibly originate from tumors but because they are detected in the circulation of LC patients without cancer, they are not useful as surveillance markers for HCC. Circulating ctmutRNA associated with LC offers a unique opportunity to facilitate early detection of cirrhosis and to reduce its global burden. In our study, analysis of circulating variants associated with both LC and HCC patients identify some variants with high recurrence and allelic frequency. Variants corresponding to EIF3G, IST1, P2RY11, EIF4EBP2, HNRNPL, CTSW, and PRDX6 showed a recurrence between 48 and 77% in LC patients and 40–80% in HCC plasma and 50–100% in tumor tissues. These overlapping variants in LC and HCC may not be critical as HCC biomarkers but one can speculate that they would serve as excellent circulating LC biomarkers because of high recurrence and high variant allele frequency.

Our objective of identifying altered circulating transcripts with SNP and indel variants correlating with HCC could be well extended to RNA fusion variants. Growing evidence suggests that RNA-level fusions can be involved in multiple cancer-related processes or can provide novel biomarkers for diagnosis and prognosis [31]. Transcriptomics and matched whole-genome sequencing data from tumors of 1,188 individuals of the Pan-Cancer Analysis of Whole Genomes Consortium revealed that 18% of fusions displayed no evidence of genomic rearrangement [32]. Though the test panel of 288 targets was not designed to characterize any fusion variants, we identified novel fusion variants in the vicinity of 288 ctmutRNA targets in both LC as well as HCC samples, which need deeper investigation (online suppl. Fig. 3A–C; online suppl. Table 10).

In conclusion, we demonstrate the detectability of HCC associated high-risk ctmutRNA variants in the circulation of HCC patients, which when validated in a larger set of patients could serve as exquisite non-invasive biomarkers for early detection of HCC. Determination of the usefulness of a cancer biomarker for early detection involves multiple stages from discovery to utility validation, as described in [33]. Our study represents an early, but significant step in this process. With further understanding of the role and physiological impacts of these high-risk variants on early cancer development, the findings are likely to have profound impact on current efforts toward development of early detection biomarkers for HCC.

We appreciate the role of Nishi Patel at UPEN in acquiring HCC and LC plasma samples and Tina Bergner from Capital Health in procuring HCC tumor tissues. We acknowledge the support of Maria Middleberg and Lee Middleberg from Serologix Inc. for help with phlebotomy procedures.

This study protocol was reviewed and approved by (1) University of Pennsylvania Institutional Review Board, Philadelphia, PA, Approval No. 828260, (2) Institutional Review Board Capital Health, Capital Health Medical Center, Pennington, NJ, Approval No. CH800, and (3) Institutional Review Board Mayo Clinic, Rochester, MN, Approval No. 21-003568. Written informed consent was obtained from patients whose specimens were used in this study under approved IRB protocols from University of Pennsylvania, Philadelphia, Mayo Clinic, Rochester, and Capital Health cancer center, NJ.

Timothy M. Block is on the board of the Directors of Hepion, a company developing medicines for liver diseases, Merlin, a company developing medicines for cancer and CIRNA, a startup company focused on development of tools for early cancer detection. Both Merlyn and CIRNA are owned by Blumberg institute. Aejaz Sayeed is the CSO for CIRNA and consults for Exocel Bio. Lewis R. Roberts was a member of the journal’s Editorial Board at the time of submission.

This study was supported by a Commonwealth University research enhancement program grant with the Pennsylvania Department of Health, SAP# 4100095321 (A.Sa.). The department disclaims responsibility for any analyses, interpretation, or conclusions. NIH NCATS grant to Medical University of South Carolina A00-2219-5010 (T.M.B.); VA Merit I01-CX-001933, I01-CX-000988 (D.E.K.). This work was also supported by private funding sources viz Blake Foundation, Steve Miller foundation (A.Si.) and Angel Investment Capital Program from Commonwealth of PA to CIRNA through Blumberg Institute.

Aejaz Sayeed (A.Sa.) conceptualized, designed the study, and drafted the manuscript. T.M.B. helped analyze the results and helped draft the manuscript. A.Sa. and D.Z. developed methodology, planned and carried out experiments, analyzed data, and provided interpretive insights. Angela Simeone carried out RNAseq analyses. A.Sa., D.Z., T.Z., and Angela Simeone created the figures. T.Z. played a role in biostatistical analyses to derive variant panels. D.E.K., M.A.H., L.R.R., C.D., and F.Y.A. provided clinical specimens, reviewed the data, and provided critical insights. All authors contributed to the article and approved the submitted version.

All processed data in our study has been submitted in the manuscript in the form of figures and online supplementary tables. Raw unprocessed FASTQ files will be uploaded and accessed via a publicly available NCBI SRA database upon acceptance of the manuscript. Further inquiries can be directed to the corresponding author.