Identification of circRNA Biomarker for Gastric Cancer through Integrated Analysis

Gastric cancer (GC) is one of the most common malignant tumors and ranks third in cancer mortality globally. Although, a lot of advancements have been made in diagnosis and treatment of gastric cancer, there is still lack of ideal biomarker for the diagnosis and treatment of gastric cancer. Due to the poor prognosis, the survival rate is not improved much. Circular RNAs (circRNAs) are single-stranded RNAs with a covalently closed loop structure that don’t have the 5′-3′ polarity and a 3′ polyA tail. Because of their circular structure, circRNAs are more stable than linear RNAs. Previous studies have found that circRNAs are involved in several biological processes like cell cycle, proliferation, apoptosis, autophagy, migration and invasion in different cancers, and participate in some molecular mechanisms including sponging microRNAs (miRNAs), protein translation and binding to RNA-binding proteins. Several studies have reported that circRNAs play crucial role in the occurrence and development of different types of cancers. Although, some studies have reported several circRNAs in gastric cancer, more studies are needed in searching new biomarkers for gastric cancer diagnosis and treatment. Here, we investigated potential circRNA biomarkers for GC using next-generation sequencing (NGS) data collected from 5 paired GC samples. A total of 45,783 circRNAs were identified in all samples and among them 478 were differentially expressed (DE). The gene ontology (GO) analysis of the host genes of the DE circRNAs showed that some genes were enriched in several important biological processes, molecular functions and cellular components. The Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis revealed that some host genes were enriched in several GC related pathways. The circRNA-miRNA-gene interaction network analysis showed that two circRNAs circCEACAM5 and circCOL1A1 were interacted with gastric cancer related miRNAs, and their host genes were also the important therapeutic and prognostic biomarkers for GC. The experimental results also validated that these two circRNAs were DE in GC compared to adjacent normal tissues. Overall, our findings suggest that these two circRNAs circCEACAM5 and circCOL1A1 might be the potential biomarkers for the diagnosis and treatment of GC.

Gastric cancer (GC) is one of the most common malignant tumors and ranks third in cancer mortality globally. Although, a lot of advancements have been made in diagnosis and treatment of gastric cancer, there is still lack of ideal biomarker for the diagnosis and treatment of gastric cancer. Due to the poor prognosis, the survival rate is not improved much. Circular RNAs (circRNAs) are single-stranded RNAs with a covalently closed loop structure that don't have the 5′-3′ polarity and a 3′ polyA tail. Because of their circular structure, circRNAs are more stable than linear RNAs. Previous studies have found that circRNAs are involved in several biological processes like cell cycle, proliferation, apoptosis, autophagy, migration and invasion in different cancers, and participate in some molecular mechanisms including sponging microRNAs (miRNAs), protein translation and binding to RNA-binding proteins. Several studies have reported that circRNAs play crucial role in the occurrence and development of different types of cancers. Although, some studies have reported several circRNAs in gastric cancer, more studies are needed in searching new biomarkers for gastric cancer diagnosis and treatment. Here, we investigated potential circRNA biomarkers for GC using next-generation sequencing (NGS) data collected from 5 paired GC samples. A total of 45,783 circRNAs were identified in all samples and among them 478 were differentially expressed (DE). The gene ontology (GO) analysis of the host genes of the DE circRNAs showed that some genes were enriched in several important biological processes, molecular functions and cellular components. The Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis revealed that some host genes were enriched in several GC related pathways. The circRNA-miRNA-gene interaction network analysis showed that two circRNAs circCEACAM5 and circCOL1A1 were interacted with gastric cancer related miRNAs, and their host genes were also the important therapeutic and prognostic biomarkers for GC. The experimental results also validated that these two circRNAs were DE in GC compared to adjacent normal tissues. Overall, our findings suggest that these two INTRODUCTION Gastric cancer (GC) is the fifth most common malignant tumors and third leading cause of cancer related deaths worldwide (Bray et al., 2018). Despite of having advanced diagnosis and treatment strategies, the prognosis of gastric cancer is still poor, and a 5-year survival rate is still less than 30% (Allemani et al., 2015). Due to poor prognosis, most of the patients are diagnosed at advance stages and are not capable of receiving the surgical therapy (Karimi et al., 2014). Therefore, new molecular diagnostic biomarker and therapeutic targets are of great interest to understand the molecular mechanism of gastric cancer diagnosis and treatment for early detection of patients and improving the overall survival rate.
Circular RNAs (circRNAs) are a class of single-stranded RNAs formed by a covalently closed loop structure without 5′-3′ polarities or 3′ polyA tail (Jiao Li et al., 2020). Although most of the circRNAs are generated from back splicing of exons, some are generated from introns and some are from both of exons and introns (Wang and Dong, 2019). Because of the circular structure, circRNAs are more stable than linear RNAs, and are abundant in different species. CircRNAs have the vast potentiality to be a molecular biomarker for diagnosis, treatment and prognosis of different cancers. Previous studies have demonstrated that circRNAs are involved in several biological processes including cell cycle, proliferation, apoptosis, autophagy, migration and invasion, in different types of cancers Sollott, 2017;Meng et al., 2018;Yu et al., 2018;Zhang et al., 2018). In recent years, circRNAs have been found to participate in several molecular mechanism like sponging microRNAs (miRNAs), protein translation and binding to RNA-binding proteins (Hansen et al., 2013;Kristensen et al., 2018;Chen et al., 2019;Zang et al., 2020). Moreover, some circRNAs can bind protein/peptides (Legnini et al., 2017;Pamudurti et al., 2017). Several studies have identified that circRNAs play key role in the occurrence and progression of a lot of malignancies like glioma (Li et al., 2018), hepatocellular carcinoma , pancreatic carcinoma (Huang et al., 2018), gastric cancer (Chen et al., 2017), colon cancer , prostate cancer (Tucker et al., 2020), breast cancer  etc. Despite several circRNAs have been reported in GC, more studies are needed to identify new biomarker for the diagnosis and treatment of this cancer.
In this study, to identify potential circRNA biomarker, we used next-generation sequencing (NGS) data from gastric cancer tissue and adjacent normal tissues. Then, differentially expressed circRNAs were identified and circRNA-miRNA-gene interaction network was constructed. GO term and KEGG pathway analyses were performed for the host genes of the DE circRNAs. GO term analysis showed that several genes functioned in some biological processes, molecular functions and cellular components. The KEGG pathway analysis showed that some genes were involved in several gastric cancer related pathways. We found two circRNAs circCEACAM5 and circCOL1A1 which were differentially expressed in all samples and these circRNAs were highly interacted with gastric cancer related miRNAs. In addition, the host genes (CEACAM5 and COL1A1) of these circRNAs were the important therapeutic and prognostic biomarker for gastric cancer (Hu and Chen, 2012;Zhou et al., 2015;Zhaoxing Li et al., 2020). The important properties of circRNAs includes sponging miRNAs and regulating gene transcription. As the host genes of the two circRNAs were important biomarker of GC and interacted with GC related miRNAs, they might have the potential to be important biomarker for GC. The expression of these two circRNAs were also validated by qRT-PCR experiments. Our results suggest that the circRNAs circCEACAM5 and circCOL1A1 might be the potential biomarkers for the diagnosis and treatment of GC.

Sample Collection
Fifteen pairs (10 for RT-PCR and 5 for RNA sequencing) of nonneoplastic gastric tissues (NT) and GCs tissues from Shenzhen Second People's Hospital (The First affiliated Hospital of Shenzhen University) were examined in the study. All tissues received no radiotherapy or chemotherapy before surgery and stored in RNAlater immediately after surgery. All patients provided written informed consent, approved by the ethics committee of Shenzhen University School of medicine.

Ribonucleic Acid Library Preparation
Total RNA was extracted from 5 paired GC samples. Ribosomal RNA was digested using ribo zero kit and linear RNA was removed using rnase R. CircRNA was fragmented and cDNA was synthesized with six base random primers and purified. Strand-specific library is prepared via incorporating a chemical label deoxy-UTP (dUTP) during production of the second-strand cDNA. The second-strand cDNA was specifically digested by UNG enzyme. The constructed library was analyzed by Agilent 2,100 Bioanalyzer and was sequenced on an Illumina Hiseq 2,500 (Chi Biotech, Shenzhen).

Ribonucleic Acid-Seq Data Analysis
We used 10 samples from 5 gastric cancer patients to investigate potential circRNA biomarker. From each patients, two samples were collected: one from the cancer tissue and other from the adjacent normal tissue. After getting the raw data in fastq format, Frontiers in Molecular Biosciences | www.frontiersin.org March 2022 | Volume 9 | Article 857320 the quality control analysis was performed. The quality control was done using the NGSQCToolkit (Patel and Jain, 2012) and the quality score 20 was used as the cutoff point. Filtering/trimming was carried out if the quality of the reads failed to reach the cutoff point. The fastq reads were aligned to the human reference genome (version GRCh38) using the BWA (Li and Durbin, 2009) aligner. CircRNAs were identified using the software CIRI (version 2) (Gao et al., 2015). The identified circRNAs were then annotated with the gene annotation file corresponding to the reference genome and the full length circRNA sequences were extracted. The full length circRNA sequences for all circRNAs were considered as the reference genome and the fastq reads were mapped using the bowtie2 (Langmead and Salzberg, 2012) aligner. Then count data was generated using bedtools (Quinlan and Hall, 2010) multiBamCov with the output of bowtie2 (converted to bam, sorted and indexed). R package DESeq (Anders and Huber, 2010) and DESeq2 (Love et al., 2014) were used to identify differentially expressed (DE) circRNAs for individual patient and combining all patients together respectively. |Log2FoldChange|>1 and p-value<0.05 were considered as the cut-off for defining significant DE circRNAs.
For the DE analysis of combination of all patients, we also considered false discovery rate (FDR) at 10%.

Real-Time-Polymerase Chain Reaction and Real-Time Quantitative Polymerase Chain Reaction
According to the manufacturer's protocol, total RNA was extracted using Trizol reagent (251,808, Invitrogen). Reverse transcription and real-time PCR were performed using GoScript ™ Reverse Transcription Mix (A2800 and A6002, Promega). Divergent primers were designed. All primers were synthesized by Sangon Biotech. Results were normalized to the 18S mRNA in each sample. We used the following primer sequences.

CircRNA-miRNA-Gene Interaction Network Analysis
We were interested with those circRNAs which were DE in four or more patients (out of five patients). Then, some top DE circRNAs were selected for miRNA interaction analysis. The circRNA-miRNA interaction was predicted using miRanda (Enright et al., 2003) software. The miRNA sequences were downloaded from mirbase (Kozomara et al., 2019). From the circRNA-miRNA interaction, a sub-network was constructed keeping only the gastric cancer related miRNAs. Then a circRNA-miRNA-gene network was constructed using Cytoscape (Shannon et al., 2003). The genes were the host genes of the circRNAs. Finally, a sub-network was constructed with the top 5 hub circRNAs.

Gene Ontology Term and Kyoto Encyclopedia of Genes and Genomes Pathway Analysis
Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/cgi/input.pl) is a database for obtaining the protein-protein interaction of the provided genes. The host genes of the DE circRNAs were mapped to STRING database, and the protein-protein interaction among the genes was obtained. Then, a PPI network was constructed using the software Cytoscape, and from the network, the top 50 hub genes were selected for GO term and KEGG pathway enrichment analyses. With the selected hub genes, GO term and KEGG pathway enrichment analyses were performed using DAVID (Huang et al., 2009) to know the function of the circRNAs. The threshold p-value<0.05 was used for the significance of the enrichment analysis.

Expression Profiles of circRNAs
A total of 45,783 circRNAs were identified in all samples and 21,652 of them were co-expressed in cancer and normal samples ( Figure 1A). Among the total circRNAs, 79% (36,218/45,783) were exonic, 1% (462/45,783) were intronic, 4% (1848/45,783) were intergenic, 15% (6,744/45,783) were sense overlapping and 1% (511/45,783) were antisense ( Figure 1C). The length distribution of the circRNAs was shown in Figure 1B and observed that the length of 74% (44 + 30) of the circRNAs was less than 1,000 nucleotides (nt). The length of only 13% circRNAs was greater than 2000 nt. The number of back-spliced reads for the 97% of the circRNAs were less than 20 ( Figure 1D). The average number of back spliced reads per circRNA was 6 whereas the minimum and the maximum number of back spliced reads were 2 and 898 respectively. The distribution of circRNAs in different chromosomes were heterogeneous, most circRNAs were originated from chromosome NC_000,001.11 (chromosome 1) ( Figure 1E). The average number of circRNAs per chromosome was 1,405. The length of the circRNAs ranged from 48 to 199,579 nt while median and average lengths were 563 nt and 3,042 nt respectively ( Table 1).
The expression values of the 10 samples from the five patients were shown in Figure 2A. The expression values were almost homogeneous across the samples. We also performed the Principal Component Analysis (PCA) of the expression values and the loading of the PC1 and PC2 were shown in Figure 2B. The distribution of circRNAs across the number of exons was shown in Figure 2C. Majority of the circRNAs were composed of 2 or 3 exons. About 76% of the circRNAs were derived from less than or equal to 5 exons. The distribution of the number of circRNAs per gene was given in Figure 2D. Total 45,783 circRNAs were originated Frontiers in Molecular Biosciences | www.frontiersin.org March 2022 | Volume 9 | Article 857320 from 8,095 genes. Majority of the cases (about 30%), the number of circRNA per gene was 1. About 55% of the cases, the number of circRNAs per gene was less than or equal to 3.

Differential Expression Analysis of circRNAs
We identified the differentially expressed circRNAs between two groups (normal vs cancer). The R package DESeq was used to find the differentially expressed circRNAs. We collected 10 samples from five patients: for each patients two samples were collected one from cancer tissue and other from the adjacent normal tissue. We also performed DE analysis between normal and cancer groups for each patient. The volcano plot was used to find the target circRNAs with the thresholds p-value<0.05 and |log2FoldChange| >1 (Figure 3). The number of DE circRNAs for all samples were given in Table 2.
Out of 45,783 total circRNAs, 478 were differentially expressed between normal vs cancer groups. From Table 2, we observed that about half of the DE circRNAs were upregulated, whereas the other half were downregulated. The expression patterns of these circRNAs were shown in the heatmap (Figure 4). The dendrograms showed that the cancer and normal samples were clearly distinguishable. Firstly, we identified DE circRNAs for all patients (5 normal samples vs 5 cancer samples) and performed DE analysis for each individual patient and isolated those circRNAs which were DE in four or more patients. 4 upregulated and 7 downregulated circRNAs were found to be DE in more than four patients. We selected 8 circRNAs (4 upregulated and 4 downregulated) from the DE analysis of individual patient. These 8 circRNAs were not only DE patient-wise but also DE while performed DE analysis combining all patients (5 normal samples vs 5 cancer samples). Finally, these 8 circRNAs were selected for circRNA-miRNA interaction analysis. The details of the 8 circRNAs were given in Table 3.

CircRNA-miRNA-Gene Interaction Analysis
We got a total of 478 circRNAs as DE between normal vs cancer groups. Among these circRNAs, we investigated those which were DE in more number of patients. We chose those circRNAs for interaction analysis, which were DE in four or more patients (out of five patients).
We selected the 4 upregulated and 4 downregulated (top 4 out of 7 downregulated) circRNAs for miRNA interaction analysis. We got a total of 1990 interactions with 8 circRNAs and 1,278 miRNAs. The miRNA sequences were downloaded from the mirbase database. From the online search we found 64 GC related miRNAs out of 1,278 miRNAs in the network. Next, we extracted the interactions of these 64 miRNAs and obtained a total of 81 interactions. Then, a circRNA-miRNA-gene network was constructed using Cytoscape.    Here, the genes were the host genes of the circRNAs. Finally, we constructed a subnetwork with top 5 hub circRNAs ( Figure 5). From the network, we observed that two genes were functionally related to GC. Two circRNAs circCEACAM5 and circCOL1A1 were interacted with GC related miRNAs and their host genes were also related to GC. This indicated that these two circRNAs might be potential biomarker for GC. The detail list of top 5 hub circRNAs, their interacted GC related miRNAs and their host genes were given in Table 4.

Differential Expression Analysis of the Targeted miRNAs Genes and Candidate circRNAs
We performed the differential expression analysis of the miRNAs in gastric cancer tissue and the adjacent normal tissue using the GEO dataset (GSE158315). R package limma (Ritchie et al., 2015) was used to find the DE miRNAs. The expression patterns of the DE miRNAs were shown in Figure 6A. The dendrograms showed that the cancer and normal samples were clearly distinguishable. The volcano plot ( Figure 6B) was used to find the DE miRNAs with the threshold p-value<0.05 and |log2FoldChange|>1. We found a total of 176 DE miRNAs of which 60 were upregulated and 116 were downregulated. The number of targeted miRNAs for circRNAs circCEACAM5 and circCOL1A1 were 9 and 5 respectively. The unique targeted miRNAs for these two circRNAs were 13 ( Figure 6C). MiRNAs hsa-miR-634 and hsa-miR-429 were common between the 9 targeted miRNAs of circCEACAM5 and the 176 DE miRNAs. Again, miRNAs hsa-miR-634, hsa-miR-21-3p and hsa-miR-145-5p were common between the 5 targeted miRNAs of circCOL1A1 and the 176 DE miRNAs. MiRNAs hsa-miR-429, hsa-miR-21-3p and hsa-miR-145-5p were found to be significantly upregulated in GC ( Figures 6D-F). Therefore, miRNA hsa-miR-429 might be the potential target of the circRNA circCEACAM5, and miRNAs hsa-miR-21-3p and hsa-miR-145-5p might be the potential targets of the circRNA circCOL1A1.
We selected top 8 circRNAs (4 upregulated and 4 downregulated) for the downstream analysis. These 8 circRNAs were originated from 6 host genes. The expression of these 6 host genes were checked by the TCGA data of 408 gastric tumor samples and 36 normal samples obtained from Gepia (version 2) (Tang et al., 2019). We found that 4 genes were upregulated in the GC tissue compared to the normal tissue We also performed the receiver operating characteristic (ROC) curve analysis for our two candidate circRNA biomarkers  circCEACAM5 and circCOL1A1 ( Figures 8A,B). The area under the ROC curves (AUC) for these two circRNAs were significant (p-value<0.05). The AUC values 0.96 and 0.94 represented respectively the accuracy of circCEACAM5 and circCOL1A1 in distinguishing GC and normal patients.

Validation of the circRNAs by qReal-Time-Polymerase Chain Reaction
We selected 8 circRNAs (4 upregulated and 4 downregulated) for miRNA interaction analysis. These 8 circRNAs were DE in four (out of five) or more patients. Among these 8 circRNAs, we selected 5 circRNAs (3 upregulated and 2 downregulated) for experimental validation by qRT-PCR. The expression of these 5 circRNAs were shown in Figure 9. From Figure 9, we observed that the expression of circCEACAM5, circMISP, circCOL1A1 was upregulated in GC, which is consistent with the sequencing result. The expression of circLIPF and circPGC was decreased in GC in line with the sequencing result. Hence, the computational results for these 5 circRNAs were validated by the qRT_PCR experiments.

DISCUSSION
The role of circRNAs in human cancer is unrevealing gradually in recent years. In the present study, we predicted some circRNA biomarkers for gastric cancer through computational approach. We found a total of 45,783 circRNAs from 5 patients (10 Frontiers in Molecular Biosciences | www.frontiersin.org March 2022 | Volume 9 | Article 857320 samples). Then we performed DE analysis and found 4 upregulated and 7 downregulated circRNAs which were DE in four or more patients. Among these, we selected 8 circRNAs (4 upregulated and 4 downregulated) for miRNA interaction analysis.
The most important property of circRNA is known to function as miRNA sponges. As miRNA plays a crucial role in cancer progression, we explored the relationship between circRNAs and miRNAs through circRNA-miRNA-gene interaction network analysis. From the interaction network, we selected 5 hub circRNAs and found that the two circRNAs CircCEACAM5 and CircCOL1A1 were interacted with GC related miRNAs and their host genes were the important therapeutic and prognostic biomarker for GC. The circRNA CircCEACAM5 was interacted with 9 GC related miRNAs and among them miR-634 inhibits the proliferation, migration, and invasion of gastric cancer cell lines , miR-143-3p inhibits GC cell growth and induces apoptosis (Wu et al., 2013), and miR-429 acts as a tumor suppressor in GC cells (Zhang et al., 2016). The Frontiers in Molecular Biosciences | www.frontiersin.org March 2022 | Volume 9 | Article 857320 9 circRNA CircCOL1A1 was interacted with 5 GC related miRNAs for instances, miR-129-5p suppresses the proliferation of GC cells , miR-145-5p is able to inhibit the proliferation, migration and invasion of GC cells (Zhou et al., 2020), miR-501-5p promotes cell proliferation and migration in GC cells (Ma et al., 2020), miR-21-3p is closely related to GC and can be used to predict the prognosis of GC (Sun et al., 2021). We also performed the differential expression analysis of the targeted miRNAs of the circRNAs circCEACAM5 and circCOL1A1. We found that miRNAs miR-429, miR-21-3p and miR-145-5p significantly upregulated in GC compared to normal tissues. Thus, miR-429 might be the potential  target of circCEACAM5 and, miR-21-3p and miR-145-5p might be the potential targets of circCOL1A1. Another important property of circRNA is to regulate gene transcription (Meng et al., 2017). We performed the differential expression analysis of the parental genes of the 8 top up and down regulated circRNAs in GC using TCGA data. We found genes CEACAM5 and COL1A1 as upregulated in GC compared to normal tissues. These two genes were the key therapeutic and prognostic biomarkers for GC (Hu and Chen, 2012;Zhou et al., 2015;Zhaoxing Li et al., 2020). Thus, circCEACAM5 and circCOL1A1 have the potentials to regulate their parental gene transcription.
The GO term analysis of the parental genes of DE circRNAs showed that some host genes were involved in several important biological processes, molecular mechanisms and cellular components. The KEGG pathway analysis showed that some host genes were involved in several significant pathways such as ECM-receptor interaction, Focal adhesion, PI3K-Akt signaling pathway, Cell cycle etc. All of these pathways are involved in important mechanism of GC. ECM-receptor interaction pathway is identified as an important pathway associated with the progression of GC (Hu and Chen, 2012). Focal adhesion plays a key role in regulating cell survival, and proliferation, migration, and invasion of GC cells (Mao et al., 2021). The PI3K-Akt signaling pathway plays an important role in the development and progression of GC (Matsuoka and Yashiro, 2014). And a stronger enrichment of cell cycle pathway is found in GC (Saberi Anvar et al., 2018).
The host genes of these two circRNAs are also functionally related to GC. Gene CEACAM5 is a promising biomarker for prewarning and prognosis of GC (Zhou et al., 2015) and COL1A1 is considered as a potential biomarker for prognosis of GC (Hu and Chen, 2012;Zhaoxing Li et al., 2020). The two circRNAs circCEACAM5 and circCOL1A1 are DE in four or more patients (out of five), are substantially interacted with gastric cancer associated miRNAs, and the host genes of the DE circRNAs are involved in several gastric cancer related pathways. Furthermore the host genes of these two circRNAs are key therapeutic and prognostic biomarker for gastric cancer. The expression of these two circRNAs are also validated by qRT-PCR experiments. Hence, these two circRNAs might be potential biomarker for gastric cancer diagnosis and treatment. Biosynthesis of amino acids 6 0.0489

CONCLUSION
In the current study, we investigated potential circRNA biomarkers for GC through integrated analyses using 10 NGS samples collected from 5 patients. From differential expression and circRNA-miRNAgene interaction analyses, we found two circRNAs circCEACAM5 and circCOL1A1 as DE and these circRNAs were interacted with GC related miRNAs. The host genes of these circRNAs were also the potential therapeutic and prognostic biomarkers for GC. The GO and KEGG pathway analyses revealed that some host genes of the DE circRNAs were enriched in several significant processes (biological, molecular and cellular) and enriched in several GC related pathways. These two circRNAs were also validated by the qRT-PCR experiments. Combining all the results, we can conclude that these two circRNAs circCEACAM5 and circCOL1A1 might play important role in the diagnosis and treatment of GC.

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 below: China National Center for Bioinformation, accession PRJCA008130.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the ethics committee of Shenzhen University School of Medicine. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
TH, SL, YP and YW designed the study. TH and SL performed the bioinformatics pipelines, analysed the data and wrote the manuscript. YP performed the experimental validation. SR, YP, SF, XZ, ZJ and YW extensively edited the manuscript. SF and YW supervised the work. All the authors read and approved the final version of the manuscript for publication.