Patients with dermatomyositis shared partially similar transcriptome signature with COVID-19 infection

Abstract

Dermatomyositis (DM) is an autoimmune disease that primarily affects the skin and skeletal muscle. Virus infection and type I interferon-related signaling pathways play an important role in the pathogenesis of dermatomyositis. In this study, we found that the skin of patients with DM and the skin of patients with COVID-19 have similar transcriptional profiles, and identified key genes involved in dermatomyositis based on bioinformatics analysis. These hub-genes might be served as potential biomarkers for the early diagnosis and therapy of DM, including MX1, ISG15, IFIT3, IFIT1, RSAD2, IFIT2, IFI6, XAF1, IRF9, MX2.

Keywords:

• Dermatomyositis
• COVID-19
• skin
• bioinformatics analysis

Introduction

Dermatomyositis (DM) is a heterogeneous diffuse connective tissue disease, characterized by skin rash, chronic muscle weakness, and skeletal muscle mononuclear cell infiltration [1]. DM is the most common subtype of idiopathic inflammatory myopathy (IIM). The incidence of IIM is low, and the annual incidence rate is about 1.16 ∼ 19/million people [2]. The onset age has two peak periods, 10–15 years old and 45–60 years old, and the average age of onset is 40–60 years old [3]. The most important clinical feature of DM is specific skin rash damage, such as Heliotrope rash, V-shaped shawl sign, and Gottron sign. Perivascular lymphocytic infiltration of activated T lymphocytes and macrophages is present in the skin of DM [4].

The etiology of DM is still unclear. Both immune and non-immune mechanisms are involved in its pathogenesis [5–8]. A dysregulated type I interferon response plays an important role in the development of DM [9–11]. Type I IFN excess may be associated with disease severity and autoantibodies [9]. In addition, some studies have suggested that the pathogenesis of dermatomyositis may be related to viral infection, especially EB virus infection [12, 13]. The current epidemic of coronavirus disease 2019 (COVID-19) is caused by infection with SARS-CoV-2. Recent study indicates that COVID19 infection and DM may share some of the same pathogenesis, and both are associated with excessive activation of the type I interferon pathway [14, 15]. And there are also studies suggesting that COVID-19-related myositis may be dermatomyositis [16]. Therefore, it is necessary to re-analyze the differential genes between dermatomyositis and healthy control. This will provide useful information for the early diagnosis and treatment of dermatomyositis and is expected to provide new targets.

Specific skin lesions are the most important clinical manifestations of dermatomyositis and one of the diagnostic criteria. Therefore, the analysis of differential genes in dermatomyositis skin lesions is extremely important. In this study, we downloaded the GSE46239 mRNA expression profiles of DM’s skin and GSE193068 of COVID-19 patients’ skin from the website of Gene Expression Omnibus (GEO). Differentially expressed genes (DEGs) analyses were carried out by using the online analysis tool GEO2R. Gene ontology and pathway enrichment analysis were performed by R package clusterProfiler. Next, we constructed the protein-protein interaction (PPI) network of DEGs and identified these genes (module genes) that appeared in two or more analyses. Then, we selected the hub genes by cytoHubba. Finally, quantitative PCR-based detection of hub genes in muscle tissue and immunohistochemistry of skin tissue were performed. In this study, we have discovered some key genes and important signal pathways related to DM, hoping to provide new targets for the early diagnosis and treatment of DM.

Materials and methods

Microarray data analysis

The microarray expression profiling datasets GSE46239 and GSE193068 were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). GSE46239 was based on the GPL570 Affymetrix Human Genome U133 Plus 2.0 Array, containing 52 samples consisting of 48 DM patients and 4 healthy controls. And the skin of patients with dermatomyositis come from the active part of the skin lesion. GSE193068 was based on GPL31171 Gene Expression Counter ® RCC – Human Immunology V2 panel, containing the skin expression of 10 COVID-19 patients and 4 healthy controls (Table 1).

Table 1. Basic information of the two microarrays.

	Subjects	Platform	Organization name	Country	Submission date/updated date	Number of samples (patients/NC)	Reference
GSE46239	The skin of active rashes of DM patients.	GPL570	Medimmune LLC	USA	2013.4.20/2020.12.09	48/4
GSE193068	COVID-19 skin lesions	GPL31171	CHUV	Switzerland	2022.1.05/2022.3.17	10/4	[17]

DEGs analysis

DEGs analysis was performed through using the online analysis tool GEO2R. The expression profiles of patients and healthy controls were compared to identify the DEGs. P-values and adjusted P-values were calculated using t-tests. Genes from each sample with the following criteria were retained: |log2 (fold-change)|≥1 and P value < 0.05. Venn diagram was used to analyze the intersection of DEGs from these two datasets.

Gene ontology and pathway enrichment analysis

The Gene ontology (GO) [18] is a description of gene function according to biological process (BP), molecular function (MF) and cell composition (CC). Kyoto Encyclopedia of Genes and Genomes(KEGG) [19] is a systematic analysis of the metabolic pathways of gene products in cells and the functions of these gene products, which helps to study genes and expression information as a whole network. We used R package clusterProfiler [20] to view the GO and KEGG enrichment of DEGs of GSE46239. And a p value <0.05 was considered statistically significant.

Protein-protein interaction (PPI) and module analysis

The Search Tool for the Retrieval of Interacting Genes (STRING) database [21] (http://string-db.org/cgi/input.pl) is to study the interaction network between proteins, which helps to mine the core regulatory genes. This database contains 24.6 million proteins and >2 billion interactions including 5090 organisms. We uploaded the DEGs of GSE46239 into the STRING database and set up the interaction score ≥ 0.900 (highest confidence) as the significant threshold with hiding disconnected nodes in the network. Then, we used the Cytoscape software [22, 23] to construct the PPI network. The Molecular Complex Detection (MCODE) [24] is a functional module for clustering construction in a huge gene (protein) network. The PPI network was imported into Cytoscape, and MCODE identified the most important two modules in the PPI network. The selection criteria were as follows: degree cutoff = 2, node score cutoff = 0.2, maximum depth = 100, k score = 2. Then, we used R package clusterProfiler to perform GO and KEGG analysis of the genes in the module.

Hub genes identification and analysis

Cytohubba [25] is a plug-in used by Cytoscape software to identify hub nodes. CytoHubba ranks nodes according to their attributes in the network and provides 11 topological analysis methods. The hub genes were calculated based on the Maximal Chique Centrality (MCC) topological analysis methods by using cytoHubba.

Muscle and skin tissues collection

We collected muscle samples from 24 patients with IIM who were admitted to the Department of Rheumatology, Xiangya Hospital, Central South University (Supplementary Table S1). All of them met the diagnostic criteria of 1975 Bohan/Peter Diagnostic Standard for dermatomyositis and polymyositis [26, 27]. Other autoimmune diseases, infections and tumors were excluded. The control group included 12 patients from Department of Neurology, Xiangya Hospital, Central South University. Skin samples from 4 patients with dermatomyositis and 3 controls were collected from the Department of Rheumatology, Xiangya Hospital, Central South University. Written informed consent has been obtained from each subject. And the study methodologies have been approved by the ethical standards of Xiangya Hospital of Central South University (No. 201212074). All experiments conform to the Declaration of Helsinki.

RNA extraction, cDNA synthesis and qRT-PCR detection

Trizol (Invitrogen/Life Technologies, Carlsbad, CA) was used to extract RNA from muscle tissue according to the instructions. It was then reversed into cDNA using the Takara reverse transcription reagent (Takara, Japan, Code No. RR036A). The mRNA level of hub genes was detected via qRT-PCR with TB Green mixture (Takara, Japan, Code No. RR420A) by using the ABI7500 system (Applied Biosystems, Foster City, CA, USA). The relative expression levels of genes in each sample were calculated by using the housekeeping gene GAPDH. The primers used in the study are shown in Supplementary Table S2.

Immunohistochemistry of skin tissue

Skin samples from patients were collected with 4% paraformaldehyde. The expression of RSAD2 was detected by immunohistochemistry.

Statistical methods

All data were analyzed by GraphPad Prism 7.0 software and expressed as mean ± S.D. Student’s unpaired T test was used to determine the statistical difference between the two groups. When the sample data were not normally distributed, Mann-Whitney test (two-tail test) was used for statistical analysis. And a p value <0.05 was considered statistically significant.

Results

Identification of DEGs in GSE46239 and GSE193068

A total of 380 DEGs were identified in GSE46239 of DM patients, and 271 genes were up-regulated and other 109 genes were down-regulated (Figure 1A, Table 2). Besides, a total of 159 up-regulated genes and 12 down-regulated genes were found in the GSE193068 dataset of COVID-19 patients (Figure 1B). And it was found that there were co-expressed differential genes in these two datasets (Figure 1C). The most significant pathways of these common genes are NOD-like receptor signaling way, viral protein interaction with cytokine and cytokine receptor, and coronavirus disease-COVID-19.

Figure 1. Volcano plots of datasets. (A) Volcano plots of GSE46239 in DM. The ordinate represents log2 (fold change), the abscissa represents log10 (p value), each point represents a gene, |log2 FC| ≥1and p < 0.05 as thresholds. (B) Volcano plots of GSE193068 in COVID-19. (C) The intersection of the two data sets. (D) KEGG pathway of these common genes.

Gene term enrichment analysis recognized similar a signature of COVID-19 infection in patients with DM

To further explore potential functions of these DEGs in DM, we performed GO and pathway analyses on DM by using the criterion of p value <0.05. As shown in Figure 2, it shows the top 10 significant terms for the BP, CC, MF and KEGG of DEGs in GSE46239. The DEGs were enriched for genes in the BP term involved in response to virus, response to interferon-gamma and type I interferon signaling pathway. In the MF term, the DEGs were enriched in receptor ligand activity, signaling receptor activator activity and glycosaminoglycan binding. In the CC term, the DEGs were enriched in collagen-containing extracellular matrix, external side of plasma membrane, and early endosome. The KEGG enrichment analysis showed that the most significant pathways were Coronavirus disease-COVID-19, influenza A and NOD-like receptor signaling pathway.

Figure 2. Gene ontology analysis and significantly enriched signal pathway of differentially expressed genes (DEGs) in DM. A. BP analysis of DEGs. B. MF analysis of DEGs. C.CC analysis of DEGs. D.KEGG analysis of DEGs.

PPI networks and module analysis of 380 DEGs in GSE46239

All of the 380 DEGs (271 up-regulated genes and 109 down-regulated genes) were uploaded onto the online website STRING (http://string-db.org) and then they were analyzed by using the Cytoscape software. The PPI network contained 173 nodes and 829 edges, and 207 genes did not fall into the PPI network (Figure 3A). Then we selected the two significant modules by using the plug-in MCODE. Module 1 consisted of 25 nodes/genes and 296 edges (Figure 3B), which are mainly associated with cellular response to type I interferon signaling pathway (BP), double-stranded RNA binding (MF), postsynaptic endocytic zone (CC), and hepatitis C(KEGG) (Figure 4, module1). Module 2 consisted of 23 nodes/genes and 117 edges (Figure 3C), which are mainly associated with cellular response to interferon-gamma (BP), immune receptor activity (MF), external side of plasma membrane (CC), and viral protein interaction with cytokine and cytokine receptor (KEGG) (Figure 4, module 2).

Figure 3. Protein–protein interaction (PPI) network of differentially expressed genes (DEGs) in DM. A. Based on the STRING online database, 173 genes were filtered into the DEG PPI network. B. The most significant module 1 (25 genes) from the PPI network. C. The most significant module 2 (23 genes) from the PPI network.

Figure 4. Functional and Pathway Enrichment of Module Genes. A. BP analysis of Module Genes. B. MF analysis of Module Genes. C.CC analysis of Module Genes. D.KEGG analysis of Module Genes.

The identified hub genes in patients with DM participates in COVID-19 infection

In this part, we used the cytoscape plug-in cytoHubba to confirm the hub gene according to the latest analysis method MCC analysis method which is strongly recommended. The top 10 genes with the highest MCC scores are considered hub genes (Figure 5A, Table 2). And the pathways enrichment of hub genes were associated with Coronavirus disease-COVID-19, influenza A and hepatitis C (Figure 5B, Table 3).

Figure 5. Hub gene and KEGG pathway analysis in DM. A. Top 10 hub genes with highest MCC score in cytoHubba. B. KEGG pathway analysis of hub genes.

Table 2. 380 DEGs were identified in GSE46239, including 271 up-regulated genes and 109 down-regulated genes.

DEGs	Gene name
Up-regulated	STEAP4, FAM46A, IFI6, IFITM1, LOC101928916, ERG, ADGRL4, CNN3, IFITM2, ANGPT2, SEMA3C, COL6A1 , IFI44L, ADAMTS1, GIMAP6, IFIT3, FAM162B, ISG15, CMPK2, GIMAP7, SIGLEC1, RSAD2, MX1, IFI44, EPSTI1, IFIT1, CYBB, SERPING1, BST2, FAM198B, IFITM3, NRP2, HERC5, GIMAP4, C1QC, PLSCR1, IFI27, GIMAP8, FLI1, NT5E, TFPI, C1QA, C1R, ADAMTS5, SEC24D, THBS1, PRRX1, TM4SF18, MX2, LY96, LY6E, DDX60L, IFIT5, IFITM2, ADM, LAP3, CD163, LOC101930489, SAMHD1, TNFAIP6, CFI, OAS2, CTHRC1, LGALS3BP, STAT1, MDK, SRGN, NID1, SLC15A3, EIF2AK2, SP110, DDX58, FGF7, C1S, OASL, HERC6, EGFL6, VCAN, SP140L, MAF, IGF1, TLR4, STARD4, IFIT2, OAS3, TRIM5, IRF9, IFI16, CP, XAF1, C1QB, SP100, TNFSF13B, APOL6, THEMIS2, RNASE1, MCTP1, GREM1, TRIM22, ACKR4, LGALS9, LUM, BLZF1, STAT2, RSPO3, SAMD9, USP18, MS4A7, NDC80, CCL8, LAIR1, IKBIP, MS4A4A, C3AR1, ICAM2, CCR1, NR2F2, ANGPTL1, APOL3, APOL2,PARP9, SOCS1, EMILIN2, TNFSF10, ISG20, CARD8-AS1, GBP3, DDX60, CLIC2, TAP2, MS4A6A, VGLL3, CYR61, PARP14, CFB, MICB, IFIH1, APOBEC3G, ADAMDEC1, TFEC, HELZ2, C19orf66, GIMAP2, SAMD9L, NCEH1, RARRES3, SLFN12, THY1, PTN, ST8SIA4, CASP1, DTX3L, NEXN, CXCL10, MMRN1, CLEC7A, GBP1, IFI35, CXCR4, IGFBP6, FCGR1, AVPR1A, SOCS3, CTGF, SNX10, AIM2, RNF213, PNPT1, OAS1, CXCL2, IRF7, UBE2L6, CTSS, PLA1A, SLFN5, SPOCK1, FCGR1B, C1orf162, CXCL11, GBP5, BATF2, RTP4, ICAM1, ITGA4, KYNU, IRF8, APOBEC3, SLC44A5, DKK1, CCL5, SLC2A3, GNLY, TMEM100, GSTT2, PCOLCE2, SULF1, PRH1-PRR4, HAVCR2, PARP12, FAM26F, ODF3B, NMI, DPP4, PLAC8, CARD16, VSIG4, GALNT15, ZBP1, KLHDC7B, TYMP, TAGAP, LYN, PSMB9, LILRB2, VCAM1, TLR3, CD36, CCL2, MMP9, CH25H, WARS, SERPINA1, TDO2, APOBEC3A, GBP4, FCER1G, MIR21, APOL1, LCP2, GPR183, KLRC, ZC3HAV1, SERPINE1, BTN3A3, 1-Mar, IFI3, IDO1, NKG7, IGSF6, MT1M, LIPA, IL18BP, MFAP5, XCL, TIMD4, GZMB, ERAP2, MPZL2, PLA2G2A, CXCL9, XCL1, XIST, SELE, FPR3, SAMSN1, RGS1, CCL3,CD86, ADORA3, CCR7, BCL2A1, IGFL1, LAMP3, ZNF595, CD274, CD69, SELL, SAA2, SLAMF7, SPP1, PTPRC
Down-regulated	WNT4, BRWD3, KMT2A, TLCD2, LYST, HECTD1, SPPL3, SDHAP2, CHD4, KANSL1, FTX, LOC158863, GNAL, GPR158, COL6A5, SLC1A2, FRMD4A, LRRC37A3, FGFR2, ADGRG6, LINC01550, LOC101930303, KIT, MYB, DPP6, CLN8, MIR205HG, MYH10, TAPT1-AS1, NFYA, MAB21L3, BAGE2, SF1, HERC1, MIR181A2HG, CD99P1, LOC101928524, PDXDC1, ABHD5, POU2F3, CD1E, WNT2B, TRPM1, PLP1, C11orf80, DDX6, XG, LOC729732, NHLH2, ABCA5, ITGB8, ZDHHC11, ALDH3A1, LOC389834, DCT, MAGI1, CXCL14, FOXQ1, GUSBP1, PCSK2, IPO5P1, ALDH3A2, SOX6, FCGBP, ARHGEF28, PAIP2B, CDHR1, TPPP, TFAP2B, SMAD5-AS1, PDE6A, SLC47A1, BPY2, MTURN, CD1A, CNTN1, APOE, ZNF254, SLC14A1, CD1C, CD207, LOC101930123///EEF2K, CHP2, MLANA, PON3, MUC1, GPLD1, GABPB1-AS1, KDM5D, AGPAT3, AATBC, OLAH, RPS4Y1, USP9Y, ACSBG1, TMEM97, TYR, TTTY14, EIF1AY, CRAT, SLC27A6, DDX3Y, FADS1, TXLNGY, HLA-DQA1, TF, GAL, FCER1A,CNTNAP3

Table 3. Functional roles of top 10 hub genes with highest MCC score.

No.	Gene symbol	Full name	MCC score	Function
1	MX1	MX dynamin like GTPase 1	2.25E + 21	Interferon-induced dynamin-like GTPase participating in the cellular antiviral response.
2	ISG15	ISG15 ubiquitin-like modifier	2.25E + 21	Chemotactic activity towards neutrophils, direction of ligated target proteins to intermediate filaments, cell-to-cell signaling, and antiviral activity during viral infections.
3	IFIT3	interferon induced protein with tetratricopeptide repeats 3	2.25E + 21	Inhibiting cellular as well as viral processes, cell migration, proliferation, signaling, and viral. Negatively regulates the apoptotic effects of IFIT2.
4	IFIT1	interferon induced protein with tetratricopeptide repeats 1	2.25E + 21	Inhibiting the expression of viral messenger RNAs.
5	RSAD2	radical S-adenosyl methionine domain containing 2	2.25E + 21	Inhibiting viral replication and translational initiation.
6	IFIT2	interferon induced protein with tetratricopeptide repeats 2	2.25E + 21	Inhibiting expression of viral messenger RNAs and promoting apoptosis.
7	IFI6	interferon alpha inducible protein 6	2.25E + 21	Induced by interferon, playing a critical role in the regulation of apoptosis.
8	XAF1	XIAP associated factor 1	2.25E + 21	A negative regulator of members of the IAP that inhibits anti-caspase activity of BIRC4.
9	IRF9	interferon regulatory factor 9	2.25E + 21	Virus-mediated activation of interferon, and modulation of cell growth, differentiation, apoptosis, and immune system activity.
10	MX2	MX dynamin like GTPase 2	2.25E + 21	Interferon-induced dynamin-like GTPase with potent antiviral activity against HIV-1.

Hub gene expression in the muscle tissue of patients with IIM

24 muscle tissue samples of IIM were obtained. 10 hub genes were selected for analysis by qPCR in IIM muscle tissues. The expression of these hub genes in IIM muscle tissues was significantly higher than that in the NC group (p < 0.05, Figure 6).

Figure 6. Hub gene expression in the muscle tissue. (IIM, idiopathic inflammatory myopathy, n = 24; NC, normal control, n = 12. **p < 0.005; ***p < 0.001; ****p < 0.0001).

RSAD2 is elevated in the skin of patients with DM

As we can see in Figure 5, RSAD2 is at the core condition. Therefore, RSAD2 was selected for immunohistochemistry for validation in skin tissues of patients with dermatomyositis. It showed that RSAD2 expression was significantly increased in patients with dermatomyositis compared with normal controls (Figure 7). It was mainly expressed in skin fibroblasts.

Figure 7. RSAD2 is elevated in skin fibroblasts of patients with DM. (A:Representative immunohistological staining of RSAD2 in the skin of DM and NC patients. B. Statistical analysis of RSAD2 immunohistochemistry. DM, dermatomyositis, n = 3; NC, normal control, n = 3. Scale bar = 50 μm.).

Discussions

Dermatomyositis is an inflammatory autoimmune disease, affecting skin, muscle, and blood vessels. The pathogenesis is complicated, involved in multiple factors, including genetic, epigenetic, and environmental factors [28]. This study was expected to find key genes and signal pathway in the skin of DM. A total of 380 DEGs of GSE46239 were identified, consisting of 271 up-regulated genes and 109 down-regulated genes. And we found that the skin of DM patients and the skin of COVID-19 patients shared a common transcriptional profile. Then, we used R package cluster Profiler to perform GO and KEGG analysis for the DEGs of GSE46239. A PPI network of DEGs in GSE46239 was constructed and 173 genes and 829 edges, and two significant modules were chosen from the PPI network. Next, 10 hub genes were identified by using plug-in cytoHubba, namely MX1, ISG15, IFIT3, IFIT1, RSAD2, IFIT2, IFI6, XAF1, IRF9, MX2. Finally, we verified these genes in the muscle and skin tissue of patients with DM.

Interferon is a cytokine produced by monocytes and lymphocytes, which has the effects of anti-viral, anti-tumor, anti-infection and immune regulation. Interferon (IFN) family can be divided into three classes, consisting of type I IFN, type II IFN, and type III IFN [29]. The etiology of DM is still unclear. However, type I IFN signal pathway is strongly associated with the pathogenesis of DM [30–32]. DM patients showed an increased expression of the type I IFN-inducible genes in blood, muscle and skin [33–36]. In addition, type I IFN-inducible genes were positively correlated with the disease activity of DM [33, 35]. Type I IFN mainly includes IFNαand β, produced by macrophages and plasmacytoid dendritic cells (pDCs). IFN can activate JAK-STAT pathway by phosphorylation of STAT 1 and 2 via binding with the receptors [29].

The involvement of RSAD2 in the pathogenesis of DM has not been reported. As an IFN-inducible gene, RSAD2 played an important role in other autoimmune diseases, such as systemic lupus erythematosus (SLE) [37–39], systemic sclerosis (SSc) [40, 41], and primary Sjogren syndrome (pSS) [42, 43]. RSAD2 was overexpressed in CD19⁺ B cell of pSS and knockdown of RSAD2 can attenuate B cell hyperactivity through suppressing NF-κb signaling pathway [42]. Besides, the expression of RSAD2 in fibroblasts was significantly increased after virus stimulation [44]. In this study, we found that skin lesions in DM were closely related to the viral infection pathway and that RSAD2 was upregulated mainly in fibroblasts of DM skin. Taken together, RSAD2 may participate in the pathogenesis of DM by affecting fibroblasts in the skin of DM. In the peripheral circulation, RSAD2 may enhance the immune response and the production of autoantibodies by enhancing the reactivity of immune cells such as B cells in DM patients. This finding was interesting and the specific mechanism needs to be further investigated.

Recent studies suggest that the type I interferon pathway plays an important role in COVID-19 [45–47]. Activation of the cGAS-STING pathway associated with type I interferon is involved in lung damage in COVID-19 patients [17]. And type I interferon exacerbates inflammation in patients with severe COVID-19 [48]. Our results suggest that the skin of patients with dermatomyositis and the skin of patients with COVID-19 have similar transcriptional profiles. And the DEGs in the skin of patients with dermatomyositis are closely related to type I interferon signaling pathway and COVID-19 signaling pathway. It shows that DM is as unfathomable as COVID-19, and the pathogenesis is similar but complex.

In summary, we found that patients with DM shared a common transcriptional profile with COVID-19. And we identified 380 DEGs and 10 hub-genes in DM, which were associated with type I IFN signaling pathway and COVID-19 closely. In-depth research on these hub-genes will help us better understand the pathogenesis of DM. These hub-genes might also become new targets for the diagnosis and treatment of DM.

Authors’ contributions

Yiying Yang conceived and designed the study, performed the experiments, analyzed the data and drafted this manuscript, as the first author. Hongjun Zhao, and Jie Song performed the experiments. Muyao Guo and Huali Zhang contributed to analysis of the data, and revising of the manuscript. All authors have read and approved the final manuscript.

Data availability

The data used to support the findings of this study is available from GSE46239 and GSE193068.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China (Grant No. 81901664, 82070018, 82270097), the Natural Science Foundation of Hunan Province (2020JJ5894).