Comprehensive lncRNA and mRNA profiles in peripheral blood mononuclear cells derived from ankylosing spondylitis patients by RNA-sequencing analysis.

ABSTRACT: The present study aimed to investigate the comprehensive expression profiles of long non-coding RNA (lncRNA) in ankylosing spondylitis (AS).The peripheral blood samples were collected from 6 AS patients and 6 age- and gender-matched healthy controls (HCs), and separated for peripheral blood mononuclear cells, followed by RNA-sequencing. Further bioinformatics analyses were performed to explore the significantly enriched biological processes, signaling pathways of differentially expressed lncRNAs (DElncRNAs) (based on cis-target and trans-target genes) and differentially expressed mRNAs (DEmRNAs).Principal component analysis plots indicated that both lncRNA and mRNA expression profiles could distinguish AS patients from HCs; heatmap diagram exhibited a relatively good consistency and tendency of lncRNA and mRNA expression profiles in AS patients and HCs, respectively; volcano plots exhibited 114 upregulated and 45 downregulated DElncRNAs, 284 upregulated and 435 downregulated DEmRNAs in AS patients compared with HCs; Gene ontology enrichment analyses indicated that DElncRNAs (based on cis-target and trans-target genes) and DEmRNAs were enriched in molecular functions (including DNA binding, protein binding, etc) and biological process (including immune response, inflammatory response, etc); Kyoto Encyclopedia of Genes and Genomes enrichment analyses revealed that these DElncRNAs (based on cis-target and trans-target genes) and DEmRNAs were enriched in immune and inflammation-related signaling, such as B cell receptor signaling pathway, TNF signaling pathway, NF-kappa B signaling pathway, etc.Our study displays the comprehensive expression profiles and functions of lncRNAs involved in AS, which provides reference for further researches discovering candidate lncRNAs with value in assisting early AS diagnosis.

Introduction

Ankylosing spondylitis (AS), as 1 type of axial spondyloarthritis, frequently occurs before age 40 with male predominance, affecting approximately 0.5% of the population.[ AS is characterized by a chronic inflammatory disease affecting the axial skeleton, meanwhile, its clinical articular-related manifestations include inflammatory back pain, limited spinal mobility, radiographic sacroiliitis, peripheral arthritis due to inflammation erosive osteopenia, and excess spinal bone formation.[ In addition, AS is also accompanied by some extra-articular clinical manifestations such as uniocular anterior uveitis, gut inflammation, osteopenia, apical pulmonary fibrosis, which often develop in active AS patients.[ Due to the gradually progressive nature and the common presentation of the onset, the recognition of early AS is often delayed for about 6 to 10 years, besides, current AS management mainly includes pharmacological treatment, such as non-steroidal anti-inflammatory drugs and tumor necrosis factor inhibitors, which effectively alleviates clinical symptoms, improves functioning, and elevates quality of life in AS patients, however, still a portion of AS patients do not respond to these medications.[ Given the diagnostic delay and insufficient medication options, it is essential to discover the underlying AS mechanism, which may provide more information concerning the disease pathogenesis, and further help in development of treatment and monitoring prognosis for AS management.

Long non-coding RNAs (lncRNAs) constitute a class of transcripts with more than 200 nucleotides which lack protein-coding capability, and existing numerous researches have demonstrated that lncRNAs disclose a wide spectrum of biological functions through epigenetic, transcriptional, posttranscriptional mechanisms and also are involved in physiology processes, such as regulating transcription, modulating mRNA processing, and interaction with proteins.[ Recently, the emerging role of lncRNAs in regulating pathological process of autoimmune diseases, including AS, has been uncovered.[ For example, 1 study shows that lncRNA USP50-2, lncRNA ZNF354A-1, and lncRNA LNI54-1 are involved in the interaction between bone morphogenetic protein 2 and Noggin, besides, they are also participated in mediating osteogenic differentiation of mesenchymal stem cells in AS.[ Furthermore, lncRNA MEG3, lncRNA TUG1, lncRNA AK001085, and lncRNA H19 are revealed to be aberrantly expressed in AS patients compared with healthy controls (HCs), and are correlated with AS susceptibility.[ However, the functions of the vast majority of lncRNAs in AS have not been investigated. In the present study, RNA-sequencing and bioinformatics analyses were performed in order to identify the comprehensive differentially expressed profiles of lncRNA and mRNA involved in AS, meanwhile, recognition of differentially expressed lncRNA (DElncRNA)-targeted genes, construction of the lncRNA-mRNA regulatory network was further conducted, which might help with discovery of candidate lncRNAs with potential for assisting early AS diagnosis.

Methods

Sample collection and processing

A total of 12 peripheral blood samples were collected from 6 AS patients and 6 HCs. All AS patients were older than 18 years and had a diagnosis of AS in accordance with EULAR criteria.[ The peripheral blood samples of AS patients were collected before treatment. HCs were age- and gender-matched healthy subjects whose healthy status was confirmed by medical examination. The current study was approved by Institutional Review Board, and all subjects provided the written informed consents. Peripheral blood samples were extracted using anticoagulant tube, then the density gradient centrifugation was carried out to separate the peripheral blood mononuclear cell, which was stored at −80°C for further use.

Determination of clinical characteristics

After enrollment, the clinical characteristics of HCs and AS patients were collected, including the age, gender, inflammatory status. Besides, for AS patients particular, their human leukocyte antigen-B27 (HLA-B27), disease activity and their signs and symptoms were also collected.

RNA-sequencing

Total RNA was separated from the peripheral blood mononuclear cell using TRIzol Reagent (Thermo Fisher Scientific, Waltham, MA) in strict accordance with the manufacturer's instructions. After isolation, the concentration and purity of total RNA were identified using NanoDropND-2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA), and RNA integrity was checked using Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA). Further purification was carried out using RNA Clean XP Kit (Beckman Coulter, Inc, Kraemer Boulevard Brea, CA) and the RNase-Free DNase Set (QIAGEN, GmbH, Germany). The ribosomal RNA was removed using Ribo-Zero Magnetic Gold Kit (Human) (Illumina, San Diego, CA). Library construction and RNA-sequencing procedures were conducted according to the methods described in previously studies.[ RNA-sequencing was carried out on Illumina HiSeq 2500 (Illumina Inc, San Diego, CA).

RNA-sequencing data analysis

Raw data processing was performed as previous study described.[ In brief, quality control for raw sequencing reads were determined by FastQC, then clean data were obtained, which were aligned to the human genome (GRCh38) using the TopHat 2.0 program,[ and the gene counts were calculated by featureCounts. Differential expression analysis was performed using DESeq2 package in R software (http://cran.r-project.org/bin/windows/base/). DElncRNA and differentially expressed mRNAs (DEmRNA) between AS and HCs were identified based on fold change > 2 and adjusted P value < .05 determined by Benjamin-Hochberg corrected multiple t test.

Bioinformatics analysis

LncRNA and mRNA expression profiles were analyzed by principal component analysis (PCA) plot and heatmap plot using Factoextra and Pheatmap packages in R software, respectively. DElncRNA and DEmRNA were illustrated by volcano plots using ggplot2 package in R software. Target genes of DElncRNA were predicted according to the trans-regulation and cis-regulation then determined by calculating the Pearson correlation coefficients and P values. The cellular component (CC), molecular function (MF), and biological process (BP) of DElncRNA (based on cis-target gene and trans-target gene) and DEmRNA were revealed by Gene ontology (GO) enrichment analysis (http://www.geneontology.org/). The signal pathways associated with DElncRNA (based on cis-target gene and trans-target gene) and DEmRNA were interpreted by Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis (http://www.genome.jp/kegg/). The lncRNA-mRNA regulatory network was built based on cis-target genes and trans-target genes of DElncRNA using igraph package in the R software. The cis-target genes were identified from the genes with a distance <20 kilobases from DElncRNA, and the trans-target genes were predicted by the LncTar software (https://www.cuilab.cn/lnctar). There were too many DElncRNAs and corresponding trans-target genes to clearly illustrate, as a result, we built the lncRNA-mRNA regulatory network based on trans-target genes of the top 30 DElncRNAs to display. The top 30 DElncRNAs including the top 15 downregulated DElncRNAs and the top 15 upregulated DElncRNAs were selected by ranking of absolute value of Log2(fold change). Ultimately, the findings by bioinformatics analysis were aggregated into an OmicCircos plot.

Results

Characteristics of HCs and AS patients

The age was 31.8 ± 8.5 years and 29.0 ± 5.6 years in HCs and AS patients, respectively (P = .510, Table S1, Supplemental Digital Content). In addition, 5 (83.3%) subjects were male, while 1 (16.7%) subject was female in HCs; in terms of the AS patients, 5 (83.3%) patients were male, while 1 (16.7%) patient was female (P = 1.000). The C-reactive protein (4.3 ± 2.4 mg/L vs 36.7 ± 16.8 mg/L, P = .005) and erythrocyte sedimentation rate (10.4 ± 6.6 mm/H vs 38.6 ± 17.0 mm/H, P = .004) were lower in HCs compared with AS patients. Furthermore, other detailed characteristics of AS patients were shown in Table S1, Supplemental Digital Content.

PCA plot and heatmap analyses

PCA plot analyses indicated that there was a clear segregation of lncRNA (Fig. 1A) and mRNA expression profiles (Fig. 1B) between AS patients and HCs. Furthermore, heatmap analyses displayed a relatively good consistency and tendency of lncRNA (Fig. 1C) and mRNA (Fig. 1D) expression profiles in AS patients and HCs, respectively.

Figure 1

LncRNA and mRNA expression files. PCA plot (A) and heatmap (B) analyses for lncRNA expression files in AS patients and HCs. PCA plot (C) and heatmap (D) analyses for mRNA expression files in AS patients and HCs. AS = ankylosing spondylitis, HCs = healthy controls, lncRNA = long non-coding RNA, mRNA = messenger RNA, PCA = principal component analysis.

Volcano plot analyses

The volcano plot analyses exhibited that there were 159 DElncRNAs (including 114 upregulated and 45 downregulated DElncRNAs) in AS patients compared with HCs (Fig. 2A). Furthermore, the volcano plot analyses illustrated that there were 719 DEmRNAs (including 284 upregulated and 435 downregulated DEmRNAs) in AS patients compared with HCs (Fig. 2B).

Figure 2

Upregulated and downregulated DElncRNAs and DEmRNAs. Volcano plot analyses for lncRNA (A) and mRNA (B) expression profiles in AS patients compared with HCs. AS = ankylosing spondylitis, DElncRNAs = differentially expressed long non-coding RNAs, DEmRNAs = differentially expressed messenger RNAs, HCs = healthy controls, lncRNA = long non-coding RNA.

GO enrichment analyses

GO enrichment analysis based on trans-target gene indicated that DElncRNAs were enriched in MFs (including RNA polymerase II core promoter proxima, DNA binding, protein binding, etc), CCs (including nucleus, external side of plasma membrane, actin cytoskeleton, etc), and BPs (including regulation of cell death, immune response, negative regulation of transcription, etc) (Fig. 3A); GO enrichment analysis based on cis-target gene revealed that DElncRNAs were enriched in MFs (including protein tyrosine/threonine phosphatase, protein binding, identical protein binding, etc), CCs (including plasma membrane, nucleus, cytoplasm, etc), and BPs (including apoptotic process, signal transduction, endoderm formation, etc) (Fig. 3B); GO enrichment analysis exhibited that DEmRNAs were enriched in MFs (including DNA binding, transcription factor activity, sequence, RNA polymerase II core promoter proxima, etc), CCs (including nucleus, integral component of plasma membrane, external side pf plasma membrane, etc), and BPs (including inflammatory response, chemokine-mediated signaling pathway, regulation of cell death, etc) (Fig. 3C).

Figure 3

GO enrichment analyses of AS-associated DElncRNAs and DEmRNAs. MFs, CCs, and BPs that DElncRNAs were enriched in according to the GO enrichment analyses based on trans-target genes (A) and cis-target genes (B). MFs, CCs, and BPs that DEmRNAs were enriched in according to the GO enrichment analyses (C). AS = ankylosing spondylitis, BPs = biological processes, CCs = cellular components, DElncRNAs = differentially expressed long non-coding RNAs, DEmRNAs = differentially expressed messenger RNAs, GO = Gene ontology, MFs = molecular functions.

KEGG enrichment analyses

KEGG enrichment analysis based on trans-target gene indicated that DElncRNAs were enriched in B cell receptor signaling pathway, TNF signaling pathway, NF-kappa B signaling pathway, etc (Fig. 4A); KEGG enrichment analysis based on cis-target gene revealed that DElncRNAs were enriched in mitogen-activated protein kinase signaling pathway, B cell receptor signaling pathway, pathways in cancer, etc (Fig. 4B); KEGG enrichment analysis exhibited that DEmRNAs were enriched in cytokine-cytokine receptor interaction, B cell receptor signaling pathway, TNF signaling pathway, etc (Fig. 4C).

Figure 4

KEGG enrichment analyses of AS-associated DElncRNAs and DEmRNAs. Signaling pathways that DElncRNAs were enriched in according to the KEGG enrichment analyses based on trans-target genes (A) and cis-target genes (B). Signaling pathways that DEmRNAs were enriched in according to the KEGG enrichment analyses (C). AS = ankylosing spondylitis, DElncRNAs = differentially expressed long non-coding RNAs, DEmRNAs = differentially expressed messenger RNAs, KEGG = Kyoko Encyclopedia of Genes and Genomes.

Regulatory network of lncRNA-mRNA

As there were too many trans-target genes of DElncRNAs, only trans-target genes of the top 30 DElncRNAs were presented in regulatory network (Fig. 5A). In brief, there were 17 DElncRNAs with multiple (≥3) trans-target genes, 7 DElncRNAs with 2 trans-target genes, 6 DElncRNAs with single trans-target gene. However, since the number of cis-target genes of DElncRNAs was limited, all cis-target genes of DElncRNAs were displayed in regulatory network (Fig. 5B). A total of 26 DElncRNAs were presented, among which 25 DElncRNAs were with single cis-target gene, and only 1 DElncRNA was with 2 cis-target genes. The top 30 DElncRNAs (including the top 15 upregulated DElncRNAs and the top 15 downregulated DElncRNAs) were shown in Table 1.

Figure 5

Merged regulatory network of DElncRNAs and their targets. The regulatory network of the top 30 DElncRNAs with their trans-target genes (A); the regulatory network of all DElncRNAs with their cis-target genes (B). DElncRNAs = differentially expressed long non-coding RNAs, DEmRNAs = differentially expressed messenger RNAs.

Table 1

Top 30 DElncRNAs (including top 15 down-regulated DElncRNAs and top 15 up-regulated DElncRNAs)∗.

Gene ID	Chr	Symbol	Biotype	Log2(FC)	P value (AS vs HC)	Padj value (AS vs HC)	Trend
ENSG00000267896	chr19	AC018766.4	Antisense	−4.271	1.5E-07	5.3E-06	DOWN
ENSG00000237090	chr1	RP11-342M1.6	Processed_pseudogene	−4.199	1.7E-03	1.1E-02	DOWN
ENSG00000245869	chr11	RP11-158I9.5	Antisense	−4.145	5.3E-05	6.8E-04	DOWN
ENSG00000272168	chr6	CASC15	LincRNA	−3.606	2.7E-03	1.6E-02	DOWN
ENSG00000247157	chr12	LINC01252	LincRNA	−3.563	3.2E-03	1.8E-02	DOWN
ENSG00000246695	chr12	RASSF8-AS1	Antisense	−3.517	1.4E-03	9.6E-03	DOWN
ENSG00000255418	chr11	RP11-266A24.1	LincRNA	−3.166	5.5E-03	2.8E-02	DOWN
ENSG00000272273	chr6	XXbac-BPG252P9.10	Antisense	−3.140	7.2E-03	3.4E-02	DOWN
ENSG00000251348	chr5	HSPD1P11	Processed_pseudogene	−3.123	9.9E-03	4.4E-02	DOWN
ENSG00000255026	chr11	RP11-326C3.2	Antisense	−2.954	3.9E-03	2.1E-02	DOWN
ENSG00000233723	chr2	LINC01122	LincRNA	−2.949	9.5E-03	4.2E-02	DOWN
ENSG00000260418	chr6	RP3-406A7.7	LincRNA	−2.646	1.1E-02	4.8E-02	DOWN
ENSG00000205056	chr12	RP11-693J15.5	LincRNA	−2.296	7.7E-05	9.1E-04	DOWN
ENSG00000257275	chr14	RP11-164H13.1	Antisense	−2.263	6.9E-03	3.3E-02	DOWN
ENSG00000258897	chr14	EGLN3-AS1	Antisense	−2.095	4.5E-05	5.9E-04	DOWN
ENSG00000267288	chr17	RP13-890H12.2	Antisense	4.794	1.60E-03	1.10E-02	UP
ENSG00000260765	chr16	CES1P2	Unprocessed_pseudogene	3.902	9.50E-03	4.20E-02	UP
ENSG00000213058	chr1	RP4-765C7.2	Processed_pseudogene	3.899	5.40E-03	2.70E-02	UP
ENSG00000255491	chr8	RP11-1082L8.4	LincRNA	3.466	2.00E-03	1.20E-02	UP
ENSG00000239569	chr7	KMT2E-AS1	Antisense	3.447	9.60E-05	1.10E-03	UP
ENSG00000214407	chr3	RP11-221J22.1	LincRNA	3.343	1.60E-03	1.10E-02	UP
ENSG00000260823	chr16	RP11-249C24.10	LincRNA	3.29	9.00E-03	4.10E-02	UP
ENSG00000229140	chr8	CCDC26	LincRNA	3.161	5.00E-07	1.40E-05	UP
ENSG00000273272	chr22	CTA-384D8.34	LincRNA	3.099	3.90E-03	2.10E-02	UP
ENSG00000254639	chr11	CTD-2589M5.5	LincRNA	3.067	2.50E-03	1.50E-02	UP
ENSG00000254704	chr11	RP11-1036E20.7	Processed_pseudogene	3.046	6.80E-04	5.30E-03	UP
ENSG00000259884	chr12	RP11-1100L3.8	LincRNA	3.022	1.70E-03	1.10E-02	UP
ENSG00000258820	chr14	RP11-293M10.2	Antisense	2.86	2.20E-03	1.40E-02	UP
ENSG00000255363	chr11	RP11-672A2.5	LincRNA	2.835	7.60E-03	3.60E-02	UP
ENSG00000260805	chr1	RP11-61J19.5	LincRNA	2.806	5.60E-03	2.80E-02	UP

AS = ankylosing spondylitis, Chr = chromosome, DElncRNAs = differentially expressed lncRNAs, FC = fold change, HC = healthy control, Padj = adjusted P value.

Top 30 DElncRNAs were selected by ranking of absolute value of Log2(FC).

OmicCircos plot for bioinformatics analysis

For the purpose to comprehensively describe bioinformatics analyses of DElncRNAs and DEmRNAs, OmicCircos plot was presented (Fig. 6). The outermost layer represented chromosome number; the second outermost layer represented upregulated and downregulated mRNAs, which were in red and green color, respectively; the third outermost layer represented upregulated and downregulated lncRNAs, which were in red and green color, respectively; the collected lines in the center of OmicCircos plot revealed trans- or cis-regulation among lncRNAs and mRNAs. This OmicCircos plot systematically illustrated the locations, expressions and regulatory network of DElncRNAs and DEmRNAs.

Figure 6

Comprehensive bioinformatics analyses of DElncRNAs and DEmRNAs. DElncRNAs = differentially expressed long non-coding RNAs, DEmRNAs = differentially expressed messenger RNAs.

Discussion

LncRNAs are considered as an abundant class of RNAs without protein-coding ability, and are divided into 5 major biotypes including antisense lncRNA, intergenic lncRNA, pseudogenic lncRNA, bidirectional lncRNA, enhancer-associated lncRNA, which all possess multiple MFs, such as involvement of mRNA processing, transcriptional regulation, etc.[ Moreover, recent emerging publications display that the dysregulation of lncRNAs play an important role in the development and progression of multiple innate immune-mediated inflammatory diseases, including AS.[ For example, 1 study reveals that lncRNA NKILA is upregulated in patients with active AS compared to HCs, and its overexpression is correlated with active disease, and elevated duration of hospitalization in AS patients.[ However, only a small portion of lncRNAs implicated in the AS pathogenesis has been identified,[ and large majority of lncRNAs are not investigated. Hence, we conducted the present study to extensively investigate lncRNA expression profiles in AS, which might provide reference for further researches discovering candidate lncRNAs with value in assisting early AS diagnosis.

In our study, initially, PCA plot and heatmap plot analyses illustrated that lncRNA and mRNA expression profiles could distinguish AS patients from HCs; and further volcano plot displayed that there were 159 DElncRNAs and 719 DEmRNAs in AS patients compared with HCs, which suggested that these differentially expressed genes might be involved in the development of AS, and have potential to be biomarkers of AS. Following that, to achieve a deeper understanding of these differentially expressed genes, we conducted GO enrichment analysis based on target genes of DElncRNAs (including trans-regulation and cis-regulation) and DEmRNAs, which exhibited that DElncRNAs were enriched in BPs related to immune and inflammatory responses. In terms of these findings, we provided the interpretations as follows DElncRNAs might regulate the expressions of their adjacent coding genes (cis-regulating) due to their co-expression, mediating BPs related to immune system, including regulating autophagy, proliferation and apoptosis of T cell, MAPK signaling pathway in AS. These data agreed with the previous evidence that dysregulated lncRNAs could regulate immune response via selecting nearby candidate mRNAs in other inflammatory-related disease.[ DElncRNAs might induce inflammatory-related targeted mRNAs at transcriptional level on different chromosomes (trans-regulation) in response to lipopolysaccharide signaling, being involved in the inflammasome secretion, macrophage activation of AS.[

In our study, further KEGG analyses disclosed that DElncRNAs were enriched in immune and inflammation-related signaling pathways, such as NF-kappa B signaling pathway, MAPK signaling pathway, B cell receptor signaling pathway, TNF signaling pathway, etc. Regarding these signaling pathways, we provided our understanding as follows: DElncRNAs might bind with the proteins in NF-kappa B signaling pathway which controlled a diverse range of genes encoding for inflammatory cytokines and chemokines, promoting the pro-inflammatory actions, thereby increasing AS risk and leading to active status of AS.[ Based on previous evidence, MAPK signaling pathway was implicated in the development of Th1 cells and the followed production of IFN-γ and macrophage activation.[ Therefore, it was speculated that DElncRNAs might correlate with AS risk via regulating its target gene-mediated MAPK signaling. For example, lncRNA RP11-164H13.1 might mediate its corresponding trans-target gene (TCL1A), activating MAPK signaling, and further enhancing the secretion of IL-10 as well as activation of T cell, mediating immunopathogenesis process of AS.[ According to existing papers, B cell receptor signaling was important for development of normal B cell and adaptive immunity, and abnormal activated B cells was correlated with AS risk.[ Given these evidences, it was speculated that DElncRNAs might trigger B cell receptor signaling activation and further lead to B-cell activation via binding of ligand to the B cell receptor, which contributes to the cascade of systematic inflammation as well as enhanced AS risk.[ Moreover, based on previous evidence, TNF signaling pathway regulated the development and survival of osteoclasts and triggered inflammatory responses in pathology of inflammatory diseases.[ DElncRNAs were speculated to influence the transcriptional level of proteins on TNF signaling pathway, mediating a number of pro-inflammatory changes in AS.[ Furthermore, the regulating action between DElncRNAs with their target genes (such as TLR10, TCL1A) was shown in our study, and these target genes were reported to play important roles in T-cell maturing progression, activation of toll like receptor signaling, which was reported to be implicated in the etiology of AS by previous publications.[ Taken together, these DElncRNAs play critical roles in initiation of AS via regulating potential inflammatory-related signaling pathways, which might have potential to be biomarkers for disease susceptibility and therapeutic targets in AS management.

Our study was a preliminary research, which still existed some limitations as follows: the sample size in our study was relatively small, which needed a larger sample size for validation. The present study was a preliminary study indicating comprehensive lncRNA expression profiles in AS, however detailed regulatory mechanism of these DElncRNAs in AS needed further studies. DElncRNAs found in our study needed to be validated by more detecting methods (such as real-time quantitative polymerase chain reaction) in larger AS populations. The HLA B27 status of controls was not clear, which might be a potential confounding factor in the present study and should be eliminated in the further study.

Conclusion

In conclusion, our study displays the comprehensive expression profiles and functions of lncRNAs involved in AS, which provides a direction for further study on lncRNAs in AS management.

Author contributions

Conceptualization: Chuangxin Li, Wa Qu.

Data curation: Chuangxin Li, Wa Qu, Xuefeng Yang.

Formal analysis: Chuangxin Li, Wa Qu.

Investigation: Chuangxin Li, Wa Qu, Xuefeng Yang.

Methodology: Chuangxin Li, Wa Qu, Xuefeng Yang.

Project administration: Chuangxin Li, Wa Qu.

Resources: Chuangxin Li, Wa Qu, Xuefeng Yang.

Software: Chuangxin Li, Wa Qu.

Supervision: Chuangxin Li, Wa Qu.

Validation: Chuangxin Li, Wa Qu, Xuefeng Yang.

Visualization: Chuangxin Li, Wa Qu, Xuefeng Yang.

Writing – original draft: Chuangxin Li, Wa Qu, Xuefeng Yang.

Writing – review & editing: Chuangxin Li, Wa Qu.