- 1Department of Gastroenterology, Affiliated Sixth People’s Hospital South Campus of Shanghai Jiaotong University, Shanghai, China
- 2Center for Bioinformatics and Computational Biology, The Institute of Biomedical Sciences, School of Life Sciences, East China Normal University, Shanghai, China
- 3Department of Gastroenterology, Affiliated Fengxian Hospital of Southern Medical University, Shanghai, China
- 4The First Affiliated Hospital, Zhejiang University, Hangzhou, China
- 5Shanghai East Hospital, Tongji University School of Medicine, Shanghai, China
Transcriptome profiling of hepatocellular carcinoma (HCC) by next-generation sequencing (NGS) technology has been broadly performed by previous studies, which facilitate our understanding of the molecular mechanisms of HCC formation, progression, and metastasis. However, few studies jointly analyze multiple types of noncoding RNAs (ncRNAs), including long noncoding RNAs (lncRNAs), circular RNAs (circRNAs), and micro-RNAs (miRNAs), and further uncover their implications in HCC. In this study, we observed that the circRNA cZRANB1 and lncRNA DUXAP10 were not only significantly upregulated in tumor tissues, but also higher expressed in blood exosomes of HCC as compared with healthy donors. From the analysis of subclass-associated dysregulated ncRNAs, we observed that DLX6-AS1, an antisense RNA of DLX6, and the sense gene DLX6 were highly expressed in S1, a subclass with a more invasive/disseminative phenotype. High correlation between DLX6-AS1 and DLX6 suggested that DLX6-AS1 may function via promoting the transcription of DLX6. Integrative analysis uncovers circRNA–miRNA, lncRNA–miRNA, and competing endogenous RNA networks (ceRNAs). Specifically, cZRANB1, LINC00501, CTD-2008L17.2, and SLC7A11-AS1 may function as ceRNAs that regulate mRNAs by competing the shared miRNAs. Further prognostic analysis demonstrated that the dysregulated ncRNAs had the potential to predict HCC patients’ overall survival. In summary, we identified some novel circRNAs and miRNAs, and dysregulated ncRNAs that could participate in HCC tumorigenesis and progression by inducing transcription of their neighboring genes, increasing their derived miRNAs, or acting as miRNA sponges. Moreover, our systematic analysis provides not only rich data resources for related researchers, but also new insights into the molecular basis of how different ncRNAs coordinately or antagonistically participate in the pathogenesis process of diseases.
Introduction
Hepatocellular carcinoma (HCC) is one of the most common malignancies and the second leading cause of cancer deaths worldwide (1). It is supposed that newly diagnosed HCC cases each year will be more than 1 million worldwide in 2025 (2, 3). Moreover, both morbidity and mortality of HCC in China are higher than other countries. Overall, the delayed diagnosis of a majority of HCC patients prohibits surgery or other effective radical treatments, resulting in poor prognosis, which makes the 5-year survival rate less than 15% (4). Currently, the widely used clinical biomarker for HCC diagnosis is alpha fetoprotein (AFP), while its sensitivity is only about 60% (5). The available targeted drug recommended by definitive guides in clinical practice (6) for advanced HCC is sorafenib, which, however, is limited in improving the overall survival (7, 8). Besides sorafenib, other targeted drugs, such as sunitinib (9), brivanib (10), and everolimus (11), were tested in clinical trials, but all failed in the third phase (4).
The latest developments in next-generation sequencing (NGS) technologies have profiled mutational spectrums, and deregulated expression and epigenetic changes of several cancers by TCGA studies (12). Whole genome or exome sequencing of HCC has identified some recurrently mutated genes, such as TERT promoter (54%–60%), p53 (12%–48%), β-catenin (11%–37%), and Axin (5%–15%) (13). Transcriptome sequencing, including RNA sequencing and small RNA sequencing, of cancers shows remarkable potential to identify both novel markers and uncharacterized aspects of tumor biology, particularly some non-coding RNAs. Studies about microRNAs (miRNAs) show that miRNAs are closely related to HCC tumorigenesis, development, and metastasis (14, 15). For example, the downregulation of miR-28-5p is related to HCC metastasis, recurrence, and poor prognosis (16). MiR-188-5p can inhibit the proliferation and metastasis of HCC by targeting FGF5 (17). Decreased expression of miRNA-122 has been frequently detected in hepatitis B virus (HBV)-related HCCs, and re-expressing of miRNA-122 conversely inhibits HCC progression (18). Moreover, long non-coding RNAs (lncRNAs), which are generally unable to encode proteins, are also involved in tumor formation, development, or metastasis. Suppression of lncRNA Hotair can significantly inhibit the expression of miR-218 and induce cell cycle arrest in G1 phase and thereby inhibit HCC progression (19). Overexpression of lncRNA HULC in liver cancer promotes HCC proliferation by downregulating tumor suppressor gene p18 (20). In addition to lncRNAs and miRNAs, circular RNAs (circRNAs), produced by non-canonical splicing events that join a splice donor to an upstream splice acceptor, are also abundant and conserved in both normal and tumor cells (21). Numerous studies show that circRNAs can act as key regulators in cancer by regulating transcription or post-transcription of driver genes (22–24).
To comprehensively profile multiple types of ncRNAs in HCC and uncover their roles in HCC tumorigenesis or progression, we performed RNA sequencing in parallel with small RNA sequencing of eight paired HCC and matched pare-cancerous tissues. From the systematic comparisons, we characterized tumorigenesis- and subclass-associated ncRNAs, and predicted their potential biological function and lncRNA/circRNA–miRNA–mRNA interactions (Figure 1), which highlighted some key functional ncRNAs at the post-transcriptional level.
Results
Transcriptome profiling of multiple types of RNAs in HCC and matched para-cancerous tissues
To comprehensively delineate the transcriptome of HCC and matched para-cancerous tissues, we performed small RNA sequencing and rRNA depletion-based total RNA sequencing on eight pairs of HCC and matched para-cancerous tissues, respectively. The workflow analyzing the small RNA and RNA sequencing datasets is illustrated in Figure 2A. We detected 581 mature miRNAs annotated by miRBase (25) (version 20) and 79 novel miRNAs annotated by miRDeep2 (26) with small RNA sequencing (Supplementary Table 1).
Figure 2 Transcriptome profiling of non-coding RNAs by RNA-seq in HCC. (A) The workflow for ncRNA detection and quantification. The circRNA and miRNA detection were performed by CIRI2 and miRDeep2, respectively. (B) The number of circRNAs originated from different genomic regions, including exonic, intronic, and intergenic regions. The novel circRNAs and miRNAs are displayed in pie charts of (C, D). The pink and blue parts represent the validated and undetermined circRNAs/miRNAs by the independent dataset, respectively.
RNA sequencing datasets were processed by the reference-based (GENCODE (27) version 19) gene expression quantification pipeline and circular RNA (circRNA) prediction pipeline, respectively. We identified a total of 25,197 genes, which consists of 18,038 protein-coding genes and 7,159 lncRNAs (read count >5 in more than four samples). In addition, RNA sequencing based on rRNA-depleted RNA library made it feasible to detect and quantify circRNAs. The circRNA prediction pipeline for the RNA sequencing data identified 13,308 circRNAs, 6,918 of which were annotated by circBase (28) (Supplementary Table 2). To annotate the origin of circRNAs, we determined whether the splice sites of circRNAs have been annotated according to the gene annotation from GENCODE. About 85.82%, 7.0%, and 7.18% of the circRNAs originated from exonic, intergenic, and intronic regions, respectively (Figure 2B). Moreover, we observed 52.81% of genes produced more than two circRNA isoforms, which indicated that alternative circularization extensively existed in HCC and para-cancerous tissues.
To further investigate whether the novel circRNAs and miRNAs could be identified in an independent HCC cohort, we performed identical data analysis on the publicly available RNA-seq and miRNA-seq datasets of 20 paired HCC and para-cancerous tissues (29). For the novel ncRNAs, including circRNAs and miRNAs, 4,871 circRNAs (76%) and 67 miRNAs (85%) were identified in the independent HCC cohort (Figures 2C, D). The high reproducibility demonstrated that the circRNA and miRNA prediction pipelines were robust, and the novel ncRNAs were highly reliable.
Identification and validation of dysregulated ncRNAs in HCC tumorigenesis
The dysregulated ncRNAs in HCC tumorigenesis should be aberrantly expressed between primary HCC and para-cancerous tissues. Prior to identifying the dysregulated ncRNAs in HCC tumorigenesis, the ncRNAs were quantified at count-based level to perform differential expression analysis. In total, we identified 844 lncRNAs, 80 circRNAs, and 114 miRNAs, which were differentially expressed between HCC and para-cancerous tissues (adjusted p-value <0.1 and fold-change >2 or < 1/2). The dysregulation of ~65% lncRNAs, ~76% circRNAs, and ~62% miRNAs was validated in the independent dataset (Figures 3A–C), including 20 novel circRNAs and one novel miRNA. Among the dysregulated ncRNAs validated by the independent dataset, PVT1, GAS5, DDX11-AS1, LINC-ROR, hsa-miR-483-5p, hsa-miR-139-5p, hsa-miR-150-5p, hsa-miR-195-5p, and hsa-miR-199a-5p were associated with HCC by previous studies (30–38).
Figure 3 Tumorigenesis-associated dysregulated ncRNAs in HCC. The expression patterns of lncRNAs, circRNAs, and miRNAs are displayed in (A-C). The red and blue points represent up- and downregulated ncRNAs in HCC. The solid points were also detected to be differentially expressed in the independent dataset. (D) The normalized expression levels of each of the top 10 ncRNAs by differential expression analysis, which are separated by three panels. HCC and para-cancerous samples were clustered by hierarchical clustering analysis, and are represented by the red and purple band on the top. The expression patterns of circRNA (hsa_circ_0000268) and lncRNA DUXAP10 in blood exosomes of HCC and healthy donors are displayed in (E, F).
We then conducted hierarchical clustering analysis to obtain the systematic comparison of ncRNA expression levels across different samples. Using each of the top 10 ncRNAs, HCC tissues were clustered into the same branch, and the para-cancerous tissues were clustered in the other branch, indicating that the ncRNAs could clearly distinguish cancerous from para-cancerous tissues (Figure 3D). To determine whether these ncRNAs were present in HCC blood exosomes, we investigated their expression patterns in the blood exosomes of 21 HCC and 32 healthy donors from the publicly available RNA-seq dataset (see Materials and Methods).
Notably, one upregulated circRNA hsa_circ_0000268 (cZRANB1) in HCC was also expressed in the exosome of HCC patients higher than the healthy donors (Figure 3E, p-value = 0.01). Moreover, the lncRNA DUXAP10 was also highly expressed in the exosome of HCC patients as compared with healthy controls (Figure 3F, p-value = 4.18e-5). In addition to potentially promoting HCC tumorigenesis, the two upregulated ncRNAs, especially the lncRNA DUXAP10, in the exosome of HCC patients may also be potential diagnostic biomarkers for HCC. The results suggested that the dysregulated ncRNAs in HCC may not only contribute to the tumorigenesis, but also have the potential to be diagnostic biomarkers for HCC.
Identification of dysregulated ncRNAs in HCC subclasses
Hoshida et al. (39) defined three robust HCC subclasses (termed S1, S2, and S3) based on a meta-analysis of gene expression profiles in datasets from eight independent patient cohorts. We classified 8 HCCs in this study and another 20 HCCs in the independent dataset into three subclasses (S1, S2, and S3) with 10, 6, and 12 cases based on the expression patterns of 619 signatures (Materials and Methods). We then performed differential expression analysis to identify dysregulated ncRNAs for each subclass (Figures 4A, B, adjusted p-value <0.25 and fold change >2 or < 1/2, Materials and Methods).
Figure 4 The subclass-associated dysregulated ncRNAs in HCC. The upregulated (A) and downregulated (B) ncRNAs in Hoshida subclasses. (C) Co-expression-based GSEA for DLX6-AS1. The highly correlated genes with DLX6-AS1 are significantly enriched in the gene set from regulation of nuclear beta-catenin signaling and target gene transcription. The expression patterns of MIR146A (D) and hsa-miR-146a-5p (E) in Hoshida subclasses. The miRNA host gene, MIR146A, and its derived miRNA, hsa-miR-146a-5p, are upregulated in the S1 subclass. (F) The hsa_circ_0007099 is downregulated in the S2 subclass.
As described by Hoshida et al. (39), S1 is a subclass with a more invasive/disseminative phenotype. We observed that DLX6-AS1, an antisense RNA of DLX6, was highly expressed in S1. A previous study (40) found that DLX6-AS1 was highly expressed in lung adenocarcinoma, and knockdown of DLX6-AS1 could significantly decrease the mRNA and protein expression of DLX6. Accordingly, high correlation (Pearson correlation coefficient, PCC > 0.8) between DLX6-AS1 and DLX6 was observed in HCC and para-cancerous tissues, suggesting that DLX6-AS1 may cis-regulate DLX6 transcription in HCC. Co-expression-based GSEA revealed that DLX6-AS1 was co-expressed with genes involved in the regulation of nuclear beta-catenin signaling and target gene transcription (Figure 4C), such as SALL4, ZCCHC12, and BCL9, suggesting that DLX6-AS1 may participate in beta-catenin-induced gene transcription in HCC.
Notably, the lncRNA MIR146A and its derived miRNA hsa-miR-146a-5p were significantly upregulated in the S1 subclass (Figures 4D, E), indicating that MIR146A may function via promoting the transcription of their derived miRNAs. To our knowledge, most of the subclass-related dysregulated circRNAs have not been reported previously. Exceptionally, hsa_circ_0007099, which was downregulated in S2, was reported to be downregulated in gastric cancer (41) (Figure 4F). As characterized by Hoshida et al. (39), MYC and PI3K/Akt activations were the features of subclass S2, which may be inversely associated with hsa_circ_0007099. In summary, the subclass-associated dysregulated ncRNAs may participate in some signaling pathways specifically activated or inactivated in a certain subclass, and were strongly associated with some clinical characteristics.
Identification of functional ncRNAs by constructing ceRNA network
From the systematic analysis above, we identified 1,598 dysregulated lncRNAs and 284 circRNAs; however, for most of them, their functionality remained unknown. To identify the functional ncRNAs responsible for HCC tumorigenesis or progression, we predicted the potential regulatory network involving lncRNAs, circRNAs, miRNAs, and mRNAs, namely, ceRNA (competing endogenous RNA) networks. As ceRNAs could regulate mRNAs by competing for the shared miRNAs, we firstly predicted the miRNA binding sites of the dysregulated ncRNAs using miRanda. Secondly, the predicted and experimentally validated mRNA–miRNA interactions were obtained from TargetScan (42) and miRTarBase (43) databases, respectively. With the threshold at −0.4 for the correlation coefficient between miRNA and target, we successfully predicted 152,881 lncRNA–miRNA pairs, 2,063 circRNA–miRNA pairs, and 8,056 mRNA–miRNA pairs, respectively.
To model the ceRNA network, we performed the hypergeometric test combined with co-expression to predict the ncRNA–miRNA–mRNA regulatory relationships. Finally, we predicted 219 lncRNA–mRNA pairs and 31 circRNA–mRNA pairs (Supplementary Table 3). Notably, the circRNA cZRANB1 (hsa_circ_0000268), which was identified as a miRNA sponge by a previous study (44), shared five miRNAs, including hsa-let-7a-5p, hsa-let-7c-5p, hsa-miR-26a-5p, hsa-miR-125b-5p, hsa-miR-195-5p, and hsa-miR-497-5p, with HMGA1 (Figure 5A). In particular, cZEANB1 was upregulated in both HCC tissues and exosomes. HMGA1 has been found to associate with tumor invasion in several cancers (45–48), which was also observed in HCC based on its significantly positive correlation with MMP9 (PCC > 0.6). The results indicated that cZRANB1 may promote tumor invasion by competing for the miRNAs with HMGA1, thereby increasing the expression level of HMGA1.
Figure 5 The ceRNA networks for dysregulated ncRNAs. The schematic diagrams for the cZRANB1, LINC00501/CTD-2008L17.2, and SLC7A11-AS1 are displayed in (A-C). The round, rhombus, rectangle, and rhomboid symbols represent the circRNAs, lncRNAs, mRNAs, and microRNAs, respectively. The symbols filled with red or blue color represent up- or downregulated in HCC.
For the lncRNA–mRNA pairs, we noticed that two lncRNAs, LINC00501 and CTD-2008L17.2, which were overexpressed in S2 subclass, may compete for miRNAs with LIN28B (Figure 5B). As described above, the S2 subclass was characterized by activation of MYC and the PI3k/Akt pathway, and the ceRNA network may be associated with these activations. In particular, the shared miRNAs belonged to the let-7 family, which function as tumor suppressor in several cancers (49, 50). Interestingly, LIN28B could also suppress the biogenesis of miRNAs including let-7 family miRNAs (51), indicating that LINC00501, CTD-2008L17.2, let-7 family miRNAs, and LIN28B constituted a complex regulatory network.
Furthermore, we also observed that SLC7A11-AS1 was predicted to regulate the driver gene CCNE1 by competing for the miRNAs, including hsa-miR-125b-5p, hsa-miR-195-5p, hsa-miR-497-5p, hsa-miR-4524a-5p, and hsa-miR-26b-5p (Figure 5C). It should be noted that the shared miRNAs were significantly downregulated, and were inversely correlated with both CCNE1 and SLC7A11-AS1. Based on the analyses, the predicted ceRNA networks provided some evidence about the post-transcriptionally regulatory roles of the ncRNAs in HCC.
Prognostic significance of the dysregulated lncRNAs
To our knowledge, there were no publicly available HCC circRNA expression datasets with clinical characteristics. Moreover, the associations between miRNAs and HCC prognosis have been extensively studied (52–55). Therefore, we only tested the prognostic significance for both tumorigenesis- and subclass-associated dysregulated lncRNAs based on two approaches. Firstly, we evaluated their associations with overall survival in the TCGA cohort (56). However, only 686 of 1,598 dysregulated lncRNAs were detected in the TCGA gene expression dataset due to different gene annotation versions, 68 of which were significantly associated with the HCC overall survival (p-value < 0.05), indicating that these lncRNAs had the potential to predict HCC overall survival. Secondly, to further test the prognostic significance of the dysregulated lncRNAs, especially the undetected lncRNAs in the TCGA cohort, we collected 86 HCC samples from four other independent RNA-seq datasets. Combined with the 28 samples used for ncRNA detection and quantification in our study, 69 of the 114 HCC samples were classified into poor or good prognosis subgroup based on the survival-associated genes by Lee et al. (57) (Materials and Methods, FDR < 0.05). Differential expression analysis was thus performed on the lncRNA expression. Finally, 498 dysregulated lncRNAs were identified as differentially expressed lncRNAs between HCC samples with poor and good prognosis (FDR < 0.05 and fold change >2 or < 1/2).
Notably, CTD-2008L17.2, which may compete for let-7 family miRNAs with LIN28B, was negatively associated with patient’s survival time in the univariate survival analysis of the TCGA dataset (Figure 6A). Twenty-nine lncRNAs were associated with the HCC prognosis based on the two approaches (Figure 6C), of which, AC068196.1, upregulated in S2 subclass, was the lncRNA most significantly associated with HCC survival (Figures 6B, C). In addition, 313 of the undetected lncRNAs in the TCGA dataset were also closely associated with HCC prognosis based on the comparison of their expression levels between HCC with poor and good prognosis. Interestingly, DUXAP10, the potential HCC diagnostic biomarker in HCC blood exosomes, was also overexpressed in HCC samples with poor prognosis (FDR < 0.05), suggesting that the expression of this lncRNA in exosome may be used to predict patient’s prognosis with validation from further clinical studies. The prognostic analysis based on the two approaches further demonstrated that the functional ncRNAs had the potential to predict patient’s survival.
Figure 6 Clinical associations between dysregulated ncRNAs and HCC prognosis. The Kaplan–Meier curves show significant overall survival between HCC samples with high and low expression of CTD-2008L17.2 (A) and AC068196.1 (B). (C) The 29 dysregulated lncRNAs that are significantly associated with HCC prognosis by both survival analysis of the TCGA dataset and differential expression analysis of HCC samples with good and poor prognosis from six HCC datasets.
Discussion
The molecular basis of HCC about protein-coding genes has largely been studied in the context of tumorigenesis, progression, and metastasis. Despite extensive research about the function of protein-coding genes in HCC, the lack of effective biomarkers for HCC diagnosis or prognostic prediction is still not thoroughly solved. Meanwhile, a majority of non-coding RNAs are characterized to act as cancer driver RNAs, and understanding their deregulation and regulatory roles can facilitate the development of new diagnostic or therapeutic strategies.
To our knowledge, this is the first study utilizing NGS data to simultaneously characterize the dysregulated mRNAs, lncRNA, circRNAs, and miRNAs in HCC and matched para-cancerous tissues. We successfully predicted some novel circRNAs and miRNAs, and quantified their expression levels based on RNA and small RNA sequencing data. The independent dataset validated about ~70% novel circRNAs and miRNAs, and more than 60% dysregulated tumorigenesis-associated ncRNAs, which demonstrated that the data analysis strategies were robust and the predicted novel ncRNAs were highly reliable. Remarkably, the expression levels of circRNA ZRANB1 and lncRNA DUXAP10 that were upregulated in HCC tissues were also observed higher in the blood exosomes of HCC than healthy donors. With the more comprehensive assessment by further experimental and clinical validation, the two ncRNAs may be used for HCC diagnosis.
As the HCC subclasses were well characterized by Hoshida et al., the subclass-associated ncRNAs could be directly linked to HCC clinical characteristics and the activated or inactivated pathways. With the well-characterized clinical phenotypes and dysregulated pathways for the HCC subclasses, we could easily speculate the impact of the dysregulated ncRNAs on HCC tumorigenesis or progression.
Nowadays, more and more studies uncovered the biological function or pathways that some ncRNAs may participate in; however, for most of them, their function still remained unknown. To further uncover the potential regulatory ncRNAs at the post-transcriptional level, we built ncRNA–miRNA–mRNA interaction networks. Non-coding RNAs, including cZRANB1, LINC00501, CTD-2008L17.2, AC092171.4, and SLC7A11-AS1, were identified as functional ncRNAs that competed for the shared miRNAs with mRNAs, and were predicted to be involved in HCC tumorigenesis or progression. In particular, cZRANB1 was upregulated in both HCC tissues and exosomes. The prognostic association analysis highlighted that CTD-2008L17.2, a potential miRNA sponge, and DUXAP10, a potential diagnostic biomarker in HCC exosome, were closely associated with HCC prognosis.
Even though this preliminary study provides an abundant data resource of non-coding RNAs in HCC, the lack of further experimental verification is the major limitation of this study. The systematic analysis revealed that the ncRNAs could function as oncogenic or tumor-suppressive RNAs by regulating gene transcription, post-transcription processes, or other mechanisms, some of which may be used as HCC diagnostic or prognostic biomarkers. Although further characterizing the molecular mechanism of the ncRNAs is extremely essential, our data analysis still provided the hint of how the dysregulated ncRNAs were involved in HCC tumorigenesis or progression for further studies, and novel insights into cancer biology.
Conclusions
In this study, we identified some novel circRNAs and miRNAs that were validated in the independent dataset. Systematic analysis identified dysregulated ncRNAs that participate in HCC tumorigenesis and progression by inducing transcription of their neighboring genes, increasing their derived miRNAs, or acting as miRNA sponges. Remarkably, two upregulated ncRNAs in the blood exosomes of HCC may be potential diagnostic biomarkers. The ncRNAs competing for miRNAs with mRNAs may be key regulators in HCC, which improved our understanding of the mechanisms about HCC tumorigenesis or progression. Moreover, the prognostic association analysis revealed that the dysregulated ncRNAs with key regulatory roles may also have the potential to predict HCC patients’ prognosis. In summary, our systematic analysis provides not only rich data resources for related researchers, but also new insights into the molecular basis of how different ncRNAs coordinately or antagonistically participate in the pathogenesis process of diseases.
Materials and methods
Patient samples
Eight pairs of cancerous and para-cancerous tissues of HCC were obtained from the Department of Hepatopancreatobiliary Surgery (in 2013), the First Affiliated Hospital of Zhejiang University (Hangzhou, China), and were frozen at −70°C.The diagnosis of HCC was confirmed by postsurgery pathology. This project was approved by the Human Research Ethics Committee of the First Affiliated Hospital of Zhejiang University. Prior to surgery, the patients were not treated with anticancer therapy.
All of the eight patients had a history of chronic hepatitis B (CHB), and a single tumor in the liver, without extrahepatic metastasis. Three patients were complicated with alcoholic liver diseases (ALDs), and another three patients had portal vein tumor thrombosis (PVTT) (Table 1).
All procedures in the experiment were in agreement with the ethical principles. Patients involved in this study have been informed and signed consent before collection of samples, and the data (which do not involve personal information) related to patients were completely anonymous.
RNA extraction and sequencing
Total RNAs from HCC and para-cancerous tissues were extracted using TRIzol following the manufacturer’s protocol. For the preparation of RNA-Seq libraries, total RNA was treated with the Ribo-Zero rRNA Removal Kit to remove rRNA according to the manufacturer’s instructions. Standard protocols recommended for small RNA-seq (Illumina TruSeq Small RNA) were used. RNA and small RNA sequencing libraries for the Illumina Hiseq 2000 platform were constructed according to the manufacturer’s instructions (Illumina).
RNA-seq read mapping, ncRNA detection, and quantification
One hundred-base-pair paired-end reads were mapped to UCSC human reference genome (GRCH37/hg19) using HISAT2 (58) version 2.0.5 with gene models from GENCODE v19 and default options. The SAM format files were compressed and sorted by samtools (59) version 1.3.1. transcripts were quantified by StringTie (60) v1.2.4 at count-based levels. The circRNA was predicted and quantified by CIRI2 (61) based on the alignment by BWA (62). The small RNA-seq data were preprocessed by Trimmomatic (63), and the miRNA prediction and quantification were implemented by miRDeep2 (26). The genes or ncRNAs, which are expressed (read count > 5) in at least 20% of samples, were retained for further analysis. The count-based expression matrix was used for differential expression analysis. The count-based expression values were normalized by a regularized log transformation in DESeq2 (64).
Differential expression analysis
The count-based expression levels of lncRNAs, mRNAs, circRNAs, and miRNAs were respectively analyzed by DESeq2 (64), which performed differential expression analysis based on negative binomial distribution.
Pathway enrichment analysis
The enrichment analysis was implemented in GSEA Java software (65). Pathway database was downloaded from MsigDB (http://software.broadinstitute.org/gsea/index.jsp). We only retained canonical pathways curated by KEGG (66) and NCI-PID (67) for enrichment analysis.
Cox regression-based survival analysis
The survival analysis based on the Cox regression model was implemented in R with packages survival and survcomp. The comparison of survival curves for high- and low-expression groups was performed using the G-rho family of tests with survdiff function.
Collection of independent datasets and identification of prognostic lncRNAs
We collected 106 HCC samples from five independent RNA-seq datasets, with accession numbers SRP039694 (68), SRP050551 (69), SRP062885 (70), SRP068976 (71), and SRP069212 (29) from the SRA (Sequence Read Archive, https://www.ncbi.nlm.nih.gov/sra) database. In particular, the RNA library of the SRP069212 dataset was constructed with the rRNA depletion protocol and small RNA data with accession number SRP068498, which was used to validate the detection and dysregulation of lncRNAs, circRNAs, and miRNAs. In addition, the accession numbers for RNA-seq data of the blood exosomes from HCC and healthy donors were SRP109668 and SRP109666. We adopted the same method to analyze the RNA-seq data of the blood exosomes and detect ncRNAs in blood exosomes.
To identify prognostic lncRNAs, the 106 HCC samples were classified into poor and good prognostic subgroups by the NTP algorithm (72). The survival signatures were obtained from the previous study by Lee et al. (57). The differential expression analysis was conducted between HCC samples with poor and good prognosis to identify the prognostic lncRNAs.
MiRNA–target prediction
MiRNA binding sites of lncRNAs and circRNAs are predicted by miRanda (73) with strict mode, that is, miRNA–target demands strict 5’ seed pairing. Moreover, the miRNAs and targets should be dysregulated with opposite directions; for example, if the miRNA is upregulated in tumor, its target gene should be downregulated in tumor, and vice versa. In addition, the interactions of mRNAs and miRNAs were downloaded from the miRTarBase database (43) (http://mirtarbase.mbc.nctu.edu.tw/php/index.php), which curated the experimentally validated microRNA–target interactions, and TargetScan (42) databases. A ceRNA network was constructed based on the hypergeometric test (p-value < 1e-6), the correlation between lncRNA/circRNA and mRNAs (>0.6), and the number of shared miRNAs (>4).
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Author contributions
In this study, JC conducted the experiment design, XX collected the samples, TS and LZ designed the data analysis, and LZ performed the data analysis. LZ and CW drafted the manuscript. TS, JC, XX, and XL revised and finalized the manuscript. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by the National High Technology Research and Development Program of China (2012AA020204 and 2015AA020108), the China Human Proteomics Project (2014DFB30010 and 2014DFB30030), the National Science Foundation of China (31671377), the Natural Science Foundation of Shanghai (No. 16ZR1429100), and Shanghai 111 Project (B14019).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.927524/full#supplementary-material
References
1. Torre LA, Bray F, Siegel RL, Ferlay J, Lortet-Tieulent J, Jemal A. Global cancer statistics, 2012. CA: Cancer J Clin (2015) 65(2):87–108. doi: 10.3322/caac.21262
2. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, et al. Cancer incidence and mortality worldwide: Sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer (2015) 136(5):E359–386. doi: 10.1002/ijc.29210
3. Simard EP, Ward EM, Siegel R, Jemal A. Cancers with increasing incidence trends in the united states: 1999 through 2008. CA: Cancer J Clin (2012) 62(2):118–28. doi: 10.3322/caac.20141
4. Dhanasekaran R, Venkatesh SK, Torbenson MS, Roberts LR. Clinical implications of basic research in hepatocellular carcinoma. J Hepatol (2016) 64(3):736–45. doi: 10.1016/j.jhep.2015.09.008
5. Forner A, Llovet JM, Bruix J. Hepatocellular carcinoma. Lancet (London England) (2012) 379(9822):1245–55. doi: 10.1016/S0140-6736(11)61347-0
6. EASL-EORTC clinical practice guidelines: Management of hepatocellular carcinoma. J Hepatol (2012) 56(4):908–43. doi: 10.1016/j.jhep.2011.12.001
7. Llovet JMRS, Mazzaferro V, Hilgard P, Gane E, Blanc JF, et al. Sorafenib in advanced hepatocellular carcinoma. N Engl J Med (2008) 359(((4):378–90. doi: 10.1056/NEJMoa0708857
8. Bruix J, Raoul JL, Sherman M, Mazzaferro V, Bolondi L, Craxi A, et al. Efficacy and safety of sorafenib in patients with advanced hepatocellular carcinoma: Subanalyses of a phase III trial. J Hepatol Oct (2012) 57(4):821–9. doi: 10.1016/j.jhep.2012.06.014
9. Cheng AL, Kang YK, Lin DY, Park JW, Kudo M, Qin S, et al. Sunitinib versus sorafenib in advanced hepatocellular cancer: results of a randomized phase III trial. J Clin Oncol Off J Am Soc Clin Oncol (2013) 31(32):4067–75. doi: 10.1200/JCO.2012.45.8372
10. Llovet JM, Decaens T, Raoul JL, Boucher E, Kudo M, Chang C, et al. Brivanib in patients with advanced hepatocellular carcinoma who were intolerant to sorafenib or for whom sorafenib failed: Results from the randomized phase III BRISK-PS study. J Clin Oncol Off J Am Soc Clin Oncol (2013) 31(28):3509–16. doi: 10.1200/JCO.2012.47.3009
11. Zhu AX, Kudo M, Assenat E, Cattan S, Kang YK, Lim HY, et al. Effect of everolimus on survival in advanced hepatocellular carcinoma after failure of sorafenib: The EVOLVE-1 randomized clinical trial. JAMA (2014) 312(1):57–67. doi: 10.1001/jama.2014.7189
12. Cancer Genome Atlas Research N, Weinstein JN, Collisson EA, Mills GB, Shaw KM, Ozenberger BA, et al. The cancer genome atlas pan-cancer analysis project. Nat Genet (2013) 45(10):1113–20. doi: 10.1038/ng.2764.
13. Zucman-Rossi J, Villanueva A, Nault JC, Llovet JM. Genetic landscape and biomarkers of hepatocellular carcinoma. Gastroenterology (2015) 149(5):1226–39.e1224. doi: 10.1053/j.gastro.2015.05.061
14. Wojcicka A, Swierniak M, Kornasiewicz O, Gierlikowski W, Maciag M, Kolanowska M, et al. Next generation sequencing reveals microRNA isoforms in liver cirrhosis and hepatocellular carcinoma. Int J Biochem Cell Biol (2014) 53:208–17. doi: 10.1016/j.biocel.2014.05.020
15. Selitsky SR, Dinh TA, Toth CL, Kurtz CL, Honda M, Struck BR, et al. Transcriptomic analysis of chronic hepatitis b and c and liver cancer reveals MicroRNA-mediated control of cholesterol synthesis programs. mBio (2015) 6(6):e01500–01515. doi: 10.1128/mBio.01500-15
16. Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell (2009) 136(2):215–33. doi: 10.1016/j.cell.2009.01.002
17. Fang F, Chang RM, Yu L, Lei X, Xiao S, Yang H, et al. MicroRNA-188-5p suppresses tumor cell proliferation and metastasis by directly targeting FGF5 in hepatocellular carcinoma. J Hepatol (2015) 63(4):874–85. doi: 10.1016/j.jhep.2015.05.008
18. Bandiera S, Pfeffer S, Baumert TF, Zeisel MB. miR-122–a key factor and therapeutic target in liver disease. J Hepatol (2015) 62(2):448–57. doi: 10.1016/j.jhep.2014.10.004
19. Fu WM, Zhu X, Wang WM, Lu YF, Hu BG, Wang H, et al. Hotair mediates hepatocarcinogenesis through suppressing miRNA-218 expression and activating P14 and P16 signaling. J Hepatol (2015) 63(4):886–95. doi: 10.1016/j.jhep.2015.05.016
20. Du Y, Kong G, You X, Zhang S, Zhang T, Gao Y, et al. Elevation of highly up-regulated in liver cancer (HULC) by hepatitis b virus X protein promotes hepatoma cell proliferation via down-regulating p18. J Biol Chem (2012) 287(31):26302–11. doi: 10.1074/jbc.M112.342113
21. Salzman J, Gawad C, Wang PL, Lacayo N, Brown PO. Circular RNAs are the predominant transcript isoform from hundreds of human genes in diverse cell types. PLoS One (2012) 7(2):e30733. doi: 10.1371/journal.pone.0030733
22. Hou LD, Zhang J. Circular RNAs: An emerging type of RNA in cancer. Int J Immunopathology Pharmacol (2017) 30(1):1–6. doi: 10.1177/0394632016686985
23. Wang F, Nazarali AJ, Ji S. Circular RNAs as potential biomarkers for cancer diagnosis and therapy. Am J Cancer Res (2016) 6(6):1167–76.
24. Zheng Q, Bao C, Guo W, Li S, Chen J, Chen B, et al. Circular RNA profiling reveals an abundant circHIPK3 that regulates cell growth by sponging multiple miRNAs. Nat Commun (2016) 7:11215. doi: 10.1038/ncomms11215
25. Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res (2014) ;42(Database issue):D68–73. doi: 10.1093/nar/gkt1181
26. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res (2012) 40(1):37–52. doi: 10.1093/nar/gkr688
27. Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, et al. GENCODE: The reference human genome annotation for the ENCODE project. Genome Res (2012) 22(9):1760–74. doi: 10.1101/gr.135350.111
28. Glazar P, Papavasileiou P, Rajewsky N. circBase: A database for circular RNAs. Rna. (2014) 20(11):1666–70. doi: 10.1261/rna.043687.113
29. Yang Y, Chen L, Gu J, Zhang H, Yuan J, Lian Q, et al. Recurrently deregulated lncRNAs in hepatocellular carcinoma. Nat Commun (2017) 8:14421. doi: 10.1038/ncomms14421
30. Wang F, Yuan JH, Wang SB, Yang F, Yuan SX, Ye C, et al. Oncofetal long noncoding RNA PVT1 promotes proliferation and stem cell-like property of hepatocellular carcinoma cells by stabilizing NOP2. Hepatology (2014) 60(4):1278–90. doi: 10.1002/hep.27239
31. Tu ZQ, Li RJ, Mei JZ, Li XH. Down-regulation of long non-coding RNA GAS5 is associated with the prognosis of hepatocellular carcinoma. Int J Clin Exp Pathol (2014) 7(7):4303–9.
32. Liao HT, Huang JW, Lan T, Wang JJ, Zhu B, Yuan KF, et al. Identification of the aberrantly expressed LncRNAs in hepatocellular carcinoma: A bioinformatics analysis based on RNA-sequencing. Sci Rep (2018) 8(1):5395. doi: 10.1038/s41598-018-23647-1
33. Takahashi K, Yan IK, Kogure T, Haga H, Patel T. Extracellular vesicle-mediated transfer of long non-coding RNA ROR modulates chemosensitivity in human hepatocellular cancer. FEBS Open Bio (2014) 4:458–67. doi: 10.1016/j.fob.2014.04.007
34. Shen J, Wang A, Wang Q, Gurvich I, Siegel AB, Remotti H, et al. Exploration of genome-wide circulating microRNA in hepatocellular carcinoma: MiR-483-5p as a potential biomarker. Cancer Epidemiology Biomarkers Prev Publ Am Assoc Cancer Research cosponsored by Am Soc Prev Oncol (2013) 22(12):2364–73. doi: 10.1158/1055-9965.EPI-13-0237
35. Hua S, Lei L, Deng L, Weng X, Liu C, Qi X, et al. miR-139-5p inhibits aerobic glycolysis, cell proliferation, migration, and invasion in hepatocellular carcinoma via a reciprocal regulatory interaction with ETS1. Oncogene (2018) 37(12):1624–36. doi: 10.1038/s41388-017-0057-3
36. Zhang J, Luo N, Luo Y, Peng Z, Zhang T, Li S. microRNA-150 inhibits human CD133-positive liver cancer stem cells through negative regulation of the transcription factor c-myb. Int J Oncol (2012) 40(3):747–56. doi: 10.3892/ijo.2011.1242
37. Xu H, Hu YW, Zhao JY, Hu XM, Li SF, Wang YC, et al. MicroRNA-195-5p acts as an anti-oncogene by targeting PHF19 in hepatocellular carcinoma. Oncol Rep (2015) 34(1):175–82. doi: 10.3892/or.2015.3957
38. Hou J, Lin L, Zhou W, Wang Z, Ding G, Dong Q, et al. Identification of miRNomes in human liver and hepatocellular carcinoma reveals miR-199a/b-3p as therapeutic target for hepatocellular carcinoma. Cancer Cell (2011) 19(2):232–43. doi: 10.1016/j.ccr.2011.01.001
39. Hoshida Y, Nijman SM, Kobayashi M, Chan JA, Brunet JP, Chiang DY, et al. Integrative transcriptome analysis reveals common molecular subclasses of human hepatocellular carcinoma. Cancer Res (2009) 69(18):7385–92. doi: 10.1158/0008-5472.CAN-09-1089
40. Li J, Li P, Zhao W, Yang R, Chen S, Bai Y, et al. Expression of long non-coding RNA DLX6-AS1 in lung adenocarcinoma. Cancer Cell Int (2015) 15:48. doi: 10.1186/s12935-015-0201-5
41. Shao Y, Li J, Lu R, Li T, Yang Y, Xiao B, et al. Global circular RNA expression profile of human gastric cancer and its clinical significance. Cancer Med (2017) 6(6):1173–80. doi: 10.1002/cam4.1055
42. Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. eLife (2015) 4. doi: 10.7554/eLife.05005
43. Chou CH, Shrestha S, Yang CD, Chang NW, Lin YL, Liao KW, et al. miRTarBase update 2018: A resource for experimentally validated microRNA-target interactions. Nucleic Acids Res (2018) 46(D1):D296–302. doi: 10.1093/nar/gkx1067
44. Wang JJ, Shan K, Liu BH, Liu C, Zhou RM, Li XM, et al. Targeting circular RNA-ZRANB1 for therapeutic intervention in retinal neurodegeneration. Cell Death Dis (2018) 9(5):540. doi: 10.1038/s41419-018-0597-7
45. Liau SS, Jazag A, Whang EE. HMGA1 is a determinant of cellular invasiveness and in vivo metastatic potential in pancreatic adenocarcinoma. Cancer Res (2006) 66(24):11613–22. doi: 10.1158/0008-5472.CAN-06-1460
46. Zhong J, Liu C, Zhang QH, Chen L, Shen YY, Chen YJ, et al. TGF-beta1 induces HMGA1 expression: The role of HMGA1 in thyroid cancer proliferation and invasion. Int J Oncol (2017) 50(5):1567–78. doi: 10.3892/ijo.2017.3958
47. Zhong J, Liu C, Chen YJ, Zhang QH, Yang J, Kang X, et al. The association between S100A13 and HMGA1 in the modulation of thyroid cancer proliferation and invasion. J Trans Med (2016) 14:80. doi: 10.1186/s12967-016-0824-x
48. Liang L, Li X, Zhang X, Lv Z, He G, Zhao W, et al. MicroRNA-137, an HMGA1 target, suppresses colorectal cancer cell invasion and metastasis in mice by directly targeting FMNL2. Gastroenterology (2013) 144(3):624–35.e624. doi: 10.1053/j.gastro.2012.11.033
49. Shell S, Park SM, Radjabi AR, Schickel R, Kistner EO, Jewell DA, et al. Let-7 expression defines two differentiation stages of cancer. Proc Natl Acad Sci USA (2007) 104(27):11400–5. doi: 10.1073/pnas.0704372104
50. Boyerinas B, Park SM, Hau A, Murmann AE, Peter ME. The role of let-7 in cell differentiation and cancer. Endocrine-Related Cancer (2010) 17(1):F19–36. doi: 10.1677/ERC-09-0184
51. Balzeau J, Menezes MR, Cao S, Hagan JP. The LIN28/let-7 pathway in cancer. Front Genet (2017) 8:31. doi: 10.3389/fgene.2017.00031
52. Ji J, Shi J, Budhu A, Yu Z, Forgues M, Roessler S, et al. MicroRNA expression, survival, and response to interferon in liver cancer. N Engl J Med (2009) 361(15):1437–47. doi: 10.1056/NEJMoa0901282
53. Xu J, Li J, Zheng TH, Bai L, Liu ZJ. MicroRNAs in the occurrence and development of primary hepatocellular carcinoma. Adv Clin Exp Med Off Organ Wroclaw Med Univ (2016) 25(5):971–5. doi: 10.17219/acem/36460
54. Fiorino S, Bacchi-Reggiani ML, Visani M, Acquaviva G, Fornelli A, Masetti M, et al. MicroRNAs as possible biomarkers for diagnosis and prognosis of hepatitis b- and c-related-hepatocellular-carcinoma. World J Gastroenterol (2016) 22(15):3907–36. doi: 10.3748/wjg.v22.i15.3907
55. Anwar SL, Lehmann U. MicroRNAs: Emerging novel clinical biomarkers for hepatocellular carcinomas. J Clin Med (2015) 4(8):1631–50. doi: 10.3390/jcm4081631
56. Cancer Genome Atlas Research Network, Electronic address wbe, Cancer Genome Atlas Research N. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell (2017) 169(7):1327–41.e1323. doi: 10.1016/j.cell.2017.05.046
57. Lee JS, Chu IS, Heo J, Calvisi DF, Sun Z, Roskams T, et al. Classification and prediction of survival in hepatocellular carcinoma by gene expression profiling. Hepatology (2004) 40(3):667–76. doi: 10.1002/hep.20375
58. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods (2015) 12(4):357–60. doi: 10.1038/nmeth.3317
59. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence Alignment/Map format and SAMtools. Bioinformatics (2009) 25(16):2078–9. doi: 10.1093/bioinformatics/btp352
60. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol (2015) 33(3):290–5. doi: 10.1038/nbt.3122
61. Gao Y, Zhang J, Zhao F. Circular RNA identification based on multiple seed matching. Briefings Bioinf (2017) 19(5):803–10. doi: 10.1093/bib/bbx014
62. Li H, Durbin R. Fast and accurate long-read alignment with burrows-wheeler transform. Bioinformatics (2010) 26(5):589–95. doi: 10.1093/bioinformatics/btp698
63. Bolger AM, Lohse M, Usadel B. Trimmomatic: A flexible trimmer for illumina sequence data. Bioinformatics (2014) 30(15):2114–20. doi: 10.1093/bioinformatics/btu170
64. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol (2014) 15(12):550. doi: 10.1186/s13059-014-0550-8
65. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–50. doi: 10.1073/pnas.0506580102
66. Kanehisa M, Furumichi M, Tanabe M, Sato Y, Morishima K. KEGG: New perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res (2017) 45(D1):D353–61. doi: 10.1093/nar/gkw1092
67. Schaefer CF, Anthony K, Krupa S, Buchoff J, Day M, Hannay T, et al. PID: the pathway interaction database. Nucleic Acids Res (2009) 37(Database issue):D674–679. doi: 10.1093/nar/gkn653
68. Gao F, Liang H, Lu H, Wang J, Via M, Yuan Z, et al. Global analysis of DNA methylation in hepatocellular carcinoma by a liquid hybridization capture-based bisulfite sequencing approach. Clin Epigenet (2015) 7:86. doi: 10.1186/s13148-015-0121-1
69. Zhang H, Weng X, Ye J, He L, Zhou D, Liu Y. Promoter hypermethylation of TERT is associated with hepatocellular carcinoma in the han Chinese population. Clinics Res Hepatol Gastroenterol (2015) 39(5):600–9. doi: 10.1016/j.clinre.2015.01.002
70. Jiang Z, Jhunjhunwala S, Liu J, Haverty PM, Kennemer MI, Guan Y, et al. The effects of hepatitis b virus integration into the genomes of hepatocellular carcinoma patients. Genome Res (2012) 22(4):593–601. doi: 10.1101/gr.133926.111
71. Liu G, Hou G, Li L, Li Y, Zhou W, Liu L. Potential diagnostic and prognostic marker dimethylglycine dehydrogenase (DMGDH) suppresses hepatocellular carcinoma metastasis in vitro and in vivo. Oncotarget (2016) 7(22):32607–16. doi: 10.18632/oncotarget.8927
72. Hoshida Y. Nearest template prediction: A single-sample-based flexible class prediction with confidence assessment. PLoS One (2010) 5(11):e15543. doi: 10.1371/journal.pone.0015543
Keywords: transcriptome profiling, noncoding RNAs, competing endogenous RNA networks, miRNA sponge, tumorigenesis
Citation: Zhang L, Wang C, Lu X, Xu X, Shi T and Chen J (2022) Transcriptome sequencing of hepatocellular carcinoma uncovers multiple types of dysregulated ncRNAs. Front. Oncol. 12:927524. doi: 10.3389/fonc.2022.927524
Received: 24 April 2022; Accepted: 05 August 2022;
Published: 05 September 2022.
Edited by:
Chunquan Li, University of South China, ChinaReviewed by:
Beenish Rahat, Eunice Kennedy Shriver National Institute of Child Health and Human Development (NIH), United StatesYuchen Liu, Shenzhen University, China
Copyright © 2022 Zhang, Wang, Lu, Xu, Shi and Chen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jinlian Chen, d3FxXzAyMTAwMkAxNjMuY29t; Tieliu Shi, dGxzaGlAYmlvLmVjbnUuZWR1LmNu; Xiao Xu, YWEwMjAyMDR4dUAxNjMuY29t
†These authors have contributed equally to this work