《1. Introduction》

1. Introduction

Hepatocellular carcinoma (HCC) is one of the most common malignant tumors and has poor prognosis [1]. The incidence of HCC is increasing more rapidly than any other cancers, which deserves a great deal of attention [2]. The mortality of HCC has increased in both women and men during the last decades. Surgery is a crucial treatment for HCC [3]. Despite the combination of targeted drug therapy [4], local ablation [5], transcatheter arterial chemoembolization [6], and even liver transplantation [7], the prognosis of HCC remains grave. Thus, it is mandatory to attempt to better understand the pathogenesis of HCC and explore its possible therapeutic targets.

MicroRNAs (miRNAs), a type of small non-coding RNA with a length of 18–25 nucleotides, play pivotal roles in regulating tumorigenesis and tumor progression via targeting tumor-related pathways [8]. Currently, miRNA clusters, which are defined as a set of precursor miRNAs (pre-miRNAs) transcribed from the same orientation or an adjacent genome [9], are particularly attractive in cancer regulation. Their regulation and expression are characterized by a high degree of consistency [10]. In addition, they can target the same genes or pathways and achieve common effects due to the same or similar seed regions [11,12]. The chromosome 19 microRNA cluster (C19MC), as the largest miRNA cluster, contains 46 pre-miRNAs [13]. In our previous study, by mining The Cancer Genome Atlas (TCGA)^↑ miRNA profiles, we demonstrated that among the top 50 and top 100 differentially upregulated miRNAs in HCC tissues, 37 (80.4%) and 45 (97.8%), respectively, belong to C19MC [14]. These results are consistent with those of other studies [15–17] and suggest that C19MC could be a novel prognostic biomarker in HCC. However, there is a lack of functional studies exploring the mechanism underlying the effects of C19MC members on HCC.

^↑  https://portal.gdc.cancer.gov/

It is possible for miRNA to be transported among different cell types, including cancer cells and adjacent cells, which can affect the biology of recipient cells [18]. Exosomes or extracellular vesicles are the most frequently studied miRNA carriers during intercellular communication and contribute to the transfer of malignant tumor behaviors [19]. Some reports have suggested that C19MC is transferred via exosomes in other areas. The exosomalmediated transfer of C19MC in human placental trophoblasts for viral resistance has received widespread interest [20,21]. Bullerdiek et al. [22] also presented the hypothesis that exosome-derived C19MC might act as a target of treatment against tumorigenesis. However, experimental evidence that tumor cells are regulated by the transfer of C19MC via exosomes remains obscure. Thus, we aim to explore the biological function transfer of the C19MC in HCC via exosomes.

In this study, we found that four oncogenic pre-miRNAs of C19MC (i.e., mir-516a-1, mir-516a-2, mir-516b-1, and mir-516b2) co-spliced an identical mature miRNA named miR-516a-3p. We assessed the clinical value of miR-516a-3p in our cohort of 82 HCC patients. We evaluated the effects of miR-516a-3p on the malignant behaviors of HCC both in vitro and in vivo and examined the intercellular transportation of miR-516a-3p by studying exosomes. We also performed multi-omics (i.e., transcriptomics, proteomics, and metabolomics) analysis to systematically explore the molecular mechanism of miR-516a-3p-mediated oncogenicity.

《2. Materials and methods》

2. Materials and methods

《2.1. TCGA data analysis》

2.1. TCGA data analysis

TCGA data were analyzed as described previously [14]. In brief, the miRNA sequence and clinical data of HCC patients and the nontumor controls were acquired from the TCGA database. The edgR package was used to normalize the miRNA sequencing and screen the differential miRNAs. According to the cut-off value of each miRNA that was calculated with the receiver operating characteristics (ROC) curve, the expression levels of the miRNAs were divided into high and low expression groups. Kaplan–Meier analysis with a log-rank test was performed to detect the correlation of miRNA expression and HCC overall survival (OS). A univariate Cox proportional hazards regression analysis was used to screen the risk factors of HCC OS. The chi-square test was performed to analyze the clinicopathologic features of HCC with the expression of the miRNAs.

《2.2. Metabolomics》

2.2. Metabolomics

Metabolomics was performed as described previously [23]. For metabolites extraction, 5 ×10⁶ (per sample) miR-516a-3ptransfected SNU-449 (n = 8) and their controls were used for the metabolomics analysis. The metabolites of all samples were extracted and transferred for high-performance liquid chromatography with tandem mass spectrometry (HPLC–MS/MS) analysis, according to the manufacturer’s instructions. For metabolite identification, the raw mass spectrometry (MS) data was transformed into mzXML by using MSconvert software (ProteoWizard, USA) and was extracted with the XCMS package from software. According to the retention time (RT) and mass-to-charge ratio (m/z) data of each ion, the MS1 spectra were identified by means of the metaX package from software. In particular, the MS2 spectra were matched with the in-house standard spectrum library. The metabolites of the MS2 spectra were then identified and annotated with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and the Human Metabolome Database (HMDB). Pathway enrichment of the identified metabolites was performed with MetaboAnalyst 4.0 [24].

 http://www.metaboanalyst.ca

《2.3. Proteomics》

2.3. Proteomics

Proteomics analysis was performed for 5 ×10⁶ (per sample) miR-516a-3p-transfected SNU-449 (n = 3) and their controls. In brief, the protein was extracted from the samples with lysis buffer and was then boiled and centrifuged at 15 000 g for 15 min. A bicinchoninic acid (BCA) protein assay (Beyotime, China) was performed to detect the protein concentration of all the samples. After 20 μg protein per sample was resolved on 12% sodium dodecyl sulphate–polyacrylamide gel electrophoresis (SDS-PAGE) gel and stained with Coomassie blue staining, filter-aided sample preparation in-gel digestion was performed. The peptides were labeled with a tandem mass tag (TMT) (Thermo Fisher Scientific, USA), according to the manufacturer’s instructions. Before the liquid chromatography with tandem mass spectrometry (LC–MS/MS) analysis, the TMT-peptides were fractionated by the hi-pH (pH 10) Agilent 1260 Infinity II reversed-phase (RC)-HPLC system (Agilent, China). The peptide fractions were separated with the Easy-nLC ultra-high performance liquid chromatography (UHPLC) system (Thermo Fisher Scientific) equipped with an Acclaim PepMap rapid separation liquid chromatography (RSLC) column (Thermo Fisher Scientific) with a flow rate of 300 nL·min^–1 . AQ Exactive Plus MS (Thermo Fisher Scientific) was used for the MS analysis. Next, the raw data were searched with Mascot 2.6 (Matrix Science, UK) and Proteome Discoverer 2.1 (Thermo Fisher Scientific), and then matched against the proteome database (Uniprot_HomoSapiens_20386_20180905). KEGG analysis of the differential proteins was performed with the online database KEGG Orthology Based Annotation System (KOBAS) 3.0^↑↑.

^↑↑ kobas.cbi.pku.edu.cn.

《2.4. Transcriptomics》

2.4. Transcriptomics

Transcriptomics was performed for 5 ×10⁶ (per sample) miR516a-3p-transfected SNU-449 (n = 3) and their controls. In brief, the total RNAs were extracted with TRIzol reagent (Invitrogen, USA), and the messenger ribonucleic acids (mRNAs) were purified with poly-T oligo-attached magnetic beads. The mRNAs were reverse-transcribed as complementary deoxyribonucleic acids (cDNA) libraries using an mRNA-seq sample preparation kit (Illumina, USA). The cDNAs were purified with AMPureXP beads and amplified using polymerase chain reaction (PCR) amplification. The total sequencing libraries were then sequenced using the Illumina Hiseq X Ten system (Illumina). The Hisat package was used to assign the reads of samples by mapping the reads to the reference genome. StringTie was performed to assemble the mapped reads. KEGG analysis of the differential mRNAs was performed with KOBAS 3.0.

《2.5. Analysis of the target genes with proteomics and transcriptomics》

2.5. Analysis of the target genes with proteomics and transcriptomics

The potential target genes of miR-516a-3p were screened with the online analysis database TargetScan^↑ . All the target candidates were overlapped with the data of downregulated proteins or transcripts from proteomics or transcriptomics. Then the overlapped proteins or transcripts were set as the potential targets of the miRNA and highlighted as nodes. The systematic analysis of the target genes’ effects was conducted as follows: ① the correlation between the target genes and the differential genes was evaluated by means of the Pearson correlation test with the corrplot R package; ② the differential genes were divided into genesets according to the functional enrichment, and the networks between the target genes and the gene sets were constructed by means of the Mantel test based on the Bray–Curtis distance.

^↑ http://www.targetscan.org/cgi-bin/targetscan/vert_71/targetscan.cgi?mirg=hsamiR-516a-3p

《2.6. Analysis of the integrated multi-omics》

2.6. Analysis of the integrated multi-omics

Unlike the traditional statistical selection of differential variables, the data integration analysis for biomarker discovery using latent components (DIABLO) model is committed to selecting variables across multiple dimensions and multiple functional levels with multivariate exploratory approaches [25]. Simultaneously selecting features via L1 penalization makes the constructed signaling network models more biologically interpretable. Thus, DIABLO has been widely used to analyze cell biology, explore mechanisms, mine pathways, and construct structures. In this study, matrices from the transcripts, proteins, and metabolites were input to DIABLO in order to identify the major combinations of biological-related features. The model was composed of two components, and each component contained 90 features that were selected jointly across all three data types (2 ×30 features per data type). A module-based approach was used to transform the omics variables datasets into pathway datasets, using the weighted-gene correlation network analysis (WGCNA) R package. The mixOmics R package was used to integrate and visualize correlations among variables. The community structure was constructed in order to visually describe the degree of density and correlation among the variables, using the igraph R package. To construct the branchconnecting networks, Metscape 3, the internal relational database from KEGG and Edinburgh human metabolic network (EHMN), was used to integrate the gene expression (i.e., transcriptomics and proteomics) and metabolomics [26].

《2.7. Cell culture》

2.7. Cell culture

The HCC cell lines—namely, SNU-449, Hep G2, HUH-7, and LM3—were purchased from the American Type Culture Collection or the Cell Bank of the Chinese Academy of Sciences. Media were used for cell culture as follows: RPMI 1640 medium (GIBCO, USA) was used for SNU-449; Eagle’s minimum essential medium (Thermo Fisher Scientific) was used for Hep G2; and Dulbecco’s modified Eagle’s medium (Thermo Fisher Scientific) was used for HUH-7 and LM3. All media were supplemented with 10% fetal bovine serum (FBS) (BioIND, China). The cells were cultured in a humidified incubator under 5% carbon dioxide (CO₂) at 37 °C.

《2.8. Cell transfection》

2.8. Cell transfection

The miR-516a-3p mimic oligonucleotides, inhibitor oligonucleotides, and their controls (RiboBio) were synthesized for the in vitro experiment. RiboFECT (RiboBio) was used to efficiently transfect these RNA oligonucleotides (100 nM for the mimics and 200 nM for the inhibitor) according to the manufacturer’s instructions.

Lentiviral vector particles (Genepharma, China) were constructed to stably express the miR-516a-3p inhibitor for the in vivo experiment. In brief, the miR-516a-3p inhibitor sequence was cloned into third-generation packaging system vectors. The vectors were transfected into the 293T cells; then, lentiviral particles were collected and concentrated from the supernatants. 2 ×10⁵ LM3 cells were transduced with 1 ×10⁷ infectious lentiviral particle units (multiplicity of infection = 50). The LM3 cells, which stably expressed miR-516a-3p inhibitor, were screened with 2 μg·L^–1 of puromycin (Thermo Fisher Scientific) for two weeks.

《2.9. Quantitative reverse-transcription PCR (qRT-PCR)》

2.9. Quantitative reverse-transcription PCR (qRT-PCR)

The miRNeasy Mini Kit (No. 217004, Qiagen, Germany) and Rneasy MinElute Cleanup Kit (No. 74204, Qiagen) were used to extract and purify the total RNAs and miRNAs from the HCC cells, according to the manufacturer’s instructions. For the isolation of exosome RNAs for the cell supernatant, anexoRNeasy Maxi Kit (No. 77164, Qiagen) was used. ABulge-Loop miRNA qRT-PCR Starter Kit (C10211, Ribobio) was used to detect the expression of pre-miRNAs and mature miRNAs. In brief, for stem-loop miRNA qRT-PCR, specific reverse-transcription primers were used for the synthesis of the pre-miRNA or miRNA cDNA, at 42 °C for 60 min and at 70 °C for 10 min, respectively. PCR amplification was then performed for 40 cycles at 95 °C for 15 s and at 60 °C for 30 s per cycle (for pre-miRNAs, a temperature of 55 °C was used). U6 was used as the normalization control, and the relative expression levels were calculated with the cycling threshold (CT) value comparison  method. The reverse-transcription primers and PCR amplification primers for miR-516a-3p (MQPS0001726-1), miR-516a-5p (MQPS0001727-1), and miR-516b-5p (MQPS0001730-1) were synthesized and purchased from Ribobio. The other primer sequences are shown in Appendix A Table S1.

《2.10. Western blot analysis》

2.10. Western blot analysis

Western blot analysis was performed as described previously [27]. In brief, after being quantified and boiled, all the collected proteins were electrophoresed with 10% SDS-PAGE gel (FD341- 100, Fudebio, China), and then transferred onto equilibrated polyvinylidene-difluoride membranes. The membranes were then blocked for 1 h at room temperature and incubated with primary antibodies overnight at 4 °C. Before detection by an enhanced chemiluminescence (ECL) system (Biotanon, China) with FDbio-FemtoECL (FD8380, Fudebio), the membranes were incubated with secondary anti-rabbit immunoglobulin G (IgG) or anti-mouse IgG horseradish peroxidase (HRP)-linked antibody. All 15 antibodies are listed in Appendix A Table S2.

《2.11. Cell counting kit-8 (CCK-8) and 5-ethynyl-20 -deoxyuridine (EdU) assay》

2.11. Cell counting kit-8 (CCK-8) and 5-ethynyl-20 -deoxyuridine (EdU) assay

A CCK-8 assay and an EdU assay were used to detect cell proliferation. The HCC cell lines that had been transfected with the miR-516a-3p oligonucleotides were seeded in 96-well plates (for CCK-8) or 24-well plates (for EdU) for 24 h. The CCK-8 reagent (Dojindo Molecular Technologies, Japan) was added and incubated for 1 h (diluted with media, 1:10), and the absorbance at 450 nm was detected using a spectrophotometric reader (Multiskan FC, Thermo Fisher Scientific). A Click-iTEdU Alexa Fluor 488 Assay Kit (Invitrogen) was used to measure cell proliferation, according to the manufacturer’s instructions. Fluorescence images were taken with a fluorescence microscope (Photometrics Prime, USA). The proportion of green fluorescent cells to blue fluorescent cells reflected the efficiency of proliferation.

《2.12. Real-time cell analysis (RTCA) assay》

2.12. Real-time cell analysis (RTCA) assay

RTCA was performed to monitor the cell index (CI), which reflects cell proliferation and invasion in a label-free, real-time manner. AnxCELLigence RTCA single plate (SP) instrument equipped with E-Plate 96 was used for cell proliferation. An xCELLigence RTCA dual purpose (DP) (ACEA Biosciences, USA) instrument equipped with CIM-Plate 16 (pre-coated with Matrigel) was used to determine cell invasion. The HCC cells that had been transfected with miRNA oligonucleotides for 24 h were seeded in the E-Plate 96 or CIM-Plate 16; the CI was then detected in realtime, according to the manufacturer’s instructions.

《2.13. Colony formation assay》

2.13. Colony formation assay

After being transfected with miR-516a-3p oligonucleotides for 24 h, the HCC cells were seeded in the 12-well plates, and the media was changed every three days. After two weeks, the cells were fixed with paraformaldehyde and stained with crystal violet (Sigma-Aldrich, Germany).

《2.14. Transwell assay》

2.14. Transwell assay

The Transwell system (Corning Inc.) was used for the cell migration and invasion assay. In brief, the HCC cells were transfected with miR-516a-3p oligonucleotides for 24 h. For the invasion assay, Matrigel (BD Biosciences, USA) was pre-coated into the upper chamber of adorning Transwell membrane and incubated with serum-free medium for 2 h. The Matrigel was not used for the migration assay. The transfected cells were resuspended with serum-free medium in the upper chamber, and the lower chamber of the plate was contained with serum culture media. After cultivation for 24 h, the chambers were fixed with paraformaldehyde and stained with crystal violet. The stained cells were counted to assess their invasion ability.

《2.15. Dual-luciferase reporter assay》

2.15. Dual-luciferase reporter assay

As described previously [27], the sequences of the wild-type or mutant-type 3' untranslated regions (3' UTR) of the six target candidates (LMBR1, CHST9, RBM3, SLC7A6, PTGFRN, and NOL12) were respectively synthesized and cloned into a pmirGLO vector (Repobio, China). The 293 T cells were then co-transfected with the constructed vectors and the miR-516a-3p mimics or control mimics using lipofectamine 3000 (Thermo Fisher Scientific). After 48 h, the dual-luciferase activity of all the cells was detected with a dual-luciferase reporter assay kit (Vazyme, China), according to the manufacturer’s instructions.

《2.16. Cy3-labeled miR-516a-3p exosomal transfer detection》

2.16. Cy3-labeled miR-516a-3p exosomal transfer detection

Two forms of culture systems were used to detect the route of miR-516a-3p exosomal transfer. One form was a co-culture assay. After being transfected with Cy3-miR-516a-3p mimics, the cells were harvested and washed three times with serum-free media. Then they were seeded in the upper chamber of a Transwell system (0.4 μm, polycarbonate membrane, Corning Inc.), with the untreated HCC cells being pre-seeded in the lower chamber. After 24 h, the cells in the lower chamber were washed three times with phosphate-buffered saline (PBS). The Cy3 red fluorescence on the cells was observed. The other form of culture system was a media transfer. After being transfected with Cy3-miR-516a-3p mimics for 24 h in a 6-cm culture dish, HCC cells were washed three times and then cultured with serum-free media for 24 h. The media were collected and centrifugated at 3000 g at 4 °C for 5 min to remove cell debris. The supernatant was then transferred and cultured into a 24-well plate, with the untreated HCC cells being pre-seeded. After 24 h, the cells were washed three times with PBS and observed with fluorescence microscopy. One of each of the two forms of culture systems was treated with an inhibitor of exosome production GW4869 (20 μM, MedChemExpress, USA) and was used as the control group.

《2.17. Exosomal purification, identification, and co-culture》

2.17. Exosomal purification, identification, and co-culture

After the HCC cells (SNU-449 stably overexpressing miR-516a3p and untreated SNU-449) were grown in a 15 cm culture dish to 80% confluence, the primary media were replaced with 10 mL of serum-free media. After 24 h, the supernatant was centrifugated at 3000 g at 4 °C for 5 min and then filtrated with a Millex GV filter unit (0.22 μm, Merck Millipore, Germany). The exosomes in the supernatant were then extracted and purified using an exoEasy Maxi Kit (No. 76064, Qiagen), according to the manufacturer’s instructions. Finally, the exosomes that were absorbed in the membrane were eluted with 400 μL of elution buffer. An exoRNeasy Maxi Kit (No. 77164, Qiagen) was used for the isolation of the exosome RNA.

To observe the exosome structure, the exosomes were attached to a copper grid and negatively stained with phosphotungstic acid (2%) for 1.5 min. After the samples were dried for 15 min, the exosomes on the copper grid were observed with a transmission electron microscope (Tecnai G2 Spirit 120 kV, Thermo Fisher Scientific). The hydrodynamic diameter distribution of the exosomes was measured by a Nano-ZS90 instrument (Malvern, UK), according to the manufacturer’s instructions. A total of five samples were tested, and each sample was tested three times.

For the exosome and cell co-culture, the exosomes were eluted with an elution buffer. Then, 125 μL of exosomes-elution was added to a 6 cm culture dish containing HCC cells with 2.5 mL of media. The cells were harvested for further experiments after 24 h.

《2.18. Patients and specimens》

2.18. Patients and specimens

HCC cells and matching adjacent non-tumor tissues were obtained from patients who underwent hepatectomy at the First Affiliated Hospital, College of Medicine, Zhejiang University, China. A total of 82 HCC patients were enrolled: 28 (34.1%) female and 54 (65.9%) male, with a mean age of (54.6 ± 12.9) years. Of these, 44 (53.7%) had a low of the American Joint Committee on Cancer (AJCC) stage and 38 (46.3%) had a high of AJCC stage; 38 (46.3%) had positive hepatitis B virus (HBV) infection, and 44 (53.7%) had negative HBV infection. All patients were confirmed as having HCC by means of postoperative histopathology. This research protocol was approved by the Ethical Committee of the First Affiliated Hospital, School of Medicine, Zhejiang University. The participants provided written informed consent to participate in this study.

《2.19. Tumor model》

2.19. Tumor model

Balb/c male nude mice (6–8 weeks old) were obtained from the Shanghai Experimental Animal Center, Chinese Academy of Science, and received care according to the criteria of the National Institute Guide for the Care and Use of Laboratory Animals. Before the mice were injected with cells, they were randomly assigned to each group. To minimize the suffering of the mice, all were anesthetized and euthanized with CO₂ inhalation before being sacrificed. The animal research protocol was approved by the Ethical Committee of the First Affiliated Hospital, School of Medicine, Zhejiang University. Written informed consent was obtained from the owners for the participation of their animals.

The tumor model was performed as described previously [27]. In brief, LM3 cell lines (2 ×10⁶ ) that had been screened for the stable expression of miR-516a-3p inhibitor and their controls were subcutaneously injected into the right flank of the mice (six weeks, n = 6 in each group). After 28 days, the mice were sacrificed, and all tumor specimens were collected for cryopreservation or fixed preservation. To further evaluate the invasion and metastasis capacity in vivo, orthotopic HCC mouse models (eight weeks, n = 5 in each group) were used. LM3 cell lines with stable expression of miR-516a-3p inhibitor or their controls (1×10⁶ , 1:1 in Matrigel, 50 μL) were respectively grafted in the right lobe of the mice liver. After 42 days, invasion and metastasis were evaluated with living bioluminescence imaging.

《2.20. Statistical analysis》

2.20. Statistical analysis

The mean ± standard deviation (SD) or median ± interquartile range (IQR) was used for the description of the quantitative variables. Student’s t-test or the Mann–Whitney test was performed to analyze the comparison of the quantitative variables between the two groups. One-way analysis of variance (ANOVA) followed by a post hoc Bonferroni test was used among more than two groups. The optimal cut-off value ( = 3.979704795) for the expression of miR-516a-3p in the 82 HCC samples was detected with the ROC curve, by calculating the best Youden index, considering both sensitivity and specificity. Kaplan–Meier survival curves were assessed with the log-rank test. The risk factors of HCC prognosis were screened with a univariate and multivariate Cox proportional hazards regression analysis. All in vitro experiments were repeated independently, at least in triplicate. Statistical analysis was performed with statistical product and service solutions (SPSS) software (version 19.0), and a P value < 0.05 was set as the significance level. *, **, ***, and **** were used to represent P values of < 0.05, < 0.01, < 0.001, and < 0.0001, respectively. The statistical analysis for multi-omics integration was performed as described above.

《3. Results》

3. Results

《3.1. MiR-516a-3 poriginating from four oncogenic pre-miRNAs of C19MC positively correlated with advanced-stage HCC and independently predicted tumor recurrence and mortality risks》

3.1. MiR-516a-3 poriginating from four oncogenic pre-miRNAs of C19MC positively correlated with advanced-stage HCC and independently predicted tumor recurrence and mortality risks

Using the TCGA database, we previously demonstrated that several C19MC members were significantly upregulated in HCC [14]. In the current research, we further searched the miRBase database^↑ and found that four members of C19MC (mir-516a-1, mir-516a-2, mir-516b-1, and mir-516b-2) co-spliced an identical mature miRNA—namely, miR-516a-3p (Fig. S1 in Appendix A). First, based on data from TCGA, we confirmed that the four pre-miRNAs were significantly upregulated in HCC tissues in comparison with nontumor tissues (Fig. S2(a) in Appendix A). High expression of the four pre-miRNAs in HCC was found to be significantly correlated with poor patient survival (Figs. S2(b) and (c) in Appendix A). The high expression of the four pre-miRNAs in HCC tissues was further verified in our Chinese cohort (Fig. S3(a) in Appendix A). We further demonstrated that miR-516a-3p was significantly more highly expressed than the complimentary chains (miR-516a-5p and miR516b-5p) in four frequently used HCC cell lines (Figs. S3(b) and (c) in Appendix A). The results of our study indicated that miR-516a-3p was the functional mature miRNA from the four pre-miRNAs (Fig. S3(d) in Appendix A).

^↑ http://www.mirbase.org/

In our cohort, the expression of miR-516a-3p was significantly upregulated in HCC tissues (Fig. 1(a)). We further correlated tumor miR-516a-3p expression with HCC features and found significantly higher tumor miR-516a-3p levels in tumor sizes greater than 8 cm (vs < 8 cm), in multiple tumors (vs single tumor), and in AJCC stages III–IV (vs I–II) (Figs. 1(b)–(d)). Moreover, patients with high tumor miR-516a-3p levels showed a significantly lower tumor-free survival rate and lower overall rate of survival, compared with those with low miR-516a-3p levels (Figs. 1(e) and (f)). In a multivariate Cox analysis, tumor miR-516a-3p expression was shown to be an independent influencing factor of HCC recurrence and HCC mortality (Figs. 1(g) and (h)). For predicting the probability of recurrence, the area under the curve (AUC) was 0.673 and 0.641 for AJCC staging and tumor miR-516a-3p expression, respectively. When the two variables were combined, the AUC showed a significant increase to 0.722 (all < 0.05) (Fig. 1(i)). We further confirmed that high tumor miR-516a-3p expression also improved the prediction of HCC mortality. When AJCC staging and miR-516a-3p expression were combined, the AUC significantly increased from 0.740 to 0.792 (all P < 0.05) (Fig. 1(j)).

《Fig. 1》

Fig. 1. MiR-516a-3p is positively correlated with HCC malignancy. (a) MiR-516a-3p is upregulated in HCC tissues in comparison with the adjacent non-tumor tissues (n = 25) (Student’s paired t-test); (b)–(d) high tumor miR-516a-3p expression is positively correlated with (b) large tumor size, (c) large tumor number, and (d) high AJCC staging (n = 82) (Mann–Whitney test); (e, f) the Kaplan–Meier and the log-rank test confirm that high tumor miR-516a-3p expression distinguishes (e) poor HCC tumor-free survival and (f) OS; (g, h) the multivariate Cox proportional hazards model demonstrates that tumor miR-516a-3p expression is an independent risk factor in HCC (g) recurrence and (h) mortality; (i, j) the AUC of the ROC shows that tumor miR-516a-3p expression improves the prediction of (i) HCC recurrence and (j) mortality. HR: hazard ratio; AUC: area under the curve.

《3.2. MiR-516a-3p increased the proliferation, invasiveness, and metastasis of HCC cells in vitro and in vivo》

3.2. MiR-516a-3p increased the proliferation, invasiveness, and metastasis of HCC cells in vitro and in vivo

We performed a set of tests to evaluate the effect of miR-516a3p on the proliferation of HCC in vitro (Fig. 2 and Figs. S4 and S5 in Appendix A). All the tests (CCK-8, RTCA, EdU, and colony formation) showed that overexpression of miR-516a-3p significantly promoted the proliferation of HUH-7 and SNU-449 cells (Fig. 2(a) and Figs. S4(a), (c), and (e) in Appendix A), while the knockdown of miR-516a-3p inhibited the proliferation of HCC cells (Fig. 2(b) and Figs. S4(b), (d), and (f) in Appendix A). We performed a Transwellassay and observed that overexpression of miR-516a-3p significantly increased the migration and invasiveness of HUH-7 and SNU-449 cells (Figs. S5(a) and (c) in Appendix A), while knockdown of miR516a-3p had the opposite effects (Figs. S5(b) and (d) in Appendix A). RTCA with the Transwell system further verified the effect of miR-516a-3p on the invasion of HCC cells (Figs. 2(c) and (d)).

Furthermore, a subcutaneous nude mice HCC model was established using LM3 cell lines, which stably expressed miR-516a-3p inhibitor or its controls. The mice were sacrificed four weeks after subcutaneous injection and the tumor volume was measured. The results revealed that the miR-516a-3p inhibitor significantly inhibited the growth of HCC (Fig. 2(e)). Compared with the control group, the miR-516a-3p inhibitor resulted insmaller tumor volume (Fig. 2(f)) and tumor weight (Fig. 2(g)). To assess the role of miR516a-3p in tumor metastasis, an HCC in situ mice model was established. ALM3 cell line stably expressing miR-516a-3p inhibitor or its controls was implanted in the right liver of nude mice for six weeks. Luciferase imaging was applied to evaluate the growth and metastasis of tumors; it showed that the miR-516a-3p inhibitor group had a lower tumor burden in both the liver (primary tumor) and lung (metastatic tumor) compared with the control group (Figs. 2(h)–(j) and Fig. S6 in Appendix A).

《Fig. 2》

Fig. 2. MiR-516a-3p promotes the malignancy of HCC cells in vitro and in vivo. (a, b) RTCA records the cell proliferation index of HUH-7 and SNU-449 after transfection with (a) miR-516a-3p mimics, (b) miR-516a-3p inhibitor, or their controls (Ctrls); (c, d) change in CI reflects the invasion of HUH-7 and SNU-449 after being transfected with (c) miR-516a-3p mimics or (d) inhibitor, by RTCA; (e)–(g) a xenograft mice model was established by using LM3 cells stably expressing miR-516a-3p inhibitor and Ctrls; (h)–(j) bioluminescence imaging of the HCC in situ mouse models reflects the abdominal tumor burden and lung metastasis after being seeded with LM3 cells stably expressing miR-516a-3p inhibitor or the Ctrls.

《3.3. MiR-516a-3p was delivered among HCC cells via exosomes and increases the oncogenic activity of recipient cells》

3.3. MiR-516a-3p was delivered among HCC cells via exosomes and increases the oncogenic activity of recipient cells

We co-cultured SNU-449 cells transfected with Cy3-labeled miR-516a-3p (donor cells) with untreated SNU-449 cells (recipient cells). We observed that the recipient cells contained intense Cy3- miR-516a-3p fluorescence after 24 h of co-culture (Fig. S7 in Appendix A), which suggested that miR-516a-3p could be transferred among HCC cells.

We further assessed the exosome-mediated miR-516a-3p delivery among HCC cells. We added GW4869 into the co-culture system and observed that the Cy3-fluorescence of the recipient cells was dramatically weakened (Fig. 3(a)). Moreover, we cultured the donor cells in a serum-free medium with GW4869 or dimethyl sulfoxide (DMSO) for 24 h, and then collected the medium and transferred it to the recipient cells. We found that the uptake of Cy3-fluorescence of the recipient cells was also weakened by GW4869 (Fig. 3(b)).

《Fig. 3》

Fig. 3. MiR-516a-3p is transferred among HCC cells via exosomes. (a) After the co-culture systems were treated with GW4869 or DMSO for 24 h, the recipient cells (green dashed line) were observed using red fluorescent signals (Cy3); (b) after the donor cells were cultured with GW4869 or DMSO for 24 h, the medium was transferred to the recipient cells (green dashed line), red fluorescent signals (Cy3) were observed; (c) real-time PCR analysis of miR-516a-3p expression level from the ^exo516a-3p exosomes, compared with the ^exoCtrl exosome; (d) real-time PCR analysis of miR-516a-3p expression level from the exosomes of miR-516a-3p-overexpressing SNU-449 cells treated with GW4869, DMSO, or mock; (e) real-time PCR analysis of miR-516a-3p expression level from the exosomes of miR-516a-3p-overexpressing SNU-449 cells treated with RNase or RNase combined with Triton X-100; (f) SNU-449 cells and HUH-7 (green dashed line) were stained with FITC-lysosome and were observed with a confocal fluorescence microscope, arrows indicate the co-localization of Cy3-labeled miR-516a-3p and lysosome; (g) schematic diagram demonstrated that exosomal transfer of miR516a-3p was achieved among HCC cells. DAPI: 4,6-diamino-2-phenyl indole; DMSO: dimethyl sulfoxide; FITC: fluorescein isothiocyanate; RNase: ribonuclease; ^exo516a-3p: isolated exosomes; exoCtrl: ^exosomal control.

Next, we isolated exosomes (^exo516a-3p) from miR-516a-3poverexpressing SNU-449 cells. Typical exosomes were identified as having a diameter of about 100 nm by transmission electron microscopy (Fig. S8(a) in Appendix A). Dynamic light scattering (DLS) was used to describe the distribution of the hydrated particles (Fig. S8(b) in Appendix A). Positive exosome markers (CD9, CD81, CD63, TSG101, HSP70, and HSP90) further confirmed the identity of the purified exosomes (Fig. S8(c) in Appendix A). We extracted exosomal miRNAs and found that miR-516a-3p was largely enriched in ^exo516a-3p, compared with in the exosomal control (^exoCtrl) group (Fig. 3(c)). The release of ^exo516a-3p was significantly inhibited by GW4869 (Fig. 3(d)). We further used Triton X-100 to destroy the membrane of the exosomes and release the contents encased in the exosomes. The exosomal miR-516a-3p was then degraded using ribonuclease (RNase) (Fig. 3(e)). Interestingly, after the recipient cells had been treated with the medium from the donor cells, we labeled the lysosomes with fluorescein isothiocyanate (FITC)-lysosomal fluorescence and found that Cy3-miR-516a-3p was co-localized with FITC-lysosome (Fig. 3(f)).

We further evaluated the role of exosome-mediated miR-516a3p transfer in HCC cells. We incubated SNU-449 and HUH-7 cells with purified ^exo516a-3p for 24 h and found a significantly increased miR-516a-3p expression in these cells (Fig. S9(a) in Appendix A). The recipient cells displayed enhanced tumor malignancy, including colony formation, cell proliferation, migration, and invasiveness (Figs. S9(b)–(e) in Appendix A). Therefore, we propose that exosome-mediated delivery plays a vital role in miR-516a-3p-induced oncogenic activity among HCC cells (Fig. 3(g)).

《3.4. A multiomics network discovery strategy revealed the relative correlation and differentiation of the three ‘‘omics”》

3.4. A multiomics network discovery strategy revealed the relative correlation and differentiation of the three ‘‘omics”

We performed transcriptomics, proteomics, and metabolomics analyses to map a multi-omics miR-516a-3-regulated molecular network. A model of the DIABLO method for omics studies—a multi-omics method for integrating correlated biological processes and discriminating among differential phenotypes—was performed to successfully select the variables representative of each omics (Figs. S10(a) and (b) in appendix A). Next, we carried out a partial least squares-discriminant analysis (PLS-DA) and created cluster heatmaps to demonstrate that the variable expression profiling from the datasets could distinguish the miR-516a-3poverexpressing SNU-449 cells from the control cells (Fig. 4(a) and Fig. S10(c) in Appendix A). Then we constructed an overview of the dataset correlation and finalized the strong correlations, which were along the ‘‘central dogma”(Fig. 4(b)). From the variable correlation analysis, the Circosplot demonstrated that the variables from the three omics were cross-correlated (Fig. 4(c)). From the correlation circle plots, which provided better visibility regarding the degree of correlation, we found a cluster-based consistency between the proteomics and metabolomics datasets. In contrast, the transcriptomics variables were relatively separated from the other two datasets (Fig. 4(d)). These results were further validated by means of a component analysis, when we transformed the features of the multi-omics variables into visualized communities. The variables in proteomics and metabolomics were clustered tightly, while the variables in transcriptomics were relatively differentiated against with the other two datasets (Fig. 4(e)).

《Fig. 4》

Fig. 4. The transcriptomics were relatively correlated and differentiated with the proteomics and metabolomics. (a) Hierarchical heatmap of the selected variables from the three omics; (b) an overview correlation construction showed a strong correlation between the proteomics, transcriptomics, and metabolomics, with the selected variables; (c) circosplot representing the mutual correlation of the selected variables in different omics; (d) circle plot highlighting the contribution of each selected variable to the relative correlation, the location of the points of each variable indicates the degree of correlation and discrimination; (e) community-connected networks of the integrated multi-omics in the variable levels, close-correlated variables were automatically clustered into various communities according to the degree of correlation.

《3.5. Deep analysis of proteomics revealed that miR-516a-3p regulates the metabolism of HCC cells by directly inhibiting six core targets》

3.5. Deep analysis of proteomics revealed that miR-516a-3p regulates the metabolism of HCC cells by directly inhibiting six core targets

We performed proteomics screening to explore the miR-516a3p-influenced molecules. There were 248 proteins with significantly differential expression between miR-516a-3poverexpressing SNU-449 cells and the control cells (Fig. 5(a)). KEGG analysis showed that the differentially expressed proteins were significantly enriched in metabolic pathways, including carbon metabolism, the biosynthesis of amino acids, pyrimidine metabolism, purine metabolism, glycolysis/gluconeogenesis, and the pentose phosphate pathway (Fig. 5(b)). It was notable that overexpression of miR-516a-3p comprehensively upregulated the proteins that were enriched in metabolic pathways (Fig. 5(c)).

To identify the direct targets of miR-516a-3p, we overlapped the downregulated proteins with the potential targets of miR-516a-3p (TargetScan) and found 16 candidates (Fig. 5(d)). We further performed Mantel’s test to correlate the potential targets with the enriched pathways. Finally, six proteins (LMBR1, CHST9, RBM3, SLC7A6, PTGFRN, and NOL12) were shown to be the direct targets of miR-516a-3p, making up a protein–protein interaction network that plays a central role in regulating the enriched metabolic pathways (Fig. 5(d)). Moreover, the six proteins were significantly correlated with the other differentially expressed proteins, indicating their central role in miR-516a-3p-interacting molecules (Fig. S11 in Appendix A). We further used experimental methods to confirm that all six proteins were the direct targets of miR-516a-3p in HCC. First, under the control of miR-516a-3p, the expression levels of the six targets were downregulated, as detected by Western blot analysis (Fig. S12 in Appendix A). Furthermore, luciferase reporter vectors containing wild-type or mutant-type 3' UTR of the six target candidates (LMBR1, CHST9, RBM3, SLC7A6, PTGFRN, and NOL12) were constructed (Fig. S13(a) in Appendix A), and a dual-luciferase reporter assay was performed. After being co-transfected with miR-516a-3p and wild-type 3' UTR vectors, the activities of luciferase were inhibited to various degrees, which did not occur with the mutant-type counterparts (Fig. S13(b) in Appendix A). Thus, proteomics combined with experimental evidence provided new insight into the investigation of miRNA targets.

《Fig. 5》

Fig. 5. Proteomics reveals that miR-516a-3p regulates the metabolism of HCC cells by directly inhibiting six core targets. (a) Differential proteomics variables depicting the differences between the miR-516a-3p group and the miR-Ctrl group in a volcano plot; (b) bubble chart showing the top 20 significant KEGG items from the proteomics; (c) expression levels of the proteins enriched in the key metabolism pathways; (d) colorful lines show the correlations between the target proteins and the enriched pathways from proteomics, the blue triangle diagram shows the correlations among the target proteins, the "P1," "P2," and "P3" highlight the strong, moderate, and weak correlations of the targets, respectively. HIF-1: hypoxia inducible factor-1.

We also compared the differential transcriptome profiles (Figs. S14(a) and (b) in Appendix A). KEGG analysis revealed the significant enrichment of differential genes in metabolic pathways (Fig. S14(c) in Appendix A). We overlapped the downregulated genes with the potential targets of miR-516a-3p (TargetScan 7.2) and found 22 candidates (Fig. S14(d) in Appendix A). However, none of them significantly correlated with the enriched pathways and other differentially expressed genes (Figs. S14(d) and (e) in Appendix A).

《3.6. MiR-516a-3p aberrantly accumulated purine and pyrimidine metabolites in HCC cells》

3.6. MiR-516a-3p aberrantly accumulated purine and pyrimidine metabolites in HCC cells

We carried out non-targeted metabolomics to verify the effects of miR-516a-3p on the metabolism regulation of HCC. Principal component analysis (PCA) and PLS-DA were used for quality control (Fig. S15(a) in Appendix A). There were 81 and 92 significant differentially expressed compounds in positive and negative ion modes, respectively, between miR-516a-3p-overexpressing SNU-449 cells and the control cells (Fig. 6(a) and Fig. S15(b) in Appendix A). Mummichog analysis using MetaboAnalyst 4.0 showed that the metabolic features induced by miR-516a-3p overexpression were associated with purine metabolism, pyrimidine metabolism, sphingolipid metabolism, the citric acid cycle, tryptophan metabolism, glycerophospholipid metabolism, linoleic acid metabolism, arachidonic acid metabolism, and the glyoxylate and dicarboxylate metabolic pathways (Fig. 6(b)).

We further overlapped the enriched pathways in proteomics and metabolomics. The co-pathway analysis revealed a potentially central role of purine and pyrimidine metabolism in the miR-516a3p-mediated network (Fig. 6(c)). Overexpression of miR-516a-3p significantly upregulated the proteins (EEF2, PAICS, UCK2, NME1, NME2, CMPK1, etc.) and metabolites in the purine and pyrimidine pathway (uridine, guanosine, inosine, guanine, deoxyguanosine, and hypoxanthine) (Fig. 6(d)).

《Fig. 6》

Fig. 6. The combination of metabolomics and proteomics demonstrates that miR-516a-3p aberrantly accumulates purine and pyrimidine metabolites in HCC cells. (a) Differential metabolites from the positive (left) and negative (right) ion modes of metabolomics both distinguish the miR-516a-3p group and miR-control group, in volcano plots; (b) bubble chart showing significant pathway items from metabolomics using Mummichog analysis with MetaboAnalyst 4.0; (c) the branch network visually exhibits the activation of the purine and pyrimidine pathways in the proteins and metabolites levels, controlled by miR-516a-3p; (d) heatmaps showing that miR-516a-3p upregulates the expression levels of the key proteins and metabolites enriched in the purine and pyrimidine metabolism pathway. FDR: false discovery rate; FC: fold change.

《4. Discussion》

4. Discussion

Many approaches have been used to explore tumor-specific miRNAs. The TCGA database has been widely investigated due to its large sample size, diverse tumor types, and sufficient clinical data. Through the TCGA miRNA profile, we screened four members of C19MC that were in the Top ranks among significantly deregulated pre-miRNAs in HCC and that predicted poor patient survival. Two previous studies have demonstrated that other members of C19MC (miR-519d, miR-519a, and miR-519c-3p in one study, and miR-517a and miR-510c in the other study) acted as oncogenes in the development and progression of HCC [16,17]. Increasing pieces of evidence have revealed that miRNA cluster members can act on the same tumor regulatory pathway and display similar effects [9,28]. Therefore, we suggest that the miRNA cluster be viewed as acting together (i.e., as a ‘‘team”), which can be likened to several ants working collectively to perform a single function.

Mature miRNAs, rather than their pre-miRNAs, play roles in cellular physiopathology. Interestingly, the four oncogenic premiRNAs co-spliced the same mature miRNA—namely, miR-516a3p. We then focused on miR-516a-3p and observed its significant clinical value. High levels of miR-516-3p in a tumor were found to be linked with an advanced stage of HCC and predicted poor patient prognosis. The incorporation of tumor miR-516a-3p levels into the AJCC staging system sharply improved the diagnostic accuracy in predicting postoperative HCC recurrence and patient survival. The in vitro and in vivo experiments further demonstrated miR-516a-3p to be oncogenic, resulting in enhanced malignant behaviors, including proliferation and invasion. Therefore, miR-516a-3p is suggested as a novel target for the treatment of HCC that improves prognostic prediction.

We further demonstrated that miR-516a-3p affects not only the host or donor cancer cells, but also others acting in transfer. It can be delivered via exosomes and thus increases the oncogenicity of the recipient cells. Blocking cellular communication by means of an exosome inhibitor decreased miR-516a-3p transportation as well as malignant behavior delivery. Notably, we observed that the transferred miR-516a-3p was incorporated into the lysosome of the recipient cells, which indicated that the lysosome participates in the release of exosome contents. These results are consistent with previous reports [29,30]. These findings are helpful for understanding the cooperation among cancer cells, and can provide a therapeutic strategy to control cellular interactions during cancer progression.

We performed a multi-omics analysis for the functional exploration of miR-516-3p. The enrichment of significantly differentially expressed genes, proteins, and metabolites revealed that miR-516- 3p-triggered oncogenic signaling played a central role in metabolic regulation. Through a deep analysis of proteomics, we constructed a comprehensive protein–protein interaction network and identified six key players in the aberrant metabolism regulated by miR-516-3p: the proteins LMBR1, CHST9, RBM3, SLC7A6, PTGFRN, and NOL12. Among them, LMBR1 [31] and RBM3 [32] have been reported to be involved in insulin signaling transduction and glucose metabolism, while PTGFRN is associated with lipid drop accumulation [33], and NOL12 has been found to be essential in RNA metabolism and DNA maintenance [34]. SLC7A6 is anamino acid transporter [35]. In this study, we confirmed that miR-516-3p directly targeted these six core proteins. We further constructed their closely linked protein network and subsequently assessed how they modulated cancer metabolisms such as carbon metabolism and the purine and pyrimidine metabolism. Therefore, the target-pathway interaction network that was generated by means of proteomics significantly demonstrates the mechanisms and functions of miRNA from a new perspective.

Metabolomics permits a better understanding of cancer metabolic dysregulation. The results of this study were highly consistent with proteomics-derived data showing that miR-516-3p activates metabolic pathways, including purine metabolism, pyrimidine metabolism, amino acid metabolism, and unsaturated fatty acids metabolism. The integrated analysis of the proteomics and metabolomics datasets further revealed a pivotal role of purine and pyrimidine metabolism in miR-516a-3p-induced metabolic dysregulation. Purine and pyrimidine metabolism involves aberrant DNA replication and RNA synthesis, which contributes to uncontrolled cancer cell proliferation and aggressive phenotype maintenance [36,37]. A boost in purine and pyrimidine metabolism acts as a hallmark in carcinogenesis and cancer development [37,38]. Ma et al. [39] confirmed that the activation of purine metabolism promotes the progress of HCC. Santana-Codina et al. [40] showed that oncogenic kirsten rat sarcoma viral oncogene (KRAS) boosts pancreatic cancer by triggering nucleotide synthesis. Feng et al. [41] demonstrated that UHMK1 promotes gastric cancer by activating purine synthesis. In addition, the enriched metabolic pathways identified in this study were closely linked with each other and had a synergistic effect on malignant behaviors. For example, the activation of one-carbon metabolism and the pentose phosphate pathway can provide essential carbon sources for the de novo synthesis of purine and pyrimidine metabolites [42]. However, in this study, the metabolites enriched in the carbon metabolism pathways could not be identified in non-targeted metabolomics, which is the inevitable bias of this approach [43].

We believe that the multi-omics strategy is a very promising approach to unravel the key biological functions of miR-516-3p. Because miRNA inhibits targets by translational repression or mRNA degradation, most of the functional exploration of miRNA started from the identification of target mRNA according to the matched seed regions and ended at a specific pathway. However, it is well known that miRNA has numerous potential targets and conducts biological activities via multiple pathways. A multi-omics study not only provides a highly efficient way to detect direct and indirect targets, but also presents a comprehensive molecular screening network. An integrative analysis can further identify the core mechanisms from the complex network. In addition, although no valuable targets were presented in the single dimension of transcriptomics, it is very interesting to find the deviation of transcriptomics data from proteomics and metabolomics data using a cluster-based algorithm, indicating a divergence between transcriptional and translational control in cellular pathophysiology. The interactions between large proteins and small metabolites play a major role in controlling a variety of cellular processes [44], and therefore deserve more attention.

In our view, it should be noted that the use of transcriptomics to explore the targets of miRNAs has a disadvantage. In this study, we failed to determine whether the target candidates from transcriptomics significantly correlated with the enriched pathways and other differentially expressed genes. This may be because the effects of miR-516a-3p depend on the inhibition of mRNA translation, rather than the degradation of mRNA. Thus, proteomics is preferred as a technology for the study of miRNA targets.

This study had some limitations. First, the clinical value of miR516-3p needs to be verified in various ethnic groups with large sample sizes. Second, our multi-omics analysis only provided a bird’s eye view of the miR-516-3p-mediated comprehensive molecular network. The details of the specific pathways targeting certain cellular metabolisms and oncogenic activities require further confirmation. Third, non-targeted metabolomics cannot fully detect the metabolites in amino acid metabolism or carbon metabolism. For analyzing a specific metabolic pathway and avoiding bias, targeted metabolomics is preferred.

In conclusion, our results provide evidence that the analysis of co-spliced mature miRNA products may be a new strategy to investigate the function of miRNA clusters. This study reveals the oncogenicity of four key members of C19MC by exploring one mature miRNA, in what might be called ‘‘killing four birds with one stone.” MiR-516a-3p can enhance tumor malignant behaviors and affect neighboring HCC cells via the exosome delivery system. This finding is expected to provide a novel molecular target for HCC therapy and prognostic prediction. Furthermore, to the best of our knowledge, this is the first study to explore the function of miRNA using a multi-omics approach, which has many advantages. This integrative analysis sheds new light on the miRNA regulation of cancer metabolism and malignant activity.

《Acknowledgments》

Acknowledgments

This study was supported by the National Natural Science Foundation of China (81771713 and 81721091) and the Zhejiang Provincial Natural Science Foundation of China (LR18H030001)

《Authors’ contributions》

Authors’ contributions

Tao Rui, Qi Ling, and ShusenZheng performed and organized the study. Shusen Zheng and Qi Ling provided funding support. Tao Rui performedmost of the experiments. Xueyou Zhang and Shi Feng participated in bioinformatic analyses. Lin Zhou and Haiyang Xie contributed to the analysis of clinical information. Tao Rui, Qi Ling, ShaoweiZhan, and Haitao Huang finished the statistical analysis. Tao Rui and Qi Ling finished the manuscript. Shusen Zheng reviewed the manuscript. All authors have approved the final manuscript.

《Compliance with ethics guidelines》

Compliance with ethics guidelines

Tao Rui, Xueyou Zhang, Shi Feng, Haitao Huang, Shaowei Zhan, Haiyang Xie, Lin Zhou, Shusen Zheng, and Qi Ling declare that they have no conflict of interest or financial conflicts to disclose.

《Appendix A. Supplementary data》

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2021.07.020.