Original Research ARTICLE
Identification and Analysis of p53-Regulated Enhancers in Hepatic Carcinoma
- 1School of Life Sciences and Engineering, Southwest Jiaotong University, Chengdu, China
- 2Beijing Institute of Biotechnology, Beijing, China
- 3Laboratory Animal Center, Academy of Military Medical Sciences, Beijing, China
- 4Xinxiang Key Laboratory of Pathogenic Microbiology, Xinxiang, China
Enhancers can act as cis-regulatory elements to control transcriptional regulation by recruiting transcription factors (TFs) in a distance and orientation-independent manner. However, it is still unclear how p53 participates in the enhancer network as TF in hepatic carcinoma under the condition of DNA damage. A total of 14,286 active enhancers were identified through the integration of stable and unstable enhancer RNAs (eRNAs) captured by CAGE and GRO-seq, respectively. Furthermore, 218 p53-bound enhancers (Enhp53) were identified by analyzing p53 ChIP-seq in HepG2 cells after DNA damage. The results showed that the enhancer expression and histone markers of enhancers (H3K4me1, H3K4me2, H3K4me3, H3K9ac, and H3K27ac) revealed significantly higher level on Enhp53 than Enhno−p53 which suggested that p53 participated in regulating enhancer activity and chromatin structure. By analyzing 124 TFs ChIP-seq from ENCODE, 93 TFs were found significantly enriched on Enhp53 such as GATA4, YY1, and CTCF, indicating p53 may co-regulate enhancers with TFs participation. Moreover, significantly differentially expressed 438 miRNAs and 1,264 mRNAs were identified by analyzing small RNA-seq and RNA-seq, and 26 Enhp53-miRNAs and 145 Enhp53-mRNA interactions were identified by the integration of 3D genome data and genomic distance. The functional enrichment analysis showed that these miRNA targets and mRNAs were significantly involved in tumor biological processes and signaling pathways such as DNA replication, p53 signaling pathway, hepatitis B, focal adhesion, etc. The above results indicated that p53 participated in regulating enhancer network in hepatic carcinoma and Enhp53 exhibited significantly different characteristics with Enhno−p53.
Enhancers can act as tissue-specific cis-regulatory elements to positively regulate gene expression by recruiting TFs and their cofactors in a distance and orientation-independent manner. Previous studies have shown that most active enhancers can transcribe RNA, namely enhancer RNA (eRNA) under the mediation of transcription factor (TF) (Wang et al., 2011). Studies have shown that the expression of eRNA which is related to the activity of enhancer and dysregulation of the expression of enhancer can lead to the occurrence of various cancers including hepatic carcinoma (Kim et al., 2010).
The structural integrity and stability of DNA are critical for cell survival and physiological functions. DNA can be damaged under various stresses to trigger DNA repair which aims to maintain the integrity of cellular function and immune response (Bartek and Lukas, 2007). The p53 tumor suppressor, as sequence-specific DNA-binding TF, plays a key role in the entire DNA damage repair process (Giono and Manfredi, 2006). A previous study found that most p53-DNA binding sites have consistent signal characteristics with enhancer regions and confirmed that these p53-bound regions had enhancer activity (Melo et al., 2013). Recently, follow-up studies confirmed that a large number of p53 participated in regulating enhancer activity by binding to enhancer regions under doxorubicin-induced DNA damage, suggesting DNA damage may induce p53, which is then involved in enhancer regulatory network (Younger and Rinn, 2017). However, it is still unclear the characteristics and mechanisms of p53-enhancer regulatory network in hepatic carcinoma under the condition of DNA damage.
In this study, we identified a series of activity enhancers through the identification of stable and unstable eRNAs based on cap analysis gene expression (CAGE) and global run-on sequencing (GRO-seq) data in hepatic carcinoma, and two categories were divided according to whether the enhancer bound to p53 or not (Enhp53 and Enhno−p53). The results showed that the enhancer expression and histone markers of enhancers revealed significantly higher level on Enhp53 than Enhno−p53. TF enrichment analysis showed that p53 regulated enhancer activity and chromosome accessibility by directly or indirectly interacting with various TFs and co-factors. Finally, a series of differentially expressed microRNAs (miRNAs) and mRNAs regulated by Enhp53 were identified and the results indicated that they were related to tumorigenesis and the development of hepatic carcinoma significantly.
Materials and Methods
Raw global run-on sequencing (GRO-seq) and p53 ChIP-seq SRA files were downloaded from GEO data set (GSM2428726 and GSE64877). Histone modification, TFs ChIP-Seq and raw Hi-C data (ENCFF419ZIV and ENCFF122SLQ) in HepG2 cells were downloaded from the ENCODE project (GRCh38). Raw miRNA-seq and RNA-seq data under conditions of DNA damage were obtained from our previous study (Yang et al., 2016).
Identification of Active Enhancers in HepG2
Enhancers with stable and unstable transcripts were identified based on extracting CAGE and GRO-seq data. Enhancers that were identified by CAGE method were obtained from the HACER database (http://bioinfo.vanderbilt.edu/AE/HACER/) (Wang et al., 2019). GRO-seq SRA files that were downloaded from the GEO data set were converted into FASTQ format by using sratoolkit (Leinonen et al., 2010). After adapter trimming and reads with a quality score below 10 being removed by Cutadapt (Martin, 2011), reads that longer than 15bp were aligned to the human genome (hg38) by using Bowie2 (Langmead et al., 2009). Reads with mapping quality <10 were removed, and enhancers were identified by using NRSA (http://bioinfo.vanderbilt.edu/NRSA/) (Wang et al., 2018).
Histone Modification, TF and Motif Analysis of Enhancer
p53-bound enhancers were identified when there appeared the intersection of p53 binding sites and the enhancer regions. Histone modification, TF binding were obtained by calculating the signal within 1 kb upstream and downstream of the enhancer. Enriched TF motifs in enhancer were identified by using AME (Bailey et al., 2009) based on known TF motifs which were obtained from the HOCOMOCO database (Kulakovskiy et al., 2016).
Differential Expression Analysis of miRNA and mRNA
Raw sequence reads with low quality were filtered by FastQC (Andrews, 2010). Adapter sequences and low quality reads were removed by using Cutadapt. The trimmed miRNA-seq and RNA-seq reads were mapped to the reference genome by using Bowtie (Langmead, 2010) and HISAT2 (Kim et al., 2015), respectively. The number of reads that were mapped to each gene was counted by using HTSeq (Anders et al., 2015). Differential expression analysis was performed by using edgeR (q-value <0.05 and log2|fold change| ≥ 1).
Hi-C Data Processing and TAD Identification
Raw paired reads were downloaded from ENCODE, then Python code was used to pre-truncate the reads that contain potential ligation junctions. The read pairs were mapped to reference genome using Bowtie2 (Langmead and Salzberg, 2012). Unmapped and low-quality mapped reads with map quality <30 were filtered out and the paired reads are selected using SAMtools (Li et al., 2009) and Unix code. Data normalization, visualization, and TAD analysis were performed by HiCtool (Calandrelli et al., 2018).
Identification of Differentially Expressed miRNAs and mRNAs Regulated by Enhancer
Three methods were used to identify enhancer-regulated miRNAs in our study. FANTOM miRNA transcriptional start site (TSS) annotation was used and miRNA promoters were defined as 1 kb upstream and 500 bp downstream of the TSS of a miRNA. (1) If the enhancer and miRNA promoter were located within the same TAD region, then miRNAs would be considered to be regulated by enhancers. (2) Enhancer-regulated miRNAs were identified if enhancer and miRNA promoter were located in interaction regions based on the 4DGenome (Teng et al., 2015), a general repository for chromatin interaction data. (3) Enhancer-miRNA regulation was identified based on a previous study which provided the linkage score for enhancer-miRNA interaction as the following formula in below, where A and B were the distance from the enhancer center to the closest miRNA gene and the closest gene, respectively. According to the research of Suzuki (Suzuki et al., 2017), miRNA genes with S value below 0.2 were categorized as enhancer-associated miRNAs.
Enhancer-regulated mRNAs were identified according to the following methods, mRNA TSS annotation was obtained from GENCODE database and promoters of mRNA were defined as 1 kb upstream and 500 bp downstream of the TSS of a gene. (1) If the enhancer and mRNA promoter were located within the same TAD region, that mRNAs were considered to be regulated by enhancers. (2) Enhancer-regulated mRNAs were identified if enhancer and mRNA promoter were located in interaction regions based on the 4DGenome. (3) If a mRNA promoter was located within 100 kb upstream or downstream of an enhancer, this mRNA was considered to be regulated by the enhancer (Chepelev et al., 2012). The final enhancer-regulated miRNAs and mRNAs is the union of the results of three methods.
Results and Discussion
Identification of Active Enhancer in HepG2 in DNA Damage
Previous studies have shown that active enhancers produce both stable and unstable transcripts. CAGE is suitable for capturing stable transcripts, while GRO-seq is more suitable for capturing unstable transcripts (Li et al., 2016). Based on these studies, more comprehensive active enhancers were identified by integrating data obtained from these two experimental methods. Firstly, 13,088 active enhancers identified by CAGE in HepG2 were download from the HACER database. Next, 1,321 active enhancers from GRO-seq were identified by obtaining raw data from GEO and NRSA, which were used to identify enhancers (Table S1). Finally, a total of 14,286 enhancers were identified through merging two types of enhancers and converted from hg19 to hg38 by using LiftOver (Table S2). As expected, the active enhancers we identified had known enhancer features, including DNase hypersensitive sites, high levels of H3K27ac, H3K4me1, low levels of H3K4me3, and enrichment of YY1 and CTCF that mediated long-range enhancer-promoter interactions (Figure 1).
Figure 1. An example genomic map of enhancer chr22: 30212552-30213117. Blue represents the histone modification, red marks the TF binding on the chromosome, and the black arc represents a Topology domain.
p53-bound Enhancer (Enhp53) Has Higher Transcriptional Activity Than Enhno-p53
A recent study showed that p53 can bind to enhancers to regulate tumorigenesis and development (Younger and Rinn, 2017). To investigate the p53-bound enhancer in hepatic carcinoma under DNA damage, we matched 14,286 enhancers with 13,723 p53 binding sites obtained from HepG2 cell in the condition of DNA damage. Finally, 218 p53-bound enhancers (Enhp53) were identified (Table S3). To deeply study the effect of p53 on enhancer function, we classified enhancers into two categories based on whether they bound p53 or not: Enhp53 and Enhno−p53. Furthermore, expression levels of these two types of enhancers were compared, and the results showed that the expression level of Enhp53 was higher than Enhno−p53 significantly (Figure 2A), suggesting that the binding of p53 further promoted the enhancer activity.
Figure 2. The differences in expression and histone modification of the two types of enhancers. (A) The expression of the two types of enhancers. p-value is calculated by Wilcoxon test. (B–F) Different histone modification on the two types of enhancers. The blue line represents Enhp53, and the red line marks Enhno−p53.
Previous studies revealed that activity of enhancers are always associated with histone modifications (Creyghton et al., 2010). To investigate if p53 could cause changes of histone modifications surrounding enhancer sequence, we downloaded nine histone modification ChIP-seq data of HepG2 from ENCODE including H3K4me1, H3K4me2, H3K4me3, H3K9ac, H3K9me3, H3K27ac, H3K27me3, H3K36me3, and H4K20me1. It was found that five histone modifications (H3K4me1, H3K4me2, H3K4me3, H3K9ac, and H3K27ac) which had been proven to be markers of identification of active enhancers (Creyghton et al., 2010) showed significant signal peaks on enhancer, and all the signals of these five histone modifications were significantly stronger on Enhp53 than Enhno−p53 (Figures 2B–F). The above results suggested that p53 could affect enhancer activity and might be regulated through changing histone modifications surrounding enhancer sequence.
Transcription-Factors Binding Analysis of p53-bound Enhancer
A previous study has shown that enhancers interact with the gene promoter by recruiting TFs to regulate gene expression (Nolis et al., 2009). To study the effect of TF binding on Enhp53/Enhno−p53 and find which TFs cooperate with p53 as a cofactor to regulate enhancer activity, we compared the TFBS on both types of enhancers. The result showed that 124 TFs could bind to the enhancers by analyzing TF ChIP-seq signals located in the region of ± 1kb from the enhancer center. 93 out of 124 TFs had higher binding signals on Enhp53 than Enhno−p53 (Figure 3A). The other 31 TFs had lower binding signals on Enhp53 than Enhno−p53. Notably, ZNF24, TBP, and JUND had the most significant signal ratio compared with other TFs (Figures 3B–D). Moreover, several TFs associated with enhancer functions, such as GATA4, HNF family, YY1 and CTCF, also presented stronger TF binding signals on Enhp53 compared with Enhno−p53 (Figures 3E–H).
Figure 3. TF binding on Enhp53 and Enhno−p53. (A) Ratios of TF binding signals on the Enhp53 and Enhno−p53 by calculating the mean value of binding signals within 1 kb upstream and downstream of the enhancer. (B–H) TF binding signals located in the region of ± 1 kb from the enhancer center.
Studies have shown that the GATA4 and HNF protein families could promote the formation of open chromatin (Rogerson et al., 2019). Another study showed that p53 could bind to “pioneering sites,” a chromosome inaccessible region in advance and p53 binding sites shifts from inaccessibility to accessibility in response to DNA damage (Younger and Rinn, 2017). In summary, p53 may be involved in the change of chromatin accessibility with the participation of GATA4, HNF4A, and HNF1A in response to DNA damage. Besides, several structural regulators of enhancer-promoter loops such as YY1 and CTCF also revealed higher signal on Enhp53 than Enhno−p53, suggesting p53 may play a key role in mediating long-range chromatin looping.
In order to find more TFs binding to enhancers, we analyzed the TF motifs on both types of enhancers using MEME. Finally, motifs of 200 TFs were obtained on Enhp53 and motifs of 492 TFs were obtained on Enhno−p53 from HOCOMOCO database (Table S4). It was found that p53 and p63 only appeared on Enhp53 compared with other TFs which were present in both two types of enhancers (Figure 4). p63 is an indispensable factor in p53-dependent apoptosis in response to DNA damage (Meek, 2004), suggesting that p53 performs regulation on enhancers with p63's assistance. In summary, p53 may cooperate with other TFs to participate in the regulation of enhancers and chromosome looping.
Figure 4. Motif analysis on p53-bound enhancer. The consensus logo represents the degree of conservation of each position using the height of the consensus character at that position.
Identification of Differentially Expressed miRNAs and mRNAs Regulated by Enhp53
Recent studies have found that enhancers (typical enhancers and super enhancers) regulate miRNA expression and participate in the synthesis of miRNAs regulated by Drosha / DGCR8 (Suzuki et al., 2017). Notably, previous study had proven that p53 also driver the Drosha/DGCR8-mediated primary miRNA processing like enhancer (Hermeking, 2007). These studies indicated that p53, enhancers, miRNAs, and mRNAs may form a complex network and play a critical role in tumorigenesis. To investigate this, we reanalyzed the miRNA-seq and RNA-seq data from previous studies in our laboratory in response to DNA damage (p53 activation) in HepG2 cells. Finally, significantly differentially expressed 438 miRNAs and 1,264 mRNAs were identified after activating p53 in HepG2. Among them, 76 miRNAs were significantly down-regulated and 362 miRNAs were significantly up-regulated. Meanwhile, 724 mRNAs were significantly up-regulated and 540 mRNAs were significantly down-regulated.
Furthermore, in order to identify the miRNA regulated by Enhp53 as comprehensively as possible, we identified the miRNA regulated by the enhancer through two regulation modes: distal and proximal enhancer-miRNA regulation. Finally, 26 Enhp53-miRNAs interactions were identified by the integration of three approaches (TAD, 4DGenome, and proximal distance, details see methods) (Table 1). Several miRNAs have shown to be post-transcriptionally upregulated in a p53-dependent and p68/p72-dependent manner with DNA damage such as miR-16-5p (upregulated 8.63-fold), miR-143-3p (upregulated 8.46-fold), and overexpression of these p53-induced miRNAs caused a decreasing rate of cell proliferation (Wan et al., 2011).
Similar to the identification of Enhp53-miRNA, 145 Enhp53-mRNA pairs of 121 differentially expressed mRNAs were identified by the integration of chromosome interaction and distance (See methods, Table S5). Notably, there were five enhancers that regulated both miRNA and mRNA, involving a total of 7 miRNAs and 11 mRNAs (Table S6), suggesting that these enhancers might play a key role in the p53 network.
To analyze the function and the potential pathway of miRNA target genes and mRNAs regulated by Enhp53, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway functional enrichment analysis were performed by using the miRTarBase database (Hsu et al., 2011) and the R package clusterProfiler (Yu et al., 2012). GO and KEGG pathway analysis revealed that these miRNA targets and mRNAs were significantly involved in tumor biological processes such as DNA replication, folding, and chromatin assembly (FDR <0.05). KEGG pathway enrichment analysis showed that miRNA target genes and mRNAs regulated by Enhp53 significantly enriched in p53 or hepatic carcinoma related pathways such as p53 signaling pathway, focal adhesion, mTOR signaling pathway and hepatitis B (FDR <0.05, Figure S1). In addition, in order to study the effect of Enhp53-regulated miRNAs target genes and mRNAs on hepatocarcinogenesis, we obtained the driver genes affecting hepatocarcinogenesis from the driverDB database (Liu et al., 2020). The results showed that 32% of the driver genes are miRNA target genes and mRNA regulated by Enhp53, such as TP53, AGO2, EGFR, FOXP1, etc. In summary, the miRNA and mRNA regulated by Enhp53 are critical for the occurrence of liver cancer.
In this study, a total of 218 p53-bound enhancers were identified by analyzing CAGE and GRO-seq data. It was found that Enhp53 had higher transcriptional activity and ability to bind TFs, and Enhp53 was also significantly different from Enhno−p53 in TF motif recognition and histone modifications. To further identify miRNAs and mRNAs regulated by Enhp53, we integrated multiple methods to identify 26 Enhp53-miRNAs and 145 Enhp53-mRNAs based on distal and proximal enhancers regulation. The results revealed that these miRNA targets and mRNAs were significantly involved in tumor biological processes and signaling pathways. However, these findings still required further finer experimental verification to prove this network. In summary, the results above indicated that p53 was involved in the hepatic tumorigenesis and development by mediating enhancers under DNA damage and provided a theoretical method basis for exploring the regulation of TFs and enhancers in the future studies.
Data Availability Statement
All datasets generated for this study are included in the article/Supplementary Material.
ZG and JY conceived and designed the experiments. QH and WW acquired the experiment data. YZ, MQ, and FT performed the study and carried out the data analysis. ZG, YL, and ZL wrote this manuscript. ZG, ZQ, and BL revised the manuscript. All authors have read and approved the final manuscript.
This work was supported by the National Science and Technology Major Project of Infectious Diseases (2018ZX10101-003-001-008) and National Natural Science Foundation of China (31671363).
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2020.00668/full#supplementary-material
Bailey, T. L., Boden, M., Buske, F. A., Frith, M., Grant, C. E., Clementi, L., et al. (2009). MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res. 37, W202–W208. doi: 10.1093/nar/gkp335
Calandrelli, R., Wu, Q., Guan, J., and Zhong, S. (2018). GITAR: an open source tool for analysis and visualization of Hi-C data. Genomics Proteomics Bioinform. 16, 365–372. doi: 10.1016/j.gpb.2018.06.006
Chepelev, I., Wei, G., Wangsa, D., Tang, Q., and Zhao, K. (2012). Characterization of genome-wide enhancer-promoter interactions reveals co-expression of interacting genes and modes of higher order chromatin organization. Cell Res. 22, 490–503. doi: 10.1038/cr.2012.15
Creyghton, M. P., Cheng, A. W., Welstead, G. G., Kooistra, T., Carey, B. W., Steine, E. J., et al. (2010). Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc. Natl. Acad. Sci. U.S.A. 107, 21931–21936. doi: 10.1073/pnas.1016071107
Hsu, S. D., Lin, F. M., Wu, W. Y., Liang, C., Huang, W. C., Chan, W. L., et al. (2011). miRTarBase: a database curates experimentally validated microRNA–target interactions. Nucleic Acids Res. 39, D163–D169. doi: 10.1093/nar/gkq1107
Kim, T. K., Hemberg, M., Gray, J. M., Costa, A. M., Bear, D. M., Wu, J., et al. (2010). Widespread transcription at neuronal activity-regulated enhancers. Nature 465, 182–187. doi: 10.1038/nature09033
Kulakovskiy, I. V., Vorontsov, I. E., Yevshin, I. S., Soboleva, A. V., Kasianov, A. S., Ashoor, H., et al. (2016). HOCOMOCO: expansion and enhancement of the collection of transcription factor binding sites models. Nucleic Acids Res. 44, D116–D125. doi: 10.1093/nar/gkv1249
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25
Liu, S. H., Shen, P. C., Chen, C. Y., Hsu, A. N., Cho, Y. C., Lai, Y. L., et al. (2020). DriverDBv3: a multi-omics database for cancer driver gene research. Nucleic Acids Res. 48, D863–D870. doi: 10.1093/nar/gkz964
Melo, C. A., Drost, J., Wijchers, P. J., van de Werken, H., de Wit, E., Vrielink, J. A. O., et al. (2013). eRNAs are required for p53-dependent enhancer activity and gene transcription. Mol. Cell 49, 524–535. doi: 10.1016/j.molcel.2012.11.021
Nolis, I. K., McKay, D. J., Mantouvalou, E., Lomvardas, S., Merika, M., Thanos, D., et al. (2009). Transcription factors mediate long-range enhancer–promoter interactions. Proc. Natl. Acad. Sci. U.S.A. 106, 20222–20227. doi: 10.1073/pnas.0902454106
Rogerson, C., Britton, E., Withey, S., Hanley, N., Ang, Y. S., Sharrocks, A. D., et al. (2019). Identification of a primitive intestinal transcription factor network shared between esophageal adenocarcinoma and its precancerous precursor state. Genome Res. 29, 723–736. doi: 10.1101/gr.243345.118
Suzuki, H. I., Young, R. A., and Sharp, P. A. (2017). Super-enhancer-mediated RNA processing revealed by integrative microRNA network analysis. Cell 168, 1000–1014.e1015. doi: 10.1016/j.cell.2017.02.015
Wang, D., Garcia-Bassets, I., Benner, C., Li, W., Su, X., Zhou, Y., et al. (2011). Reprogramming transcription by distinct classes of enhancers functionally defined by eRNA. Nature 474, 390–394. doi: 10.1038/nature10006
Wang, J., Dai, X., Berry, L. D., Cogan, J. D., Liu, Q., Shyr, Y., et al. (2019). HACER: an atlas of human active enhancers to interpret regulatory variants. Nucleic Acids Res. 47, D106–D112. doi: 10.1093/nar/gky864
Wang, J., Zhao, Y., Zhou, X., Hiebert, S. W., Liu, Q., Shyr, Y., et al. (2018). Nascent RNA sequencing analysis provides insights into enhancer-mediated gene regulation. BMC Genomics 19:633. doi: 10.1186/s12864-018-5016-z
Yang, Y., Liu, W., Ding, R., Xiong, L., Dou, R., Zhang, Y., et al. (2016). Comprehensive expression profiling and functional network analysis of p53-regulated microRNAs in HepG2 cells treated with doxorubicin. PLoS ONE 11:e0149227. doi: 10.1371/journal.pone.0149227
Keywords: p53, enhancer, hepatic carcinoma, transcription factors, microRNA
Citation: Zhang Y, Qian M, Tang F, Huang Q, Wang W, Li Y, Li Z, Li B, Qiu Z, Yue J and Guo Z (2020) Identification and Analysis of p53-Regulated Enhancers in Hepatic Carcinoma. Front. Bioeng. Biotechnol. 8:668. doi: 10.3389/fbioe.2020.00668
Received: 16 March 2020; Accepted: 28 May 2020;
Published: 30 June 2020.
Edited by:Feng-Biao Guo, University of Electronic Science and Technology of China, China
Reviewed by:Dong Li, Beijing Proteome Research Center, China
Xiaobin Zheng, Carnegie Institution for Science (CIS), United States
Sujun Zheng, Capital Medical University, China
Copyright © 2020 Zhang, Qian, Tang, Huang, Wang, Li, Li, Li, Qiu, Yue and Guo. 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.
†These authors have contributed equally to this work and share first authorship