Background: The purpose of this study was to investigate the relationship between Long non-coding RNAs (lncRNAs) expression and epigenetic changes, and explore the prognostic value of these findings in patients with prostate cancer.
Methods: In this study, we systematically compared the histone modification and methylation regions on lncRNAs promoter and enhancer elements by multiomics. We analyzed the distribution characteristics of histone modified apparent lncRNAs promoters and enhancers in the genome. Single sample Gene Set Enrichment Analysis (ssGSEA) and Kyoto Encyclopedia of Gene and Genomes (KEGG) were used to analyze the functional enrichment of epigenetic dysregulated lncRNAs. The prostate cancer associated lncRNAs were downloaded from lnc2cancer V3.0 and survival analysis was performed.
Results: We found 1124 epigenetic-dysregulated lncRNAs (epi-lncRNAs) and 9734 epigenetic-dysregulated protein coding genes (epi-PCGs) in prostate cancer, and the abnormal frequency of lncRNAs was much lower than that of PCGs. H3K4me3, H3K4me1, H3K27ac and H3K27me3 abnormally modified histones accounted for a large proportion and were mainly concentrated in the promoter region. Epigenetic abnormal lncRNAs affects the biological specific molecular function of prostate cancer by changing the abnormal modification of histone. The abnormal expression of lncRNAs caused by abnormal modification of H3K36me3 may be an important factor in the occurrence of prostate cancer. Finally, based on lncRNAs analysis of the prognosis of prostate cancer samples, three lncRNAs (HAR1A, SNHG12 and SNHG17) were identified as prognostic markers of prostate cancer.
Conclusions: These findings may help to better understand the abnormal epigenetic regulation of lncRNAs expression in prostate cancer.
Prostate cancer (PCa) remains the second leading cause of deaths due to cancer in the United States in men . Siegel et al.  demonstrated that there will have 164,690 newly diagnosed PCa patients and 29,430 deaths in 2018 in the United States. Thus, the diagnosis of PCa in early stage is vitally important . Currently, the serum prostate-specific antigen (PSA) screening remains the primary way for early diagnosis of PCa. However, the sensitivity and specificity of PSA test remains low . Therefore, it is important to find notable biomarkers for the diagnosis of PCa.
Previous studies have shown that epigenetic modifications plays a crucial role in the diagnosis and prognosis of PCa and proving additional options for PCa diagnosis and treatment strategies [5, 6]. Distinct PCa molecular subtypes have been characterized based on gene expression, clinical, and prognostic features of PCa. Long non-coding RNA (lncRNA) has been linked to many human cancers, including genetically heterogeneous PCa [7, 8]. lncRNAs are frequently differentially expressed in PCa subtypes . Aberrant lncRNA expression could be caused by alteration of epigenetic patterns, such as changes in DNA methylation and post-translational histone modifications. For example, lncRNA DLX6-AS1 promotes malignant phenotype and lymph node metastasis in PCa by inducing LARGE methylation . Long non-coding RNA LINC00673 silencing inhibits proliferation and drug resistance of PCa cells via decreasing KLF4 promoter methylation . LncRNA HOXD-AS1 recruited WDR5 to directly regulate the expression of target genes by mediating histone H3 lysine 4 tri-methylation (H3K4me3), Thus promotes proliferation, castration resistance, and chemo-resistance in PCa . Although global epigenetic abnormalities have been identified as prominent cancer hallmarks, attempts to characterize the relationship between epigenetic alterations and expression of lncRNAs and their prognostic value in prostate cancer subtypes have been limited.
In this study, we screened out the peaks specific to prostate cancer based on the physical location of the histone modified peaks, and identification of lncRNAs and protein coding genes (PCGs) with epigenetic disorders. Subsequently, we extracted the expression profiles of lncRNAs caused by different histone modifications, and calculated the enrichment scores of each sample in these lncRNAs using single sample Gene Set Enrichment Analysis (ssGSEA). Collected prostate cancer-related lncRNAs gene were further loaded from Lnc2Cancer v3.0, and were performed survival analysis. Finally, three epigenetic-dysregulated lncRNA (epi-lncRNA) genes specific for prostate cancer were screened out, and a prognostic-related model was established.
The gene expression profile of prostate cancer gene expression profile and its corresponding clinical information were downloaded from the Cancer Genome Atlas (TCGA) database (https://confluence.broadinstitute.org/display/GDAC/Home/). Through the gene annotation file of GENCODE, the expression profiles were divided into lncRNAs and PCGs. Based on the sample number of TCGA, the expression profiles of lncRNAs and PCGs were divided into 495 prostate cancer samples and 52 normal samples. The 450K array data of prostate adenocarcinoma (PRAD) were downloaded from TCGA database. The array data was split into 498 prostate cancer samples and 50 normal samples based on the PRAD sample number. The K-Nearest Neighbors (KNN) method was used to fill in missing values in the prostate cancer data. The replicated narrowPeak data of 6 types of histone modifications (H3K9me3, H3K36me3, H3K27me3, H3K27ac, H3K4me1, and H3K4me3) was downloaded from the hg38 version of the prostate cancer cell line PC-3 (https://www.ncbi.nlm.nih.gov/assembly/GCF_000001405.26/) and normal cell line RWPE2 in the Encyclopedia of DNA Elements (ENCODE) database (https://www.encodeproject.org/). The data of 218 prostate cancer samples with prognostic information were download from the cBioPortal database (https://www.cbioportal.org/, MSKCC, Nat Med 2017).
To explore the correlation between epigenetic variation and prostate cancer, we converted the normalized gene expression profile count into integer Count, and then used the R library “DESeq2” to identify differentially expressed lncRNAs and PCGs, The Benjami-hochberg method was used to determine P value. IncRNAs or PCGs with false discovery rate (FDR) <0.05 were considered significant. We screened out the peaks specific to prostate cancer based on the physical location of the histone modified peaks, and only retained the peaks with p-value <0.05 as the differential peaks, and then we combined the GENCODE file to obtain the differential histone modifications genes. We downloaded the human enhancer database from FANTOM5(https://fantom.gsc.riken.jp/5/) to identify the enhancer of the gene; the promoter of the gene is defined as 2 kb upstream of the transcription start site (TSS) and 0.5 kb downstream. We used the R library “ChIPseeker”  and the human lncRNAs Transcripts database(https://lncipedia.org/) to identify gene promoters. Finally, we used the Bumphunter method of the R library ‘CHAMP’ [16, 17] to identify the differentially methylated region (DMR). The area where Bumphunter DMR p.value <0.001 is considered a significant DMR. Finally, we defined lncRNAs and PCGs as epigenetic abnormally regulated lncRNAs or PCGs. The criteria were as follows: (1) lncRNAs and PCGs are differentially expressed in prostate cancer and normal samples; (2) promoters or enhancers at least overlap a region of differential histone modification or differential methylation (named epi-lncRNA, non-epi-lncRNA, epi-PCG, non-epi-PCG). To compare the genomic characteristics of epigenetic and non-dysregulated lncRNA/PCG, we compared the exons and transcripts of epi-lncRNA, non-epi-lncRNA, epi-PCG, and non-epi-PCG genes, and the number and length of genes.
To explore the epigenetic characteristics of lncRNAs caused by histone modifications, we analyzed the distribution characteristics of the promoters and enhancers of epi-lncRNAs with different histone modifications on the genome.
The "clusterProfiler" package in R software was used to implement the Kyoto Encyclopedia of Gene and Genomes (KEGG) pathway . The enrichment score of each sample in these lncRNAs was calculated by ssGSEA algorithm in the "GSVA"  package.
We extracted the m6a, m5c, and m1A gene expression profiles from the TCGA prostate cancer expression profile, and calculated the correlation between these genes and 12 epi-lncRNA enrichment score.
To understand the relationship between epigenetic disorders lncRNA and cancer, we used Lnc2Cancer v3.0 (http://www.bio-bigdata.net/lnc2cancer/) to analyze the enrichment of epi-lncRNA and known cancer lncRNAs. We screened out epi-lncRNA related to prostate cancer. In addition, we combined the transcripts per kilobase million (TPM) expression profile data of TCGA to identify epi-lncRNAs with significant changes in expression in prostate cancer.
To understand the potential prognostic value of epi-lncRNA dysregulation, we divided the samples into two subgroups of high expression and low expression based on the expression value of the above epi-lncRNA in cancer samples. The Kaplan Meier survival curve was used to estimate the survival difference between the two groups by “survminer” R package. The Receiver Operating Characteristic (ROC) curve was applied to evaluate the predictive accuracy of the risk model by the “timeROC” package in R . The R software package max_stat was used to truncate genes to obtain high and low expression groups.
We aim to analyze the association between differential expression of lncRNA and epigenetic alterations in prostate cancer. Deseq2 identified significantly differentially expressed genes, including 4836 lncRNAs (Additional file1: Table S1) and 12410 PCGs (Additional file2: Table S2). Binding histone modification data and 450k methylation chip data, we identified 1124 epi-lncRNAs, 12130 non-epi-lncRNAs, 9734 epi-PCGs and 9626 non-epi-PCGs (Additional file3: Table S3). The aberrant frequency of lncRNAs (8%) was much lower than that of PCGs (50%) in prostate cancer (Fig. 1A).
To characterize genomic signatures of epi-lncRNAs, we compared the number and length of genes, exons of epi-lncRNAs, other lncRNAs (non-epi-lncRNAs), epi-PCGs, and other PCGs (non-epi-PCGs). Epi-lncRNAs had a larger number of transcripts but the transcript length was lower than that of non-epi-lncRNAs (Fig. 1B, C); epi-PCGs had a larger number of transcripts and the transcript length was higher than that of non-epi-PCGs; at the same time, epi- lncRNAs have longer exon length, epi-PCGs has more exon number and exon length (Fig. 1D, E).
We systematically analyzed epi-lncRNAs in prostate cancer, revealing the landscape of epi-lncRNAs caused by different histone modifications and differentially methylated regions (Fig. 2A). Abnormal histones modifications in lncRNA were mainly H3K4me3, H3K4me1, H3K27ac, and H3K27me3, and most epi-lncRNAs were accompanied by a variety of abnormal histones modifications (Fig. 2B). We also observed that the abnormally modified regions of these histones were mainly concentrated in the promoter region (Fig. 2B).
To characterize the potential function of lncRNAs dysregulation caused by these histone modifications, we systematically analyzed the relationship between epi-lncRNAs expression and pathways in prostate cancer. The enrichment scores in adjacent samples were significantly higher than those in tumor samples, including H3K27ac_promoter, H3K4me1_enhancer, H3K4me3_enhancer, H3K4me3_promoter, H3K9me3_promoter, H3K27me3_enhancer, and H3K27me3_promoter (Fig. 3A). In addition, the enrichment scores of H3K9me3_enhancer, H3K36me3_enhancer, and H3K36me3_promoter in tumor samples were significantly higher than those in adjacent samples (Fig. 3A). This indicates that H3K9me3_enhancer, H3K36me3_enhancer, and H3K36me3_promoter might have cancer-promoting effects, while H3K27ac_promoter, H3K4me1_enhancer, 3K4me3_enhancer, H3K4me3_promoter, H3K9me3_cerpromoter, H3K9me3_promoter, and H3K9me3_promoter might have protective effects.
In addition, we also evaluated the KEGG pathway score of each sample, and calculated the relationship between the enrichment score of each epi-lncRNA and KEGG pathway. A total of 39 KEGG pathways were related to 12 types of epi-lncRNA (Fig. 3B). This suggested that different types of epi-lncRNA-related pathways have a certain consistency, and were closely related to the occurrence, development and metabolism of tumors.
RNA modification is an important epigenetic feature, which is related to a variety of important biological processes. We have observed that these enrichment scores were significantly correlated with M6a, M5c and M1A modified genes (Fig. 4). Meantime, we found significant differences in the correlation between H3K36me3 _ enhancer, H3K36me3 _ promoter and M6A, M5c, M1A modified genes. The former was significantly correlated with M6A RNA modification related genes.
To further ensure the regulatory role of candidate lncRNAs in prostate cancer, we collected prostate cancer-related lncRNA gene set from Lnc2Cancer v3.0, and found that 14 epi-lncRNAs have been reported to be directly related to prostate cancer (Fig. 5A; Additional file4,5: Table S4,5).
Then we calculated the difference between 14 epi-lncRNAs in prostate cancer samples and normal samples. The expression of HOTTIP, LINC00663, LINC00844, and MIR222HG in tumors was significantly lower than that of normal samples, while the expression of the remaining 10 genes (GLIDR, HOXA11−AS, LINC00308, PCAT14, PCAT18, PCAT5, SNHG12, SNHG17, SNHG6, and TMPO-AS1) in tumors was significantly higher than that of normal samples (Fig. 5B). Among these 14 epi-lncRNAs, only HOXA11-AS was related to both histone-modified promoters and DMR promoters, and the remaining 13 epi-lncRNAs were related to histone-modified promoters (Fig. 5C, D; Additional file6: Fig S1).
Overall, the 78.57% of these 14 genes on the TCGA dataset are consistent with the performance of the Lnc2Cancer v3.0 database (Additional file5: Table S5). Among them, the expression of 11 genes (GLIDR, HOXA11−AS, LINC00308, LINC00663, LINC00844, PCAT14, PCAT18, PCAT5, SNHG17, SNHG6, and TMPO-AS1) in the tumors of the TCGA data was consistent with the performance in the Lnc2Cancer v3.0 database.
To gain insight into the potential prognostic value of epi-lncRNAs dysregulation, we used 91 epi-lncRNA genes in the above results to perform survival analysis on TCGA prostate cancer samples. In the univariate COX analysis, there were four genes (HAR1A, LINC00629, SNHG12 and SNHG17) related to the survival of prostate cancer patients (Fig. 6A). We found that except for the high and low expression survival curves of gene LINC00629, the survival curves of the other 3 genes were significant (Fig. 6B-E).
To establish a prognostic model of epi-lncRNAs dysregulation, the 3-EpiLncRNA signature was established using multivariate survival analysis in the TCGA data set. The higher 3-EpiLncRNA score has a higher mortality rate (Fig. 7A). ROC analysis showed that the 3-EpiLncRNA score of 1-, 3-,5-year AUC was 1, 0.78,0.79, which had good prognostic performance (Fig. 7B). The prognosis of high 3-EpiLncRNA score was significantly worse than that of low 3-EpiLncRNA score (Fig. 7C).
To evaluate the portability of 3-EpiLncRNA score, we obtained the data set MSKCC with prognostic information from the cBioPortal database, and extracted the expression profiles of 3 epi-lncRNAs (the gene SNHG12 does not exist in the data set MSKCC). The 3-EpiLncRNA score of the sample was evaluated by the above method, and ROC analysis showed that the 1-year AUC of 3-EpiLncRNA score was 0.69, which had good prognostic prediction performance (Fig. 7D). The prognosis of high 3-EpiLncRNA score was significantly worse than that of low 3-EpiLncRNA score (Fig. 7E), which was consistent with the previous results.
In this study, we found a large number of epi-lncRNAs by using public databases such as TCGA, ENCODE and lncRNAs Transcripts, combined with the joint analysis of genomics, transcriptomics and epigenetic genomics. Many of these epi-lncRNAs have been reported to be associated with cancer in previous literature.
It can be clearly seen from our results that the rate of epi-lncRNA in prostate cancer (8.5%) is much lower than that of epi-PCG (50.2%) (Fig. 1), which is basically consistent with previous literature reports. Previous studies have shown that in specific tissues and cells, the overall occupancy of histone marks across lncRNA genes is in the range 27% - 38%, while that on the transcriptional start sites of protein-coding is 65%-73% . LncRNA transcription is generally regarded as a genome-wide monitoring mechanism and plays an important role in RNA quality control. The relatively low aberration rate of LncRNA in prostate cancer also indirectly reflects its high stability in the process of gene transcription or expression, which may be more suitable as a molecular marker of prostate cancer. In addition, from characterizing genomic signatures, it is also observed that epi-lncRNA had similar structural characteristics to epi-PCG except for the shorter length of transcripts. Compared with non-epi-lncRNA and non-epi-PCG, they have more transcripts, exons and longer exon length (Fig. 2). Based on the multiple relationships between epi-lncRNAs and PCGs (inter-gene, overlap, partial overlap, intron or exon) [23, 24], the more complex the splicing pattern of lncRNA is, the more likely it is to be regulated by abnormal epigenetic modification . Our results also coincide with the above view.
In the landscape of epi-lncRNA, it is not difficult to find that H3K4me3, H3K4me1, H3K27ac and H3K27me3 account for a high proportion of abnormally modified histones in the causes of lncRNA epigenetic disorders, and these abnormally modified histones are mainly concentrated in the promoter region (Fig. 2). Combined with the ssGSEA analysis of epi-lncRNAs (Fig. 3), we further found that the abnormal expression of H3K4me3, H3K4me1, H3K27ac and H3K27me3, whether located in the enhancer or promoter region, they almost consistently showed protective effects. H3K9me3_enhancer, H3K36me3_enhancer and H3K36me3_promoter showed significant cancer promoting effect. Although the latter accounts for only a small proportion in the epi-lncRNA genome landscape, it eventually leads to tumorigenesis. It shows that abnormal modification of H3K36me3 plays an extremely important role in the development of prostate cancer. This seems to reflect a message that under the gene expression pattern of prostate cancer, the final outcome of epigenetic aberrant regulation does not depend on the number of histone aberrant modifications. On the contrary, the type, location or function of histone aberrant modification may be more important. Inhibition of abnormal modification of H3K36me3 may bring new hope for the treatment of prostate cancer.
We analyzed 39 KEGG pathways most related to 12 kinds of epi-lncRNA and found that different types of epi-lncRNA-related pathways had certain consistency and included a large number of metabolism-related pathways (HISTIDINE_METABOLISM, TYROSINE_METABOLISM, BETA_ALANINE_METABOLISM, PHENYLALANINE_METABOLISM) and tumor-related pathways (COLORECTAL_CANCER, ENDOMETRIAL_CANCER, PROSTATE_CANCER). Thus, it can be seen that epigenetic abnormal LncRNA affects the biological specific molecular function of prostate cancer by changing the abnormal modification of histone and plays an important role in the pathogenesis of prostate cancer.
Previously, RNA modification was not considered to be a cancer driver, but accumulated evidence gradually suggests that abnormal RNA modification contributes to cancer cell proliferation, self-renewal, migration, stress adaptation and survival, and was named epitranscriptome in 2015. The functional study of these modifications is now on the rise and has shown great significance in human pathology [27, 28]. We visualized the correlation between 12 categories of epi-lncRNA enrichment scores and M6A, M5C and M1A modified genes. Among them, enhancer and promoter of H3K9me3, H3K27me3, H3K27ac, H3K4me1 and H3K4me3 have both similar and unique correlations with these genes. This suggests that there may be different regulation modes of lncRNA disorders caused by histone modification in enhancer and promoter, and these epi-lncRNAs are closely related to RNA modification. These results have confirmed previous findings. These findings [29–31] emphasize that DNA, RNA and histone and their modification work together in a synergistic manner, leading to a special chromatin state that determines the important function and ultimate direction of the genome. m6A is an important post-transcriptional regulation mechanism of genes [32–34]. It is also the most abundant and well-characterized internal modification in mRNA. It plays an important role in a variety of normal and pathobiological processes, including the regulation of self-renewal of embryonic stem cells and cancer cells, and survival after DNA damage [35–37] However, how N6-Methyladenosine m6A is accurately and dynamically deposited in the transcriptome has always been a mystery. Until a recent study found , H3K36me3 was regarded as a transcriptional extension marker that could guide m6A modification globally. M6A modification was enriched near the H3K36me3 peak. When intracellular H3K36me3 was exhausted, m6A modification decreased as a whole. H3K36me3 can be directly recognized and bound by RNA methyltransferase METTL14, promoting the binding of m6A-MTC to adjacent RNA polymerase II, thus transferring m6A-MTC to newly born RNA with active transcriptional activity, resulting in m6A co-transcriptional deposition. Our results further support the above findings.
One of the 14 prostate cancer epi-lncRNAs screened (HOXA11-AS) was not only related to the histone modified promoter but also related to the DMR promoter. The remaining 13 related epi-lncRNAs (GLIDR, LINC00308, PCAT14, PCAT18, PCAT5, SNHG12, SNHG17, SNHG6, MPo -AS1, HOTTIP, LINC00663, LINC00844, MIR222HG) was only associated with histone modified promoters. These results were consistent with our previous findings in the landscape of epigenetic-dysregulated lncRNAs of prostate cancer, that is, the areas with abnormal histone modification were mainly concentrated in the promoter region (Fig. 3B). By comparing the performance of the 14 genes on the TCGA dataset, 78.57% of them were consistent with those in the Lnc2Cancerv3.0 database, which once again verified the reliability of the research results.
We identified three lncRNAs (HAR1A, SNHG12, SNHG17) to construct the prostate cancer 3-EpiLncRNA signature. These three lncRNAs are independent prognostic markers of prostate cancer and have been identified as risk factors for a variety of cancers. HAR1A, also known as LINC00064, is widely expressed in brain, testis, spleen, lymph nodes and prostate tissues. Previous studies have shown that HAR1 is related to the development and evolution of the brain [39, 40]. A recent study suggests that HAR1 is also associated with tumorigenesis, regulating the progression of oral cancer through α-kinase-1, bromodomain-7 and myosin IIA axis.SNHG12 is usually overexpressed in tumor cells and may promote tumorigenesis and metastasis by acting as a sponge of microRNA[42, 43]. There is evidence that SNHG12 promotes prostate tumorigenesis and progression through AKT regulation , and promotes renal cell carcinoma progression and resistance to sunitinib through up-regulation of CDCA3. LncRNASNHG17 is a non-coding RNA widely expressed in human body, with the highest expression in bone marrow, followed by high expression in ovary, prostate, thyroid and other organs, which is closely related to the progression of osteosarcoma  and promotes the proliferation of gastric and colon cancer cells through epigenetic silencing of P57[47–49]. In order to evaluate the predictive value of the model, we visualize the model and draw the survival curve and ROC curve (Fig. 8). The high 3-EpiLncRNA score of prostate cancer began to show a significant survival disadvantage in the fifth year, and the number of survivors in the 10th year was only 1/3 of the low 3-EpiLncRNA score. ROC analysis showed that the 1-year AUC, 3-year AUC and 5-year AUC of 3-EpiLncRNA signature were 1, 0.78 and 0.79, respectively, indicating that the model has a good ability to evaluate the prognosis and can be used to predict the survival of prostate cancer.
In summary, our study analyzed the relationship between abnormal expression of lncRNA and epigenetic disorders and related genomic characteristics, and screened three specific epi-lncRNA for the construction of prostate cancer prognosis-related models. Mining data from a new perspective and supporting previous experimental research. At the same time, it is also found that the abnormal expression of LncRNA caused by the abnormal modification of enhancer-related H3K36me3 may be an important factor in the occurrence of prostate cancer. Reversing the abnormal modification of H3K36me3 may provide a new way for the treatment of prostate cancer, but further experimental research is needed.
At present, unexpected progress has been made in the analysis of the molecular mechanism of epigenetic regulation, which has a far-reaching impact on a better understanding of normal development and the treatment of human diseases. Different from genetic changes, the errors in epigenetic characteristics will be reversible, which has important guiding significance for disease treatment and drug research. With the deepening of epigenetic and other multi-group studies, the veil of prostate cancer will be lifted layer by layer, and it is expected to achieve an important breakthrough in future research.
lncRNA: Long non-coding RNA
ssGSEA: Single sample Gene Set Enrichment Analysis
KEGG: Kyoto Encyclopedia of Gene and Genomes
epi-lncRNA: epigenetic-dysregulated lncRNA
epi-PCGs: epigenetic-dysregulated protein coding genes
PCa: Prostate cancer
PSA: prostate-specific antigen
H3K4me3: H3 lysine 4 tri-methylation
PCGs: protein coding genes
TCGA: the Cancer Genome Atlas
PRAD: prostate adenocarcinoma
KNN: K-Nearest Neighbors
ENCODE: Encyclopedia of DNA Elements
FDR: false discovery rate
TSS: transcription start site
DMR: differentially methylated region
TPM: transcripts per kilobase million
ROC: Receiver Operating Characteristic
Ethics approval and consent to participate
Consent for publication
Availability of data and materials
The datasets generated and analysed during the current study are available in the Cancer Genome Atlas (TCGA) database [https://confluence.broadinstitute.org/display/GDAC/Home/]and the cBioPortal database [https://www.cbioportal.org/, MSKCC, Nat Med 2017].
The authors declare that they have no competing interests.
This work was supported by the ZheJiang Province Public Welfare Technology Application Research Project [number 2021KY782] and the Wenzhou Municipal Science and Technology Bureau [number Y2020151].
Conceptualization, J.Z, X.X, H.S; methodology and software, J.Z, Z.X, Y.X, X,H; validation, J.Z, Y.X, X.H, H.S; formal analysis, X.X, Z.X, Y.X, X.H; data curation, all authors; writing—original draft preparation, J.Z, X.X, H.S; writing—review and editing, J.Z, X.X, H.S; visualization, J.Z, Z.X; supervision, H.S; project administration and funding acquisition. All authors have read and agreed to the published version of the manuscript.
Figure 8 is not available with this version.