RNA sequencing analysis reveals the competing endogenous RNAs interplay in resected liver cancer patients who received interferon-alpha therapy

Interferon-alpha (IFN-α) is a general therapeutic regimen to be utilized in HCC. However, regulatory mechanisms of IFN-α on competing endogenous RNAs (ceRNAs) level in HCC are rarely understood. HCC patients with and without IFN-α treatment were calculated to analyze the expression prole of mRNA, lncRNA, miRNA, and circRNA by RNA sequence, and signicant differential expression (DE) of these types of RNAs were screened out for the following analysis. A ceRNA regulatory network was constructed to explore effects of IFN-α intervention on HCC. Furthermore, the potential prognostic factors among these DE RNAs were identied.

. Thus, adjuvant therapy is increasing attractive and becomes the research focus. Interferon-alpha (IFNα) is an antiviral cytokine which is with known anti-proliferative and immunomodulatory properties against malignant tumors [5], and has been published that could inhibit tumor metastasis in nude mice bearing human HCC xenografts [6], and in our clinical trials [7]. However, the characteristics of patients who can bene t from IFN-α are not consistent.
Yang and his team have reported that higher IFIT3 expression in tumor tissue than normal tissue can be regarded as a factor to predict response of IFN-α therapy [8]. Retinoic acid-inducible gene-I (RIG-I), another interferon-stimulated gene, has also proven as a predictive factor. Namely, patients with low RIG-I expression always presented poorer response and shorter survival to IFN-α therapy [9]. Meanwhile, we found that dihydropyrimidine dehydrogenase (DPYD), a pyrimidine catabolic enzyme, was dosedependently downregulated by IFN-α in HCC tissues, and might be a potential prognostic biomarker and a therapeutic target for HCC [10]. Although many research have screened the potential patients bene ted from IFN-α therapy, no agreement is reached.
Competing endogenous RNA (ceRNA) network is an entrant which can indirectly regulate mRNAs via competitively binding miRNAs. And in this theory, ceRNA refers to all the transcripts which can be used as the targets of miRNAs, namely including circular RNA (circRNA), pseudogene RNA, and long non-coding RNA (lncRNA) [11]. The mechanisms on regulating IFN-α mRNA level by ceRNA network have been widely reported. Speci cally, Tominori and his colleagues have published that according to binding with miR-1270, the ceRNA including IFN-α1 antisense RNA (AS) and cellular mRNAs precisely maintains the level of the type-I IFN and functions on the innate immune system [12]. In vitro experiments also have indicated that several lncRNA were differentially expressed to induce the IFN pathway under the IFN-stimulation [13]. In addition, ceRNA has also been researched in the area of liver cancer. Based on The Cancer Genome Atlas (TCGA) database, Zhao et al has discovered that a HCC-related ceRNA network which can be used as the prognostic factor to reveal cancer progression [14]. Another larger-clinical-sample research further indicated that the ceRNA network might aid to elucidate the mechanism of HCC pathogenesis [15].
Although as previously mentioned, less results have been announced to probe the mechanism on ceRNA regulating the IFN-associated pathway in HCC based on clinical samples.
Currently, utilizing the clinical samples, we try to probe the possible mechanism on how the IFN-α functions on cancer cells based on RNA sequencing. TIMER and GEPIA was used to further analyze the correction between MARCH3 expression and tumor-in ltrating in HCC. The results may lay the foundation on better comprehending the functional process of IFN-α treatment on HCC, and screening potentially bene cial patients who received the IFN-α therapy.

Patients
Totally, 4 patients received surgical treatment at Fudan University Shanghai Cancer Center (Shanghai, China) were enrolled in this study. Two weeks after surgery, four samples from these patients were collected before interferon-alpha treatment (group A). While the other four samples were immediately collected from these resected patients injected with interferon-alpha (30 mg twice a week; recombinant human interferon a-2b; Shenzhen Kexing Bioengineering Co., Ltd.) after two weeks (group B). The clinicpathological characteristics of these patients are shown in Supplementary Table 1. Prior patient consent were obtained, and the study was approved by the Ethics Committee of Fudan University Shanghai Cancer Center.
Isolation of peripheral blood mononuclear cell (PBMC) and sequencing PBMCs were respectively isolated from the peripheral blood of all the participants based on the previous method with some modi cations [16]. Firstly, the peripheral blood was layered and centrifuged (950g, 30min). Then the Ficoll-Histopaque layer was collected and stored for the following experiments. The trypan blue exclusion test was utilized to assay the cell viability.
Total RNA was isolated from the extracted PBMC by Trizol reagent (Invitrogen, USA) following the manufacturer's protocol. With the aid of 2100 Bioanalyzer (Agilent, USA) and NanoDrop Spectrophotometer (NanoDrop, USA), the RNA quality was detected.

RNA library construction and sequencing
The TrueSeq small RNA library prep kit (Illumina, USA) was utilized to construct the RNA (RNA) library followed the manufacturer's protocol. Brie y, after removing the ribosomal RNA, mRNA, miRNA and circRNAs were respectively fragmented and synthesized cDNAs. After purifying, connecting with adaptors, they was used as the templates, and the single-end sequencing was worked on the HiSeq4000 sequencer (Illumina, USA).
Date processing and differentially expressed RNAs screening Raw reads were screened with an in-house pipeline which had been internally validated. To acquire the clean reads, adaptor, low quality, and high unknown bases were trimmed. Genome GRCh38 was utilized as the human reference background and the analyzing tool was Bowtie [17] or Tophet [18]. The expression levels of RNAs were counted by using the reads per kilobase transcriptome per million mapped reads (RPKM) method. The differentially expressed (DE) RNAs (miRNAs, mRNAs, or circRNAs) between samples were identi ed by the DEGseq or Cu inks package under the R environment. The cut-off cirteria was designed as log2 (Fold_change) > 2 with the p value < 0.01. After that, Cytoscape software (version 3.7.0) was utilized to construct the ceRNA network.

Bioinformatics analysis
Firstly, the miRWalk (http://mirwalk.umm.uni-heidelberg.de/) and TargetScan (http://www.targetscan.org/mamm_31/) software were respectively utilized to predict the target genes and target circRNAs for the DEmiRNAs. The starbased website (http://starbase.sysu.edu.cn/) was utilized to predict the target genes for DEcircRNAs. Then Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway were analyzed on these DERNAs, including DEmRNAs and the target genes of DEmiRNAs and DEcircRNAs. Merely when P value < 0.05, the analysis results were regarded as owning the statistical signi cance.
Secondly, the possible functional relationships between DERNAs were further probed. Signi cantly DEmRNAs and DEcircRNAs were analyzed the Pearson correlation coe cient to determine whether presented co-expressed, and the Pearson correlation coe cient index was designed over than 0.9 with an adjusted p value < 0.1. The whole ceRNA network including the miRNA-circRNA regulatory network, the miRNA-target gene regulatory network, and the circRNA-mRNA co-expression network were constructed.
TIMER database analysis TIMER is a comprehensive resource for systematic analysis of immune in ltrates across diverse cancer types (https://cistrome.shinyapps.io/timer/). It includes 10,897 samples across 32 cancer types from The Cancer Genome Atlas (TCGA) to estimate the abundance of immune in ltrates. We used TIMER to investigate the correlations between MARCH3 expression with the abundance of immune in ltrates, including B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells, via gene modules. In addition, we also analyzed the correlations between MARCH3 expression and gene markers of the different tumor in ltrating immune cells. The relative gene markers are ever reported in literature.
Gene correlation analysis in GEPIA GEPIA (http://gepia.cancer-pku.cn/) is an online tool for analyzing the RNA sequencing expression data of 9,736 tumors and 8,587 normal samples from the Cancer Genome Atlas and the Genotype Tissue Expression (GTEx) projects. We further con rmed the correlation of MARCH3 with T cell exhaustion via analyzing GEPIA dataset. The Spearman method was used to determine the correlation coe cient.

Statistical analysis
Baseline characteristics of the two groups were compared using the Chi square test or Fisher's exact test. Survival curves were depicted by the "Survival" package under the R environment with the Kaplan-Meier method. The survival data originated from the cBio Cancer Genomics Portal (https://www.cbioportal.org/). The p value < 0.05 was regarded as the cut-off threshold for signi cant meaning. All statistical analyses were performed using SPSS 20.0 software (IBM, USA).

Patient's characteristics
The clinical characteristics of the enrolled population were listed in Supplementary Table 1. All the 4 patients were diagnosed as the hepatocellular carcinoma (HCC). The median age for all patients was 52 years. All of the patients had HBV viral infection, and the sex was all male. All the patients had a normal BMI and good ECOG performance status. No patients had received neoadjuvant therapy, and all the patients presented negative surgical margin.

Functional annotation of the ceRNA network
The GO and KEGG analysis were performed to understand the biological function in the DE-ceRNA network. The result of GO analysis was shown in Figure 3. Firstly, for circRNA-miRNA-mRNA network, the genes were mainly enriched in synapse, external side of the plasma membrane, and juxtaparanode region of axon; the most enriched GO terms in molecular function (MF) were chemokine receptor activity, lipoprotein transporter activity and P2Y1 nucleotide receptor binding; in additional, humoral immune response, vocalization behavior, chemokine-mediated signaling pathway and immune response were the mainly enriched biological process (BP). For lncRNA-miRNA-mRNA network, the enriched terms of cellular component (CC), MF and BP were basically similar with that of the circRNA-miRNA-mRNA network.
Then, we performed the pathway enrichment analysis of the ceRNA network (Fig. 4a, b). The up-regulated RNAs networks including the circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA were both enriched in the ve pathways, namely hematopoietic cell lineage, cell adhesion molecules (CAMs), chemokine signaling pathway, cytokine-cytokine receptor interaction, and in uenza A. And the down-regulated RNAs network mainly mapped on the following pathways, including nicotinate and nicotinamide metabolism, asthma, and allograft rejection.

Survival analysis of pivotal genes
We utilized the data from TCGA to perform the survival analysis of pivotal genes (Fig. 5a, b, c). Only three genes among the DE-ceRNA regulatory network were related with the prognosis. Speci cally, the results indicated that patients with FAM20A high expression presented a signi cantly longer survival time than the rest (P=0.0012), and the same results were also observed in patients with IGFBP4 (P=0.0083) and MARCH3 (P=0.049).

MARCH3 expression is correlated with immune in ltration level in HCC
To further analyze the effect of MARCH3 gene on HCC occurrence, we used TIMER to examine the expression lever of MARCH3 in liver cancer. The data showed that MARCH3 mRNA was highly expressed in normal tissues than tumor tissues (Fig. 6a). Tumor-in ltrating lymphocytes are an independent predictor of sentinel lymph node status and survival in cancers. Tumor purity is an important factor that in uences the analysis of immune in ltration in clinical tumor samples by genomic approaches, and TIMER and GEPIA have most of the homologous data from TCGA. By analyzing data from TIMER, we found that MARCH3 had a negative correlation with tumor purity in LIHC. And MARCH3 expression level also correlated with poorer prognosis and high in ltration in LIHC, as the results showed that MARCH3 had positive related with B cell (r=0.365, p=2.64e-12), CD8+ T cell (r=0.312, p=3.53e-9), CD4+ T cell (r=0.417, p=3.23e-21), macrophage (r=0.514, p=2.09e-24), neutrophil (r=0.437, p=1.72e-72) and dendritic cell (r=0.44, p=9.85e-16) (Fig. 6b). Next, we explored the effects of the copy number variation (CNV) of MARCH3 on the immune in ltration, the results showed that high amplication of MARCH3 could increase the immune in ltration of B cell (P 0.001) CD4+ T cell (P 0.01) macrophage (P 0.01), neutrophil (P 0.0001) and dendritic cell (P 0.01), while arm-lever gain of MARCH3 could increase the immune in ltration of CD4+ T cell (P 0.001) macrophage (P 0.01), neutrophil (P 0.01) (Fig. 6c). This nding strongly indicates MARCH3 plays a crucial role in immune in ltration in HCC.

Correlation analysis between MARCH3 expression and Immune marker sets
In order to analyze the correlation between MARCH3 and the diverse immune in ltrating cells, we further investigate the relationships between MARCH3 and immune marker sets of various immune cells of LIHC in the GEPIA and TIMER databases. We studied the associations between MARCH3 expression and immune marker genes of different immune cells, such as CD8+ T cells, T cells (general), B cells, monocytes, TAMs, M1 and M2 macrophages, neutrophils, NK cells and DCs in LIHC, and the different functional T cells included Th1 cells, Th2 cells, Tfh cells, Th17 cells, Tregs and exhausted T cells. After the correlation adjustment by purity, the results showed that MARCH3 expression was signi cantly correlated with most immune marker sets of various immune cells and different T cells in HCC (Table 2).
Interestingly, we also found that MARCH3 expression had strong positive correlations with the expression levels of most marker sets of T cell exhaustion (PD-1 and CTAL4) in HCC (Fig. 6d, P 0.0001). We also examined the correlation by analyzing a HCC TCGA data from GEPIA database, and our data showed that MARCH3 expression was positively associated with the marker sets of PD-1 (R=0.43; P 0.0001) and CTAL (R=0.42; P 0.0001) ( Table 3), whose results were consistent with TIMER analysis. These ndings suggest that MARCH3 may regulate the immune function to affect the HCC occurrence.

Discussion
In this study, to clarify the regulatory mechanism of IFN-α therapy on HCC, we analyzed the RNA expression pro le by RNA sequencing. According to comparing RNA pro le of patients with IFN-α or not, totally 60 (0.37%) mRNAs, 23 (0.08%) circRNAs and 23 (0.91%) miRNAs were differentially expressed. Then a ceRNA network included the circRNA-miRNA-mRNA and lncRNA-miRNA-mRNA were constructed.
Besides, the biological and functional attributes of these signi cantly differentially expressed RNAs were analyzed by GO and KEGG databases. Immune-related biological progresses were signi cantly more enriched. Survival analysis showed that three potential prognostic associated genes, namely FAM20A, IGFBP4 and MARCH3. Finally, we further analyzed the correlation between MARCH3 and immune in ltrate in HCC. To our knowledge, this is the rst publication to probe the ceRNA network of IFN-α therapy in HCC and identify the effect of MARCH3 associating with immune function on LIHC.
Of the genes which were involved in the ceRNA network, GO enrichment analysis indicated that these RNAs were mainly enriched in the synapse, external side of the plasma membrane, and juxtaparanode region of axon. IFN-α, a proin ammatory cytokine, which is commonly produced by natural killer cells and T-lymphocytes [19]. As under the normal conditions, these cells minimally entry into the central nervous system, IFN-α is generally not detected in the brain [20]. But several research have revealed that the high expression level of IFN-α can cause hippocampal development and abnormal cerebellar [21,22]. The effects of IFN-α on neurons is increasingly drawn attention. McHugh has published a comment that IFN-α can stimulate synapse loss [23]. And based on a mouse model, IFN-α inhibits the dendritic outgrowth which may lead to decreasing the rate of synapse formation [24]. Besides, IFN-α has also been proven to provide the enhanced axon protection [25]. And for localizing in the plasma membrane, this result has been widely accepted as many plasma membrane proteins, like tetherin and IFITM3, are identi ed in the interferon-treated cells [26]. In addition, the biological process of GO term and pathway have mapped the immune-related signaling. As we all know, the IFN-α is a family of cytokine mediators which involves in mediating the cellular immune system [27]. The normal interferon pathways always start with binding ligand to form a ternary complex, and then followed by activating the downstream signaling, like mediating the Janus tyrosine kinases [28]. Then the Toll-like receptors, tissue-destructive cytokines, in ammatory factors, cytokines are regulated to activate the STAT. Subsequently, the function of immune effector genes, like chemokines, phagocytic receptors, antiviral proteins, antigen-presenting molecules, are activated [29]. In current research, we discovered that the signi cantly differential expression of ceRNA mainly mapped on the immune-related pathway, and this result further proven the importance of the relationship between IFN-α and immune pathway.
FAM20A, encoding one of the "family with sequence similarity 20" (FAM20) proteins, function in the secretory signaling to promote protein phosphorylation. In human mesenchymal stem cells, previous paper had proven that this protein presented high expressed level following the IFN-α treatment [30]. Although this protein is originally identi ed from the hematopoietic cell, FAM20A presents a quite restricted high expression in liver tissue [31]. In addition, based on the data from the Human Protein Altas (HPA, https://www.proteinatlas.org/ENSG00000108950-FAM20A/pathology/liver+cancer#ihc), FAM20A displayed low expressed level in LIHC. And what' more, FAM20A is a prognostic factor in LIHC. However, another largely clinical cohorts are needed to further test the above results. IGFBP4 is a speci c insulinlike growths factor (IGF) binding protein which binds to IGF or not to tune the cellular activities, like cell proliferation and migration [32,33]. In vitro experiment has revealed that the IFN-α can increase the mRNA and protein level of IGFBP4 [34]. And this protein, functioning as a tumor suppressor, participates in driving epigenetic reprogramming of the hepatic carcinogenesis [35]. Although currently there is no research on the relationship between IGFBP4 and prognosis in liver cancer, it has been reported to be associated with prognosis in various cancer types, like lung cancer [36]. Our work rst revealed that patients with higher IGFBP4 expression had a longer survival time than that with lower IGFBP4 expression. Membrane-associated RING-CH-3 (MARCH3) is one of the membrane-associated MARCH family members which functions as a negative regulator of adaptive immunity [37]. The whole MARCH family includes 11 members, and 9 of them are transmembrane proteins. Many of these family members have been revealed to associate with prognosis. By researching on a mouse model, targeting MARCH1 exhibited signi cant inhibition of the growth of HCC [38], and MARCH8 has also played a crucial role in NSCLC against carcinogenesis and progression by validating in clinical samples [39]. In this work, we rst found that MARCH3 was highly expressed in HCC and patients with higher MARCH3 expression suffered shorter survive, which suggested MARCH3 as an unfavorable prognostic factor in HCC. Next, we found that MARCH3 expression was related with diverse immune in ltration levels in HCC. Our results demonstrated that there was a strong positive relationships between MARCH3 expression level and CD4 + T cells, macrophages, neutrophil and DCs, and signi cantly positive correlations between B cells, CD8 + T cells and MARCH3 expression in HCC. Moreover, the correlations between MARCH3 expression and the marker genes of immune cells implicated the role of MARCH3 in regulating tumor immunology in HCC.
We found that MARCH3 expression was highly correlated with the marker sets of T cell exhaustion, namely PD-1 and CTAL4. These correlations could be indicative of a potential mechanism where MARCH3 regulates T cell functions in HCC. Together these ndings suggest that the MARCH3 plays an important role in recruitment and regulation of immune in ltrating cells in HCC.
Although the current research constructed a ceRNA network in HCC for probing the regulation mechanism of IFN-α for the rst time, potential limitations existed. First, merely 8 HCC blood samples were obtained from the clinical patients which might lead to the selection bias. Second, in the current research, we haven't provided the validation data on the discovered results based on the downstream experiments. Anyway, a clinical cohort in larger patient sizes and a following validated experiment are unmet needed in the future.

Conclusions
We constructed a ceRNA network and analyzed the biological function of RNAs with signi cantly differential expression in HCC on researching the effect from IFN-α. The results indicated that immunerelated pathways played crucial role in participating in IFN-α treatment. Three genes (FAM20A, IGFBP4 and MARCH3) were identi ed as the prognostic markers, respectively. Finally, we identi ed MARCH3 as a vital factor in recruitment and regulation of immune in ltrating cells in HCC. These results laid the foundation on understanding the regulatory mechanism of IFN-α treatment.

Availability of data and materials
The data sets generated and analyzed during the current study are available from the corresponding author on reasonable request.