The transcriptome expression levels related to ovulation induction and acupuncture protection therapy in rats through gene microarray

ABSTRACT

This study aimed to explore the potential molecular mechanisms of acupuncture in the adjuvant treatment of infertility. A rat model of ovulation induction was constructed using horse serum gonadotropin and human chorionic gonadotropin, and then acupuncture was used to treat the model rats at Guanyuan and Sanyinjiao points. Hematoxylin eosin (HE) staining, the endometrial thickness, and enlargement of uterine cavity were determined to evaluate the effects of acupuncture on the uterus of ovulation induction rats. Meanwhile, the uterus tissues were sent for gene microarray analyses. Acupuncture enhanced the endometrial receptivity of ovulation induction rats. Gene microarray showed that 189 overlapped differential expressed genes (DEGs) were identified, and these overlapped DEGs were divided into four clusters by Mfuzz algorithm. Afterwards, interaction networks containing 842 interaction pairs were established, and the DEGs in the interaction networks were significantly enriched in 20 BP and 10 KEGG pathways, including wound healing, response to glucocorticoid, TNF signaling pathway, terpenoid backbone biosynthesis, and glutathione metabolism. By calculating topological parameters, Anxa5, Casp8, Tgm2, Frk, Mmp7, Timp1, Hmgcr, Cth, Serpinal, and Abcbla were the important hub genes of the interaction networks. Our findings revealed the changes at transcriptome levels related to ovulation induction and acupuncture protection therapy.

KEYWORDS:

• Infertility
• ovulation induction
• acupuncture
• transcriptome

Introduction

The World Health Organization (WHO) defines infertility as a disease of the reproductive system, characterized by the failure to clinically conceive after 12 months or more of regular unprotected sexual intercourse (Boltz et al. 2017). Infertility is a major public health problem that has not received attention, and affects about 12%–15% of childbearing couples worldwide (McLaren 2012). Infertility remains a social burden on women (Inhorn and Patrizio 2015), as well as lifestyle factors (e.g. diet, condition, and behavior) and organic diseases (such as fallopian tube disease, pelvic adhesion, ovulation disorders, and endometriosis) are the primary cause of female infertility (Homan et al. 2007). Nowadays, the conventional treatment for female infertility is primarily drug therapy, including human chorionic gonadotropin (HCG), clomiphene, follicle-stimulating hormone, gonadotropin-releasing hormone analogue, and melbine (Stanford et al. 2013; Weiss and Clapauch 2014). However, the use of these drugs may have some potential side effects, such as an increased risk of cancer and unknown pharmacological effects (Weeg et al. 2012; Diergaarde and Kurta 2014). Therefore, there is an urgent need to search and develop novel therapeutic strategies to assist the treatment of female infertility.

With the continuous development of science and technology, assisted reproductive technology (ART) including in vitro fertilization and artificial insemination has been used to treat female infertility (Lindsay and Vitrikas 2015). Although ART can help women to improve pregnant rate, there are some undesirable effects to the technology. On the one hand, the live birth rate through ART decreased with the maternal age and pregnancy period (Abuzeid et al. 2014), and on the other hand, ART may have the following adverse effects: antepartum hemorrhage, premature delivery, cesarean, congenital abnormality, low body-weight newborn, gestational diabetes, and perinatal mortality rate (Pandey et al. 2012). Due to the problems of cost, low live birth rates and adverse effects, an increasing number of infertile couples have considered to use complementary therapies (Miner et al. 2018).

Acupuncture is one of the important methods of complementary therapy, which can be combined with traditional drugs to achieve better effects. Acupuncture is an indispensable part of traditional Chinese medicine, which involves inserting needles into the skin along the meridians to restore the normal circulation of the veins and achieve the function of treating patients (Yun et al. 2019). A systematic review and meta-analysis of Yun et al. (2019) showed that acupuncture and its combination therapy may be effective in the treatment of female infertility. Wang et al. (2018) used acupuncture combined with cupping to treat an infertile woman, and found that after 28 treatments, the woman was successfully pregnant, suggesting that acupuncture can help treat infertility. Another study has demonstrated that acupuncture is a non-drug therapy with minimal side effects, and is beneficial in increasing the success rate of in vitro fertilization (Djaali et al. 2019). Based on these reports, acupuncture can be used as an adjunct treatment for infertility. However, the molecular mechanisms of acupuncture in the adjuvant treatment of infertility remain unclear.

In this study, a rat model of ovulation induction was constructed, and then acupuncture was used to treat the model rats at Guanyuan and Sanyinjiao points. Afterwards, the underlying molecular mechanisms of acupuncture protection therapy were determined using gene microarray. These results would improve our understanding of acupuncture in treatment of infertility, and may provide new insights and therapeutic targets for infertility therapy.

Materials and methods

Grouping and treatment of experimental animals

Eighteen SPF female Sprague Dawley (SD) rats weighing 250 ± 20 g, and 18 SPF male SD rats weighing 550 ± 20 g were obtained from Shanghai Jiesijie Experimental Animal Co., Ltd. (Shanghai, China). All rats were fed under controlled temperature (24 ± 2°C) and humidity (50 ± 5%) conditions, with a 12 h light/dark cycle. The rats were free access to food and water during the experiment. After acclimatization for seven days, the 18 female rats were randomly and equally divided into three groups (n = 6 for each group): control group, model group, and treatment group. The female rats in the model and treatments groups were used to establish an ovulation induction model, and the rats in the control group were without treatment. The female rats in the treatment group before modeling were electro-acupunctured at Guanyuan and Sanyinjiao points for 15 min each time, and for seven consecutive days; and after treatment, the rats were used for modeling (continue to acupuncture).

The establishment method of an ovulation induction model was shown as follows: the female rats were injected with pregnant horse serum gonadotropin (PMSG, 0.5U/rat, Ningbo Second Biological Pharmaceutical Factory, Zhejiang, China) intraperitoneally, and then after 48 h, the female rats were injected with the same amount of HCG (Ningbo Second Biological Pharmaceutical Factory). After that, the rats in the groups were caged with male rats in a ratio of 1:1. All the animal experiments were conducted in accordance with the National Medical Advisory Committee (NMAC) guidelines using approved procedures of the Institutional Animal Care and Use committee of Nanjing University of Chinese Medicine.

Samples collection and hematoxylin eosin (HE) staining

After the experiments, all female rats were killed by cervical dislocation, and the uterus tissues of all female rats were collected for the measurement of endometrial thickness, and the enlargement of uterine cavity. Furthermore, the midsections on the same side of the uterus tissues in all the female rats were used for HE staining, and the rest uterus tissues were sent to Yanzai Biotechnology (Shanghai) Co., Ltd. (Shanghai, China) for gene microarray analyses.

The methods of HE staining were described as previously (Feldman and Wolfe 2014). Briefly, the uterus tissues were washed by PBS, fixed in 4% paraformaldehyde, and then embedded in paraffin. After that, the 5-µm sections were cut, and stained with HE, and the images were taken under a microscope (Olympus Corporation, Tokyo, Japan) at 40× magnification.

Screen of differentially expressed genes (DEGs) by gene microarray data

The sequencing data were deposited in a public database (NCBI GEO database), and the number is GSE190472. The samples were divided into two comparison groups: control vs. model, and control vs. treatment. DEGs between control and model groups as well as control and treatment groups were screened using Limma version 3.34.0 package in R 3.6.1 (https://bioconductor.org/packages/release/bioc/html/limma.html) (Ritchie et al. 2015). The threshold values for selection of DEGs were |log₂ fold change (FC)|>1 and false discovery rate (FDR) <0.05. Based on the obtained DEGs, pheatmap package version 1.0.8 in R (https://cran.r-project.org/package=pheatmap) (Wang et al. 2014) was used to perform bidirectional hierarchical clustering analysis (Eisen et al. 1998) based on Euclidean distance of expression values of DEGs (Mostaco-Guidolin et al. 2017). Afterwards, the DEGs obtained from each comparison group were compared, and the overlapped DEGs were selected for further study. These overlapped DEGs were then submitted for biological process (BP) of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses using DAVID version 6.8 (https://david.ncifcrf.gov/) (Huang da et al. 2009a, 2009b). P-value less than .05 was selected as the threshold for enrichment significance.

Screen of DEGs in the same expression pattern

Mfuzz version 2.42.0 (http://www.bioconductor.org/packages/release/bioc/html/Mfuzz.html) (Kumar and Futschik 2007) in R was used to analyze the expression patterns of the overlapped DEGs, and the gene clustering of expression modules was obtained to observe the expression change tendency of the target genes in the development process of different groups.

Construction and analysis of interaction networks

STRING database (version 10.0, http://string-db.org/) (Szklarczyk et al. 2017) was used to search the interaction relationships between proteins of the aforementioned DEGs, and interaction networks were constructed. The interaction networks were visualized using Cytoscape version 3.6.1 (http://www.cytoscape.org/) (Shannon et al. 2003), and the DEGs in the interaction networks were further submitted for BP of GO terms and KEGG pathways enrichment analyses (Huang da et al. 2009a, 2009b).

After that, we analyzed the topology of the interaction networks, and calculated the four important network topology parameters: degree, betweenness centrality (BC), closeness centrality (CC), and path length (Jeong et al. 2001; Brown et al. 2004).

Results

Rat uterus morphology and endometrial thickness

The morphology of rat uterus was observed to determine the protective effects of acupuncture on the uterus of ovulation induction rats. It is clear that the rat uterus in the control group was normal. In the model group, the rat uterine cavity was full without obvious enlargement; the endometrium was thickened; no obvious changes were observed in the epithelial cells; and the interstitial cells were closely arranged without obvious decidualization. However, in the treatment group, acupuncture could enlarge the uterine cavity and thicken the wall and the endometrium; as well as the epithelial cells were columnar; and the interstitial cells were loosely arranged, with large cells, and decidua-likes changes (Figure 1(a)).

Figure 1. The protective effects of acupuncture on the uterus of ovulation induction rats. (a) The uterus morphology in the rats of the control, model, and treatment groups under a magnification of 40×. (b) The endometrial thickness of the rats with different treatments. (c) The uterine cavity percentage of the rats with different treatments. *: P < .05, compared with the control group; ^#: P < .05, compared with the model group.

After that, the endometrial thickness and enlargement of uterine cavity were further measured. Compared with the control rats, the endometrial thickness in the induction model rats was significantly decreased (P < .05); while the thickness in the acupuncture-treated rats was evidently higher than that in the model rats (P < .05, Figure 1(b)). Besides, the uterine cavity percentage in the control, model, and treatment groups were 15.23 ± 4.17%, 4.67 ± 1.72%, and 7.24 ± 1.19%, respectively (Figure 1(c)). These showed that after the establishment of the induction model, the uterine cavity percentage was significantly reduced compare to the control rats (P < .05); and after acupuncture treatment, the uterine cavity percentage was increased, but had no significant difference between the model and treatment groups (P > .05, Figure 2(c)). Taken together, all the results indicated that acupuncture could improve the endometrial receptivity of ovulation induction rats.

Figure 2. The results of differentially expressed genes (DEGs) between model and control groups, as well as between treatment and control groups. (a) The volcano (left) and heatmap (right) figures of top 15 down-regulated and up-regulated DEGs between model and control groups. (b) The volcano (left) and heatmap (right) figures of top 15 down-regulated and up-regulated DEGs between treatment and control groups. The blue and orange dots represent the down-regulated and up-regulated genes, respectively. The black horizontal line indicates FDR >0.05. The two vertical lines mean |log₂ fold change| > 1. (c) Venn diagram of DEGs by comparing the DEGs between model and control groups, and between treatment and control groups.

Screen of DEGs and functional analyses

After comparison, a total of 446 and 404 DEGs were, respectively, screened between model and control groups, and treatment and control groups based on the thresholds of |log₂ FC| > 1 and FDR < 0.05 (Figure 2(a,b)). The bidirectional hierarchical clustering analysis showed that the expression values of these DEGs could well distinguish the model and control groups, as well as the treatment and control groups, and suggested that the screened DERs were characteristically expressive (Figure 2(a,b)). Afterwards, by comparing the DEGs between model and control groups, and treatment and control groups, a total of 189 overlapped DEGs were obtained (Figure 2(c)). The obtained overlapped DEGs were then submitted for BP of GO terms and KEGG pathways analyses. It was found that the 189 overlapped DEGs were significantly enriched in 15 BP of GO terms and 9 KEGG pathways. As shown in Figure 3(a), in the involved BP, these overlapped DEGs were related to negative regulation of apoptotic process, wound healing, aging, response to glucocorticoid, protein hormooligomerization, response to lipopolysaccharide, response to estradiol, extracellular matrix organization, hippocampus development, circadian rhythm, response to nutrient, female pregnancy, proteolysis, cell adhesion and response to hypoxia. Furthermore, these overlapped DEGs also played important roles in vasopressin-regulated water reabsorption, thyroid hormone synthesis, TNF signaling pathway, complement and coagulation cascades, biosynthesis of antibiotics, metabolic pathways, terpenoid backbone biosynthesis, propanoate metabolism, and glutathione metabolism (Figure 3(b)).

Figure 3. Functional analyses of the overlapped DEGs. (a) Biological process (BP) of these overlapped DEGs. (b) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment of the overlapped DEGs. The horizontal axis represents gene count, and the vertical axis represents the name of items. Additionally, the color represents significance. The closer the color is to red, the more significant it is.

Selection of DEGs in the same expression pattern using Mfuzz algorithm

We used Mfuzz algorithm to analyze the expression patterns of the overlapped DEGs, and observe the expression change tendency of the target genes under different treatments. It was found that the above overlapped DEGs could be clustered into four expression patterns: cluster 1, cluster 2, cluster 3, and cluster 4 (Figure 4). In the cluster 1, the expression levels of DEGs were continuously decreased from the ovulation induction (model) to acupuncture protective ovulation induction (treatment). The expression levels of DEGs in the cluster 2 were steadily increased from the ovulation induction (model) to acupuncture protective ovulation induction (treatment). However, the expression levels of DEGs were firstly decreased and then increased in the cluster 3 from ovulation induction (model) to acupuncture protective ovulation induction (treatment); whereas in the cluster 4, the expression levels of DEGs were firstly increased and then decreased. The DEGs in each cluster were displayed in Supplementary Table S1.

Figure 4. Clustering diagram of expression patterns of the overlapping DEGs based on Mfuzz algorithm.

Analyses of the interaction networks

We used STRING database to search the interaction relationships between proteins of the aforementioned overlapped DEGs, and obtained 842 interaction pairs (842 gene nodes) through retaining the interaction pairs with an interaction score greater than 0.4 (Figure 5). After that, the DEGs in the interaction networks were further sent for BP and KEGG pathways enrichment analyses. These DEGs in the interaction networks were significantly enriched in the 20 BP terms (Figure 6(a)), including the aforementioned 15 BP terms, and other five terms (defense response to bacterium, response to drug, heart development, actin filament organization, and response to peptide hormone). In addition, these DEGs in the interaction networks were also associated with 10 significant KEGG pathways, including the aforementioned nine KEGG pathways, and toxoplasmosis (Figure 6(b)).

Figure 5. The interaction networks of the overlapped DEGs product proteins. Different colors represent different clusters, and the size of nodes indicates the degree in the networks.

Figure 6. Functional analyses of the DEGs in the interaction networks. (a) BP of these overlapped DEGs. (b) KEGG pathways enrichment of the overlapped DEGs.

Topology analysis of the interaction networks

Finally, we analyzed the network topology parameters of the interaction networks, including degree, BC, CC, and path length. According to the degree, the top 10 hub genes were Anxa5, Casp8, Tgm2, Frk, Mmp7, Timp1, Hmgcr, Cth, Serpinal, and Abcbla. The specific information of the top 10 hub genes were shown in Table 1. It is clear that Anxa5, Frk, Mmp7, and Timp1 belonged to cluster 1; Cth and Serpinal belonged to cluster 2; Casp8, Tgm2, and Hmgcr belonged to cluster 3; as well as Abcbla belonged to cluster 4.

Table 1. The four important network topology parameters of link gene nodes (Top 10).

	Average shortest	Betweenness	Closeness
Symbol	Path length	Centrality	Centrality	Degree	Cluster
Anxa5	1.98830	0.05790	0.50294	33	1
Casp8	2.05263	0.05912	0.48718	29	3
Tgm2	2.14035	0.02563	0.46721	25	3
Frk	2.11696	0.05971	0.47238	25	1
Mmp7	2.22807	0.01714	0.44882	24	1
Timp1	2.09942	0.02007	0.47632	24	1
Hmgcr	2.17544	0.03782	0.45968	23	3
Cth	2.21637	0.03465	0.45119	23	2
Serpina1	2.19883	0.03017	0.45479	23	2
Abcb1a	2.22222	0.02803	0.45000	22	4

Discussion

Infertility is a global health problem, seriously affecting mental and physical health of infertile couples (Rouchou 2013). Acupuncture is well tolerated, and has been reported to be a part of an effective treatment for infertility (Wang et al. 2018). However, the potential molecular mechanisms are still unclear. This study found that acupuncture could improve the endometrial receptivity of ovulation induction model rats, and may improve the success rate. Gene microarray analysis showed that a total of 446 and 404 DEGs were, respectively, screened between model and control groups, and treatment and control groups. After comparison, 189 overlapped DEGs were obtained, and were significantly enriched in 15 BP and 9 KEGG pathways. By Mfuzz algorithm, these overlapped DEGs were clustered into four clusters. After that, interaction networks containing 842 interaction pairs were proposed, and the DEGs in the interaction networks were significantly enriched in 20 BP and 10 KEGG pathways, including wound healing, response to glucocorticoid, TNF signaling pathway, complement and coagulation cascades, terpenoid backbone biosynthesis, and glutathione metabolism. By analyzing degree, BC, CC, and path length, Anxa5, Casp8, Tgm2, Frk, Mmp7, Timp1, Hmgcr, Cth, Serpinal, and Abcbla were closely related to acupuncture protection therapy. This finding was the first time to reveal the changes at the transcriptome level associated with ovulation induction and acupuncture protection.

Previous studies have reported that acupuncture has positive effects on infertility treatment (Djaali et al. 2019; Yun et al. 2019). The follicles are essential for maintaining the endocrine function of the female gonads and fertility (Sanchez et al. 2017). Therefore, sustaining or improving the quality of oocytes, and improving the survival rate of ova can help to maintain female ovarian function, thus enhancing the success rate of pregnancy and treating infertility (Tilly and Sinclair 2013). A meta-analysis of Manheimer et al. (2008) showed that acupuncture combined with embryo transfer can improve the survival rate of ova, pregnancy, and live birth rates in women undergoing in vitro fertilization. Furthermore, a previous case report reported that acupuncture could have protective effects on infertility through improving sperm quality and ova survival, as well as balancing the endocrine system and hormones (Zhu et al. 2018). Our study established an ovulation induction model, and found that acupuncture enhanced the endometrial receptivity and may promote survival rate, which were in consistence with previous findings (Zhong et al. 2019). All these findings indicated that acupuncture can be used as an adjunct treatment for infertility.

Further to explore the molecular mechanisms of acupuncture affecting endometrial receptivity, gene microarray was performed. It was found that Anxa5, Casp8, Tgm2, Frk, Mmp7, Timp1, Hmgcr, Cth, Serpinal, and Abcbla were closely related to acupuncture protection therapy. By Mfuzz algorithm, Anxa5, Frk, Mmp7, and Timp1 were down-regulated in the ovulation induction model, and unceasingly decreased after acupuncture; while Cth and Serpinal were continuously increased in the ovulation induction model, and after acupuncture treatment. For Casp8, Tgm2, and Hmgcr, their levels were firstly reduced in the ovulation induction model, and then increased after acupuncture; whereas the tendency of Abcbla level was opposite to Casp8, Tgm2, and Hmgcr level.

Anxa5 has been reported to be highly expressed in normal placenta, and is considered to have anticoagulant properties (Krikun et al. 1994). A study of Aranda et al. (2018) showed that decreased placental Anxa5 expression may contribute to the development of vascular complications in pregnancy, such as intrauterine growth restriction and preeclampsia. Frk, a non-receptor tyrosine kinase, has been reported to play an important role in many tumors, such as breast cancer (Ogunbolude et al. 2017) and non-small-cell lung cancer (Zhang et al. 2020). Daigle et al. (2019) demonstrated that ScFRK1, ScFRK 2, and ScFRK 3 participated in the gametophyte development of a wild potato species Solanum chacoense. Mmp7, a kind of matrix metalloproteinases, can degrade all components of the basement membrane and the extracellular matrix. Timp1 is usually secreted by the same type of cells that secrete Mmp9, and is a tissue suppressor corresponding to MMP-9 (Zhang et al. 2019). A previous study indicated that Mmp7 and Mmp9 were the key mediators of trophoblast cell invasion into decidua region and angiogenesis during placenta formation (Hamutoglu et al. 2020). Another study reported that Cth inhibition in pregnant mice could decrease plasma hydrogen sulfide levels and increase blood pressure/placental fetal vascular branching irregularities (Akahoshi et al. 2019), which indicated that Cth may be important for infertility treatment. Casp8 is related to apoptosis, and its methylation has been demonstrated to be associated with the increased risk for neural tube defects of fetuses (Huang et al. 2019). Tgm2 has been reported to be highly expressed in the mid-trimester amniotic fluid of pregnant women with systemic lupus erythematosus (Jeon et al. 2020). In addition, Serpinal, Hmgcr, and Abcbla, the other three genes, have not been reported in pregnancy and infertility. When these studies combined with our results, it can be inferred that the changes of Anxa5, Casp8, Tgm2, Frk, Mmp7, Timp1, Hmgcr, Cth, Serpinal, and Abcbla expressions may be involved in the ovulation induction and acupuncture protection.

Further, these DEGs were submitted for functional analyses, and were found to be significantly enriched in 20 BP and 10 KEGG pathways, including wound healing, complement and coagulation cascades, response to glucocorticoid, TNF signaling pathway, terpenoid backbone biosynthesis, and glutathione metabolism. It has been reported that angiogenesis, the formation of new blood vessels from pre-existing vascular beds, occurs naturally during wound healing, a woman’s menstrual cycle and pregnancy (Rizzi et al. 2017). Qiao et al. (2016) found that complement and coagulation cascades were related to pregnancy loss through whole exome sequencing. Another study demonstrated that high levels of progesterone during pregnancy could induce selective T cell death by binding to glucocorticoid receptors, which suggested that response to glucocorticoid may participate in pregnancy process (Hierweger et al. 2019). Terpenoid backbone biosynthesis was reported to play an important role in the lipid metabolism of obese boys (Liang et al. 2020). TNF signaling pathway is associated with inflammatory response, and increased TNF-α level could intensify the inflammatory response in trophoblast cells, thereby affecting pregnancy progression (Huang et al. 2018). Besides, glutathione, the mother of all antioxidants, can be directly involved in the neutralization of free radicals and reactive oxygen species, as well as plays crucial roles in infertility and adverse pregnancies (Adeoye et al. 2018). Taken together, we speculated that wound healing, complement and coagulation cascades, response to glucocorticoid, TNF signaling pathway, terpenoid backbone biosynthesis, and glutathione metabolism may be closely related to ovulation induction protected by acupuncture treatment. However, the specific effects of these pathways after acupuncture treatment should be further investigated.

Conclusions

In conclusion, acupuncture can improve the endometrial receptivity and may promote survival rate of pregnancy. By gene microarray analyses, Anxa5, Casp8, Tgm2, Frk, Mmp7, Timp1, Hmgcr, Cth, Serpinal, and Abcbla expressions were found to be the key genes in acupuncture protective treatment. In addition, some pathways like wound healing, complement and coagulation cascades, response to glucocorticoid, TNF signaling pathway, terpenoid backbone biosynthesis, and glutathione metabolism, may play essential roles in ovulation induction protected by acupuncture treatment. However, the survival rate of ova and the ovulation-related data after acupuncture treatment need to be investigated in the further experiments. Our findings reveal the changes at transcriptome levels related to ovulation induction and acupuncture protection therapy, and provide a basis for Anxa5, Casp8, Tgm2, Frk, Mmp7, Timp1, Hmgcr, Cth, Serpinal, and Abcbla as potential targets in infertility therapy.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available in NCBI GEO database at https://www.ncbi.nlm.nih.gov/geo/info/linking.html, reference number GSE190472 and in Science Data Bank at http://doi.org/10.57760/sciencedb.02017.

Ethics approval

All the animal experiments were conducted in accordance with the National Medical Advisory Committee (NMAC) guidelines using approved procedures of the Institutional Animal Care and Use committee of Nanjing University of Chinese Medicine.

Additional information

Funding

This study was supported by National Natural Science Foundation of China (No. 81873371), and Natural Science Foundation of Jiangsu Province (No. BK20181420).