Integrated analysis of long noncoding RNA interactions reveals the potential role in progression of human papillary thyroid cancer

Abstract Recent scientific evidence has suggested that long noncoding RNAs (lncRNAs) play an important part in tumorigenesis as an important member of competing endogenous RNAs (ceRNAs). Hundreds of RNA sequence data and relevant clinic information are freely accessible in The Cancer Genome Atlas (TCGA) datasets. However, the role of cancer‐related lncRNAs in papillary thyroid cancer (PTC) is not fully understood yet. In this study, we identified 461 RNA sequencing data from TCGA. Subsequently, 45 lncRNAs, 21 miRNAs, and 78 mRNAs were chosen to construct a ceRNA network of PTC. Then, we analyzed the correlation between these 45 PTC‐specific lncRNAs and clinic features and patient outcome. Thirty‐seven of these lncRNAs were found to be closely related to age, race, gender, lymph node metastasis, TNM staging system, and patient outcome. Additionally, three of them were linked to PTC patient overall survival. Eventually, we selected eight lncRNAs randomly and performed quantificational real‐time polymerase chain reaction (qRT‐PCR) in 28 newly diagnosed patients with PTC to verify the reliability of the above results. The results of qRT‐PCR are totally in agreement with the bioinformatics analysis. Additionally, it was found that HAND2‐AS1 was negatively related to tumor size (P < 0.05). The results were consistent with the bioinformatics analysis in TCGA. Taken together, we identified the differentially expressed lncRNAs and constructed a PTC ceRNA network. The study provides a new perspective and supplement for our understanding of lncRNAs in PTC development and reveals potential diagnostic and prognostic markers in PTC.


| INTRODUCTION
Thyroid cancer remains the most prevalent malignant endocrine disorder with low mortality but increased incidence and recurrence. According to the United States cancer statistics, approximately 726 646 people lived with thyroid cancer in 2014, and there would be 56 870 new cases diagnosed as thyroid cancer in 2017. 1 Moreover, by 2030, thyroid cancer was expected to be the fourth leading cancer diagnosis. 2 Papillary thyroid cancer (PTC) is the most common pathological type, accounting for 80% of thyroid cancer. In the long run, the prognosis of PTC is excellent and the 5-year survival rate could reach 98.2%; however, its recurrence rate is relatively high. 3,4 On the other hand, anaplastic thyroid cancer (ATC) is less frequent even if it has more aggressive biological behavior and worse prognosis. Recent studies showed that genetic and epigenetic alterations were involved in PTC pathogenesis and could be markers for prognosis, such as BRAF mutations, 5-7 RAS mutations, 8,9 and RET rearrangements. 10,11 LncRNAs are described as transcripts with longer than 200 nucleotides without translation function. It plays an essential part in the regulation of gene transcription, posttranscription, and epigenetic modulation. [12][13][14] Amounting evidence indicate that lncRNAs are crucial in oncogenesis and tumor development. [15][16][17] It was also found that plenty of lncRNAs was differentially expressed in thyroid cancer compared with normal adjacent tissues. These lncRNAs were expected to be correlated with diagnosis and prognosis of PTC. 18 The competing endogenous RNA (ceRNA) hypothesis was first proposed about how messenger RNAs (mRNAs), transcribed pseudogenes, and lncRNAs connect to each other using microRNA response elements (MREs). 19 Recently, ceRNA hypothesis was proved to be essential in physiologic and pathological conditions such as cancer. 20 Thousands of lncRNAs were found aberrantly expressed in PTC samples compared with normal thyroid tissues using genomewide analysis. 21 With further research, more and more lncRNAs were found to be closely related to the initiation and development of PTC. For instance, 218 aberrantly expressed lncRNAs were identified in genomewide expression screening, and two lncRNAs (XLOC-051122 and XLOC-006074) were found significantly related to lymphatic metastasis. 22 Another study performed by Li et al 23 described four lncRNAs signature (RP11-536N17.1, RP11-508M8.1, AC026150.8, and CTD-2139B15.2) predicting PTC prognosis. Differentially expressed RNAs in 348 cases were identified in three PTC variants and were selected to construct the ceRNA network by Yanjing Zhao. 24 However, the ceRNA network of PTC is still poorly investigated in a larger cohort, and the validation analysis needs to be performed. In this study, dysregulated RNAs (including lncRNAs, miRNAs, and mRNAs) were detected in different stages and lymph node status, and these RNAs were used to construct the ceRNA network of PTC. Next, quantificational real-time polymerase chain reaction (qRT-PCR) validation was performed to confirm the credibility of bioinformatics analysis results. GO terms and KEGG pathways were performed to identify the differentially expressed mRNAs. Our results indicate that ceRNA network plays an important part in the development of PTC, and the study also provides a novel perspective for the better understanding of lncRNA-miRNA-mRNA interactions in PTC initiation and development.

| RNA sequence data analysis
Firstly, we obtained all RNA expression profiles (level 3) of patients with PTC from the TCGA data portal (up to May 05, 2016). The raw data of lncRNAs and mRNAs were processed and normalized by TCGA RNASeqv2 system. The PTC level 3 miRNA data were normalized and provided by Illumina HiSeq 2000 microRNA sequencing platforms (Illumina Inc., Hayward, CA, USA). No further normalization was needed as all the RNA sequencing data have already been normalized already. Then, all samples collected in TCGA were divided into three groups to perform the differential analysis of lncRNAs, miRNAs, and mRNAs. Thus, all collected samples were grouped as: (a) PTC tumor samples and normal thyroid samples; (b) PTC samples with lymphatic metastasis and PTC samples without lymphatic metastasis; (c) PTC samples with stages I-II and PTC samples with stages III-IV. Four groups were listed as follows: stages I-II with/without lymphatic metastasis and stages III-IV with/without lymphatic metastasis, and these groups contained 112, 195, 107, and 47 cases, respectively. The intersected lncRNAs, miR-NAs, and mRNAs were then chosen to construct the ceRNA network which shows the genetic regulation related to PTC. The bioinformatics analysis is illustrated in the flowchart ( Figure 1).

| Gene functional enrichment analysis
The intersected mRNAs were selected following the above methods and were analyzed to detect the potential biologic functions and pathways. The Database for Annotation, Visualization, and Integrated Discovery (DAVID) Bioinformatics resources (http://david.abcc.ncifcrf.gov/) were employed to analyze the potential functions and pathways. The biological functions and pathways were considered significant when an enrichment score of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis exceeds 1.5 (P < 0.05).

| The ceRNA network in PTC
The interactions among lncRNAs, miRNAs, and mRNAs contribute to the post-transcriptional regulation of mRNAs. The ceRNA network is built based on miRNA sponge 25 in our study. This theory indicates that lncRNAs conjunct with miRNAs to regulate mRNAs expression. Furthermore, the miRNA-lncRNA and miRNA-mRNA relationships were collected in the starBase v2.0 database (http://starbase.sysu. edu.cn/), 26  lncRNAs, miRNAs, and mRNAs with fold change (FC) >2.0 or FC <0.5 and P < 0.05 were selected. Subsequently, mi-Randa tools (http://www.microrna.org/microrna/home.do), the miRNAs target prediction tool, were employed to detect the lncRNA-miRNA interactions. Then, the mRNA-miRNA interactions were analyzed by miRTarBase (http://mirtarbase. mbc.nctu.edu.tw/) and TargetScan (http://www.targetscan. org/). The summary results were identified to determine the specific interacted lncRNAs and mRNAs. Then, according to ceRNA theory, the negatively regulated lncRNA-miRNA and mRNA-miRNA pairs were analyzed in the ceRNA network in PTC. 27 Finally, the ceRNA network was built using Cytoscape v3.0. 28 The ceRNA network construction is illustrated in a flowchart ( Figure 2).

PCR verification
Specific lncRNAs involved in ceRNA network and bioinformatics analysis were selected to assess the clinical characteristics including age, race, gender, lymphatic metastasis, TNM staging system, and patient outcome in TCGA. In addition, qRT-PCR validation of eight randomly selected lncRNAs in 28 paired frozen tissues was performed to examine the credibility of the bioinformatics analysis results.
Following the manufacturer's protocol, total RNAs were extracted and purified using TRIzol reagent (Pufei, Shanghai, China). Then, concentration and purity of total RNA were detected by the spectrophotometer from BioTek (Shoreline, WA, USA). According to the manufacturer's protocol, a total of 1000 ng RNA was used to perform cDNA synthesis, and reverse transcription reaction was conducted using PrimeScript™ from Takara (Dalian, China). Candidate ln-cRNA expression was amplified with TAKARA SYBR ® Premix Ex Taq™ (Dalian, China) on an Applied Biosystems ViiA™ 7 System (Foster City, CA, USA) using specific primers in duplicate. All primers were designed and synthesized by Generay Biotech (Shanghai, China), and the sequences of selected lncRNAs primers were presented in Table 2. qPCR primers were determined to anneal at 60°C according to the instructions. The comparative CT method was used to measure the relative expression of candidate lncRNAs. 29 β-actin is used as the internal control in our experiment. FC values were shown as 2 −∆∆Ct , where ∆∆Ct = (Ct RNAs − Ct β-actin ) tumor − (Ct RNAs − Ct β-actin ) adjacent normal samples .

| Statistical analysis
The Student's t test was performed to detect the statistical significance for two groups, and ANOVA for multiple groups using the GraphPad Prism 6.0. All data are expressed as mean ± SD (standard deviation). The threshold of false discovery rate (FDR) was set as 0.001 to reduce the false positives in multiple tests in bioinformatics analysis, and the threshold of P value was set as 0.05 to evaluate the null hypothesis. Box plots were used to present the FC and clinical relevance of the specific lncRNAs in donor samples. The upper and lower boundaries of the boxplots are at 25th and 75th percentiles. The vertical lines below and above the box

| PTC-specific lncRNAs
Totally, 1806 lncRNAs were identified from TCGA database (Data S1). We compared the expression of lncR-NAs in 461 patients with PTC tumor tissues with adjacent normal tissues and identified 199 differentially expressed lncRNAs (FC > 2, P < 0.01). Further analysis was carried out between tumor tissues and adjacent normal thyroid tissues on tumor stages and lymph node status. For patients without lymphatic metastasis, there are 99 aberrant expression lncRNAs between stages I-II PTC tumors samples and normal thyroid samples, and 122 aberrant expression lncRNAs between stages III-IV PTC tumors samples and normal thyroid samples. For patients with lymphatic metastasis, 160 aberrant expression lncRNAs were identified between stages I-II PTC tumors samples and normal thyroid samples, and 177 aberrant expression lncRNAs were identified between stages III-IV PTC tumors samples and normal thyroid samples ( Figure 3A). The intersected 83 specific intersected lncRNAs were selected initially for further analysis on ceRNA network to improve the credibility of bioinformatics analysis results (Table 3).

| Functional enrichment analysis
DAVID bioinformatics resources were employed to predict the biological functions and pathways of aberrantly expressed intersection mRNAs. A total of 18 633 mRNAs was identified from TCGA database, and 3531 mRNAs were aberrantly expressed in PTC tumor samples compared with normal thyroid samples (FC >2, P < 0.01). For patients without lymphatic metastasis, 2209 aberrant expression mRNAs were identified between stages I-II PTC tumors samples and normal thyroid samples, and 2457 aberrant expression mRNAs were identified between stages III-IV PTC tumors samples and normal thyroid samples. For patients with lymphatic metastasis, 2903 aberrant expression mRNAs were identified between stages I-II PTC tumors samples and normal thyroid samples, and 3128 aberrant expression mRNAs were identified between stages III-IV PTC tumors samples and normal thyroid samples. Finally, we identified 1822 mRNAs from the intersection of the four groups for further functional analysis. The upregulation and downregulation of these intersected mRNAs were further analyzed, respectively. GO analysis revealed 531 biological processes corresponding to these upregulated genes and 692 biological processes corresponding to these downregulated genes. Cell adhesion (GO:0007155) was the most enriched function for both upregulation and downregulation transcripts (Figure 4). KEGG pathway indicated 100 pathways matching with upregulated genes, and the top pathway was neuroactive ligand-receptor interaction. Furthermore, pathway analysis showed 101 pathways matching with downregulated genes and the top pathway was pathways in cancer ( Figure 5). Among these pathways, "MAPK signaling pathway" has already been demonstrated as the main cause of genetic alterations and PTC development. Targeting "MAPK signaling pathway" has shown an effective antitumor effect in preclinical studies and ongoing clinical trials. 30

| Prediction of miRNAs targets and construction of ceRNA network
A total of 1030 miRNAs were obtained from TCGA database to identify potential differentially expressed miR-NAs. And it was found that 87 miRNAs were aberrantly expressed in PTC tumor samples compared with normal thyroid samples (FC > 2, P < 0.05). According to the comparisons of four groups described above, 36 miRNAs were selected from the intersections of 87 PTC-associated miR-NAs to build the ceRNA network. Consequently, miRNAs target prediction including lncRNAs and mRNAs were conducted following the instructions described above. Firstly, 83 intersected lncRNAs and 36 specific miRNAs were collected to analyze the potential relationship in PTC. The potential miRNA-lncRNA interactions were predicted using StarBase v2.0 which could explore the candidate MREs. Finally, according to the analysis results, there were 29 specific miRNAs interacting with 54 specific lncRNAs (Table 4).
In the further analysis, the miRNA-targeted mRNAs were predicted using TargetScan and miRTarBase. It was found that 30 miRNAs interacted with 130 mRNAs  (Table 5). Then, we built a ceRNA network based on the above results (Tables 4 and 5). Cytoscape 3.0 was performed to depict the ceRNA network ( Figure 6). The ceRNA network contains 45 lncRNAs, 21 miRNAs, and 78 mRNAs finally.
The univariate Cox proportional hazards regression model was performed to detect the prognostic value of 45 specific ln-cRNAs. We found that three specific lncRNAs closely related with PTC patient overall survival (OS, log-rank P < 0.05). There are two lncRNAs (GGT3P and KRTAP5-AS1) positively relating to OS, and only one lncRNA (DPY19L2P4) negatively relating to OS (Figure 7).

| qRT-PCR validation
Evaluation of randomly selected eight lncRNAs (LOC100130238, HAND2-AS1, MIR9-3HG, LOC143666, EGFEM1P, LINC00284, TINCR, and ABCC6P1) was performed to prove the reliability of the above bioinformatics analysis (Table 7). Firstly, the expression status of the eight specific lncRNAs was detected in 28 newly diagnosed patients with PTC using qRT-PCR (Data S2). Compared with paired normal thyroid tissues, four lncR-NAs (LOC100130238, HAND2-AS1, MIR9-3HG, and LOC143666) were downregulated in PTC tumor tissues, whereas the remaining four lncRNAs (EGFEM1P, LINC00284, TINCR, and ABCC6P1) were upregulated in PTC tumor tissues ( Figure 8). All the eight lncRNAs were aberrantly expressed with the same trend and reached statistical significance (P < 0.05). Thus, the results were coherent with the above bioinformatics analysis. Subsequently, we performed the analysis of these eight lncRNAs and clinic features, and we found that HAND2-AS1 was significantly related to tumor size (Figure 9, P < 0.05). The result was identical with the analysis in TCGA database. Thus, our bioinformatics analysis is convincing based on the evaluation results.

| DISCUSSION
Thyroid cancer is a common malignant endocrine disorder worldwide with increasing incidence in last few decades. According to the National Cancer Center, PTC is the most frequently diagnosed cancer among Chinese women before 30 years old. 31 We cannot overlook the rapid rising in PTC, although the use of ultrasonography and fine needle aspiration in clinical may lead to overdiagnosis and overtreatment of PTC. In general, PTC has a relatively favorable prognosis; however, there are still about 5%-10% of patients suffering from an aggressive form of PTC. 32 Under such circumstances, thyroid cancer medicine is facing a huge challenge nowadays since concerns may cause overtreatment of PTC. Conventionally, risk evaluation for PTC is only based on clinicopathologic factors, such as age and lymphatic metastasis, which are definitely insufficient to accurately distinguish high-risk patients from low-risk patients. Thus, it is urgent to explore the underlying genetic background of PTC in order to identify the aggressive PTC and avoid overtreatment clinically.
Research concerning the genetic regulatory mechanism in thyroid cancer has shown sustained development. It is evident that aberrant expression of lncRNAs is closely related to oncogenesis in various cancer types. 33  is notable that lncRNAs are crucial in the genetic regulation and are responsible for cellular hemostasis. 34,35 Moreover, many studies showed that lncRNAs play an important part in tumorigenesis and are regulated by tumor suppressors or oncogenes transcriptionally. However, there are only a few studies focusing on lncRNAs in PTC using whole genome T A B L E 5 miRNAs targeting PTC-specific mRNAs miRNAs mRNAs hsa-miR-1179  HDAC4, NCAM1, NRCAM, RUNX1, RUNX1T1, SGCD   hsa-miR-138-1-3p  CCND1, CCND2, EFNB3, ENPP1, FUT9, GPR83   hsa-miR-138-5p  CTSH, EFNB3, ERBB4, EZH2, PLXNA4, PPARGC1A, RELN, SHANK2,  SLC16A2, SLC17A7, SOX4, TCF7L1, TRPC5, UNC5A   hsa-miR-144-3p  ERBB4, FGF7, FOSB, GABRB2, IRS1, ITPR1, LRP2,   and transcriptome sequencing technologies, 18,36,37 which demonstrated that dysregulated lncRNAs are essential in PTC carcinogenesis and may serve as diagnostic, prognostic biomarkers, and even therapeutic targets. Genomic analysis of PTC was firstly conducted using microarray in 2015 by Lan Xiabin et al, 21 which identified thousands of differentially expressed lncRNAs and mRNAs in PTC tissues compared with normal thyroid cancer. Additionally, genomewide expression screening in 12 paired PTC tumor tissues and normal thyroid tissues identified 218 differentially expressed lncRNAs in PTC, and two lncRNAs (XLOC-051122 and XLOC-006074) were found closely related to lymphatic metastasis. 22 Du et al 38 have performed a genomewide analysis in 18 PTC patients' tissues and 4 healthy donors' tissues, and build the coexpression network of lncRNA-mRNA. By now, much is unknown about lncRNAs and mRNAs, or ln-cRNAs and miRNAs in PTC based on ceRNA hypothesis. The ceRNA network in cancer research has drawn increasing attention as it was proposed in 2011. 19 It presented a hypothesis about how ceRNAs regulated other RNA transcripts by competing for miRNAs. It could be seen that ceRNA network has been built in many cancer types, such as colorectal cancer 39 and lung cancer. 40 These results indicate that lncRNAs prove to be significant in tumorigenesis. TCGA is the largest cancer genetic information database to study the genomic profiles of PTC. The ceRNA (lncRNA-miRNA-mRNA) network of PTC was once constructed based on the 348 PTC samples in TCGA. 24 In this study, dysregulated RNAs were identified in three PTC variants, including classical PTC, follicular PTC, and tall-cell PTC. However, the study missed part of the data in TCGA and lacked laboratory-based studies to validate the bioinformatics results. We aimed to identify the specific ln-cRNAs with diagnostic and prognostic roles in the study. Firstly, we retrieved the lncRNAs, miRNAs, and mRNAs profiles correlated with PTC stage and lymphatic metastasis in 461 PTC cases from TCGA database. Subsequently, the ceRNA network was built to predict the genetic interactions among the specific lncRNAs, miRNAs, and mRNAs. Next, we selected the intersected lncRNAs, miRNAs, and mRNAs in ceRNA network to assess their values in clinical relevance and prognosis based on the RNA sequencing data of 461 PTC samples and 55 normal samples. To validate the analysis results, the expression level of eight randomly selected lncRNAs (LOC100130238, HAND2-AS1, MIR9-3HG, LOC143666, EGFEM1P, LINC00284, TINCR, and ABCC6P1) was detected in 28 newly diagnosed tumor tissues and paired normal tissues in PTC patients using qRT-PCR. The evaluation displayed the consistent trend of up-and downregulation of selected lncRNAs with the expression level in TCGA. Eventually, the clinical information of donor samples was collected and compared with the data of TCGA. Similarly, the clinical relevance of lncRNAs was consistent with that of TCGA.

F I G U R E 6
The lncRNA-miRNA-mRNA ceRNA network. Red diamonds represent upregulated miRNAs, red balls, upregulated mRNAs, and red cones, upregulated lncRNAs. Green diamonds represent downregulated miRNAs, green balls, downregulated mRNAs; green cones, downregulated lncRNAs We classified PTC patients in TCGA into four groups based on TNM staging and lymph node status. Then, we compared the RNA sequence data of tumor tissues with that of normal thyroid tissues. The intersected mRNAs were finally decided accordingly. GO enrichment analysis and KEGG pathway analysis of intersected mRNAs were performed to reveal the potential mRNA functions. The top 25 GO enrichment terms indicated metastasis-associated functions, such as cell adhesion, proteolysis, and extracellular matrix organization. As for KEGG pathway analysis, top 25 pathways of intersected mRNAs included many cancer-specific pathways, such as MAPK, PI3K/Akt, and p53 signaling pathway. The results were consistent with the integrated genomic analysis in TCGA which demonstrated in detail about the importance of somatic genetic mutations, such as BRAF and RAS mutation, in the MAPK and PI3K pathways in PTC. 41 Many studies have shown that BRAF gene mutations increased the MAPK signaling pathway activation which could promote the initiation and development of PTC. [42][43][44] Generally, BRAF gene mutations are linked to larger tumor size, multifocality, extrathyroid extension, and lymphatic metastasis and predict poor prognosis. 7,45,46 It was reported that AZD6244, a MAPK pathway inhibitor, was strongly promising to treat thyroid cancer. 47 As the second most common gene mutation in thyroid cancer, RAS  48,49 Moreover, many drugs, such as temsirolimus and everolimus, targeting PI3K/Akt have been investigated in phases I to III clinical trials and proven to be effective in treating thyroid cancer. 50 In addition, abnormal p53 activation is closely related to the development of thyroid cancer. 51,52 In conclusion, the functional analysis was closely related to the development of PTC, which was also coherent with the classification standards. A great number of studies have presented that lncRNAs may function as ceRNAs and play an essential role in regulating gene expression. [53][54][55][56] For instance, lncRNA H19 and HULC may act through a ceRNAs manner by regulating let-7a/let-7b and miR-372/miR-373 which play an essential part in cholangiocarcinoma development. 57 We constructed ceRNA network of PTC with aberrantly expressed lncRNAs, miRNAs, and mRNAs in TCGA. Moreover, it was also found some cancer-specific lncRNAs, such as LINC01140, LINC00261, ABHD11-AS1, and TINCR, 58-61 served as important biomarkers in other different cancer types as well. We also found certain lncRNAs in ceRNAs played important roles in carcinogeneses, such as FOXD2-AS, DLEU2, FER1L4, HAND2-AS1, and TPTE2P1. [62][63][64][65][66] Next, we analyzed the biologic functions and pathways of identified mRNAs in ceRNA network. The results revealed many important functions and pathways, including MAPK, PI3K/Akt, and p53 signaling pathway.
We analyzed the correlation of 45 PTC-specific lncRNAs and clinic features, including age, race, gender, lymph node metastasis, TNM staging system, and patient outcome. The result showed that 37 specific lncRNAs were correlated with the above clinic characteristics. Eleven of these 37 lncRNAs have been reported to be closely related to various cancer types. For instance, lncRNA MIR31HG might function as an endogenous "sponge" for competing miR-193b binding to regulate gene expression and promote tumor progression in pancreatic ductal adenocarcinoma. 67 So far, there are few reports on the clinical relevance of these specific lncRNAs and PTC. More laboratory-based work needs to be done to verify the interaction of lncRNAs, miRNAs, and mRNAs in PTC. Eventually, we also analyzed the association between these 45 lncRNAs with patient OS. It was found that six lncRNAs (LOC100130238, PLEKHA8P1, GGT3P, MIR9-3HG, LINC00284, and EGFEM1P) were linked to PTC OS.
Subsequently, qRT-PCR validation of 8 lncRNAs in 28 paired samples was performed to assess the credibility and accuracy of the bioinformatics results. The results of qRT-PCR were consistent with the expression data in TCGA. Next, we assessed the clinical relevance of these eight lncRNAs based on our collected clinical information. The results suggested that HAND2-AS1 was closely related to tumor size, which is in coherence with the results of TCGA analysis. Hence, the bioinformatics analysis results of our study are credible and convincing.
Our study intended to identify the specific lncRNAs in PTC and investigate their genetic regulation by constructing the ceRNA network. Moreover, the expression status and clinical relevance of identified lncRNAs from TCGA database were assessed and compared with the data of donor samples. Finally, the credibility of bioinformatics analysis results was demonstrated and confirmed by qRT-PCR. Taken together, our results suggested that lncRNArelated ceRNA network might play an important role in the initiation and development of PTC. We also hope that our study can inspire researchers in this field to carry out further work.

ACKNOWLEDGMENT
The present study was supported by the technology development project of Nanjing (YKK10080). We appreciate Mr. Cheng donglin for his technical assistance.