The Effect of Botulinum Toxin Type A on Expression Profiling of Long Noncoding RNAs in Human Dermal Fibroblasts

Objective. This study was aimed at analyzing the expressions of long noncoding RNAs (lncRNAs) in Botulinum Toxin Type A (BoNTA) treated human dermal fibroblasts (HDFs) in vitro. Methods. We used RNA sequencing to characterize the lncRNAs and mRNAs transcriptome in the control and BoNTA treated group, in conjunction with application of GO (gene ontology) analysis and KEGG (kyoto encyclopedia of genes and genomes) analysis to delineate the alterations in gene expression. We also obtained quantitative real time polymerase chain reaction (qRT-PCR) to confirm some differentially expressed genes. Results. Numerous differentially expressed genes were observed by microarrays between the two groups. qRT-PCR confirmed the changes of six lncRNAs (RP11-517C16.2-001, FR271872, LOC283352, RP11-401E9.3, FGFR3P, and XXbac-BPG16N22.5) and nine mRNAs (NOS2, C13orf15, FOS, FCN2, SPINT1, PLAC8, BIRC5, NOS2, and COL19A1). Farther studies indicated that the downregulating effect of BoNTA on the expression of FGFR3P was time-related and the dosage of BoNTA at a range from 2.5 U/106 cells to 7.5 U/106 cells increased the expression of FGFR3P and COL19A1 in HDFs as well. Conclusion. The expression profiling of lncRNAs was visibly changed in BoNTA treated HDFs. Further studies should focus on several lncRNAs to investigate their functions in BoNTA treated HDFs and the underlying mechanisms.


Introduction
Botulinum Toxin Type A (BoNTA) is the most effective one among the seven neurotoxins secreted by Clostridium botulinum [1]. BoNTA causes muscle relaxation and this concept nowadays is widely being used in the cosmetic treatment of wrinkles [2,3]. In 2005 and 2008, some researchers found a face-lifting effect after intradermal injection of BoNTA to the mid and lower face [4,5]. However, some other researchers reported that needle pricks themselves without BoNTA can make the skin become smoother as well [6]. In order to find out whether BoNTA can affect the human dermal fibroblasts (HDFs) directly, in 2012, Oh et al. studied the in vitro effects of BoNTA on normal HDFs and found that BoNTA has a notable effect in increasing the level of collagen production and downregulating its degradation [7]. Collagen is the most abundant basic element of fibrous components in the dermis and is responsible for maintaining the structural integrity of the skin by joining cells together and to the extracellular matrix (ECM) [8,9]. These studies not only showed the positive effects of BoNTA on HDFs for remodeling skin but also implied the importance of HDFs. In 2016, Zhu et al. proved that topical BoNTA application could enhance the rejuvenation effect of fractional CO 2 laser, further indicating that BoNTA can refine skin texture via improving the activity of HDFs [10]. But until now, the molecular mechanisms through which BoNTA could affect HDFs are still not completely understood.
Long noncoding RNAs (lncRNAs) are a group of noncoding RNA transcripts longer than 200 nucleotides which cannot encode proteins [11]. In comparison with proteincoding genes, lncRNAs have limited coding potential and show little evolutionary conservation in sequence. Furthermore, some researchers have detected that lncRNAs expression is more tissue specific and at apparently lower levels [12]. LncRNAs, which were previously thought to be transcriptional "noise," are now proved to have some functions by regulating gene expression at the epigenetic, transcriptional, and posttranscriptional levels and participating in some biologic functions, such as genomic imprinting, chromosome modification, intranuclear transport, transcriptional activation, and interference [13]. Therefore, the understanding of cellular processes in physiological conditions will not be complete without analyzing the contributions made by lncRNAs. Until now, no information is available regarding the effect of BoNTA on expression profiling of lncRNAs in HDFs.
In this study, we investigated on lncRNA expression signature together with messenger RNA (mRNA) expression profile in BoNTA treated HDFs and confirmed the changing of some differentially expressed lncRNAs and mRNA using qRT-PCR. In conjunction, we also conducted functional analysis using Gene Ontology (GO) analysis and pathway analysis, in which genes are mapped to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Cell Separation and Culture.
Normal human skin samples were obtained from the prepuce of young healthy individuals in accordance with the ethics committee approval process of The First Affiliated Hospital of Nanjing Medical University (Nanjing, China). The acquirement of HDFs can be divided into two procedures. Initially dispase enzyme was used to separate the dermis and epidermis, and then collagenase enzyme was used to extract the HDFs. HDFs were grown in Dulbecco's modified Eagle medium (DMEM) with 1% penicillin-streptomycin and 10% fetal bovine serum in an environment of 5% CO 2 at 37 ∘ C. The cells used in our study were from passages 8-11.

Group Divisions and Botulinum Toxin Type A (BoNTA)
Treatment. In order to study differentially expressed lncR-NAs and mRNAs, we separated the cells into two groups, control group and BoNTA group: (1) control group: HDFs were grown in DMEM with 1% penicillin-streptomycin and 10% fetal bovine serum for 5 days and then serum-starved for 4 days, without receiving BoNTA treatment; (2) BoNTA group (48 h): HDFs were grown in DMEM with 1% penicillinstreptomycin and 10% fetal bovine serum for 5 days, serumstarved for 2 days, and then were grown in serum-free DMEM with BoNTA at a dose of 5 U/10 6 cells for 2 days.
In order to determine whether the changes of RNAs expression in BoNTA treated HDFs were time or dosage dependent, the cells were divided into 4 groups: (1) BoNTA group (24 h): HDFs were grown in DMEM with 1% penicillinstreptomycin and 10% fetal bovine serum for 5 days, serumstarved for 2 days, and then were grown in serum-free DMEM with BoNTA at a dose of 5 U/10 6 cells for 24 h; (2) BoNTA group (72 h): HDFs were grown in DMEM with 1% penicillin-streptomycin and 10% fetal bovine serum for 5 days, serum-starved for 2 days, and then were grown in serum-free DMEM with BoNTA at a dose of 5 U/10 6 cells for 72 days; (3) BoNTA group (48 h 2.5 U): HDFs were grown in DMEM with 1% penicillin-streptomycin and 10% fetal bovine serum for 5 days, serum-starved for 48 h, and then were grown in serum-free DMEM with BoNTA at a dose of 2.5 U/10 6 cells for 48 h; (4) BoNTA group (48 h 7.5 U): HDFs were grown in DMEM with 1% penicillin-streptomycin and 10% fetal bovine serum for 5 days, serum-starved for 2 days, and then were grown in serum-free DMEM with BoNTA at a dose of 7.5 U/10 6 cells for 48 h.
All groups were rinsed with Phosphate Buffer Solution (PBS) and the medium was changed every day, except during BoNTA treatment. Each method of detection consists of these six groups of cell culture, containing about 1 × 10 7 cells in each group, and there were at least 3 samples in each group. The BoNTA used in this study was manufactured by Lanzhou Institute of Biological Products Co., Ltd., Lanzhou, China.

Isolation of RNA and Preparation of Array Hybridization.
The RNA extraction was conducted using the TRIZOL reagent and then was dissolved in RNase-free water. The purified labeled genomic DNA was used for quantification and the RNA quantity was decided spectrophotometrically as A260/A280 ratio (1.9-2.1 were obtained). The collected RNAs were stored at −70 ∘ C for microarray analysis and qRT-PCR.

Microarray Analysis of lncRNAs and mRNAs Expression.
The Agilent Human lncRNA (8 * 60 K) arrays were designed in this experiment for measuring the expression profiles of lncRNAs and mRNAs. The lncRNA sequences were acquired from the following databases: NCBI-RefSeq, NONCODE v4, Ensembl, broad lincRNA, and frnadb v3.4. We conducted the sample labeling, microarray hybridization, and washing according to the manufacturer's proposals. Briefly speaking, we initially transcribed the total RNA to double stranded cDNA, synthesized them into cRNA, and labeled them with Cyanine-3-CTP and then hybridized the labeled cRNAs onto the microarray. We applied the Agilent Scanner G2505C (Agilent Technologies) to scan the arrays after washing.

Quantitative Reverse-Transcription Polymerase Chain
Reaction (qRT-PCR). In order to confirm the results obtained from microarray, we applied qRT-PCR to remeasure the abundance of differentially expressed lncRNAs and mRNAs selected from microarray analysis. RNAs in the medium were determined based on the protocol of KeyGen Biotech Co., Ltd., Nanjing, Jiangsu, China. In brief, we applied TRIzol to extract the total RNAs from HDFs and then synthesized the cDNAs from the separated RNA using SuperScript III Reverse Transcriptase (KeyGen, China). QRT-PCR was performed on ABI Prism 7700 Sequence Detector (Applied Biosystems). The reactions were performed at 9 ∘ C for 10 min, then at 40 cycles at 95 ∘ C for 15 s, and followed by at 60 ∘ C for 30 s. The 2 (−ΔΔCt) methods were used to evaluate relative quantification of lncRNAs and mRNAs expression. Primer sequences of lncRNAs and mRNAs for qRT-PCR are listed in Table 1. The results were expressed as mean ± SD (standard deviation) of each independent experiment.

Data Analysis.
We applied the Feature Extraction software (version 10.7.1.1, Agilent Technologies) to analyze array images to get raw data and used GeneSpring (version 13.1, Agilent Technologies) to complete the fundamental analysis with the raw data. First of all, the raw data was normalized with the quantile algorithm. The probes which have signed with P ( value ≤ 0.05 is recommended) were selected for further data analysis. Deferentially expressed lncRNAs and protein-coding RNAs were then identified through fold change. The critical value set for up-and downregulated genes was a fold change ≥ 2.0. Afterwards, GO analysis and KEGG analysis were applied to depict alterations in the gene expression. SPSS 17.0 software (SPSS Inc., Chicago, IL, United States) was used to analyze the statistical data. The differences in RNA expression between the two groups and more than three groups were analyzed using the Student's t-test and oneway ANOVA, separately. < 0.05 was considered to indicate a statistically significant difference.  Table 2 (fold change > 5) and the distinctly expressed mRNAs are summarized in Table 3 (fold change > 3).

GO and Pathway Analysis of Deferentially Expressed
RNAs. The upregulated genes were involved in 898 biological processes, 171 cellular components, and 236 molecular functions separately. The response to cell-cell signaling (−log 10 (p value) = 3.7114) was the most significant term among the biological process category. The extracellular region (−log 10 (p value) = 2.8128) was the most represented GO term in   the cellular component category. Heparin binding (−log 10 (p value) = 3.3714) was the most highly represented term within the molecular component category. The downregulated genes were involved in 880 biological processes, 179 cellular components, and 240 molecular functions, respectively. The most significant term was the response to mitotic cell cycle (−log 10 (p value) = 31.3859) among the biological process category. The most represented GO term was the condensed chromosome kinetochore (−log 10 (p value) = 13.0822) in the cellular component category. Heparin binding (−log 10 (p value) = 6.2714) was the most highly represented term with in the molecular component category. The top ten of the number of deferentially expressed genes and the significance analyzed by GO term are demonstrated in Figure 4. KEGG (Kyoto Encyclopedia of Genes and Genomes) was used to analyze the pathway enrichment. The pathway analysis showed that changed genes participated in the DNA replication, cell cycle, P53 signaling pathway, pathway in cancer, and melanoma. The upregulated genes participated in 157 pathways while downregulated genes participated in 169 pathways. The top ten pathways among them are shown in Figure 5.

The Regulation of BoNTA on the Expression of FGFR3P
and COL19A1. We also performed qRT-PCR to validate the altered lncRNA expression at different time-points and different dosage of BoNTA treatment. We found that the level of FGFR3P reached the lowest at 48 h but indicated an upward trend at 72 h when treated by BoNTA at the dosage of 5 U/10 6 cells ( Figure 6). The level of FGFR3P was gradually increasing with the increasing dosage of BoNTA (at a range from 2.5 U/10 6 cells to 7.5 U/10 6 cells). At the same time, the level of COL19A1 in HDFs was gradually increasing with the extension of time treated by BoNTA at the dosage of 5 U/10 6 cells (Figure 6), and its level showed a trend of increasing when the dosage of BoNTA ranges from 2.5 U/10 6 cells to 7.5 U/10 6 cells ( Figure 6).

Discussion
Recently, lncRNAs have come into the extent of probing into human biological functions and diseases [13][14][15][16][17]. Until now, only a few lncRNAs had been reported about HDFs, and studies on the expression of lncRNAs in BoNTA treated HDFs were still lacking. In this study, we identified some differentially expressed lncRNAs and mRNAs between the BoNTA treated group and control group by analyzing gene expression profiles. The number of changed lncRNAs was greater than that of mRNAs. Although currently the precise roles of the changed lncRNAs in BoNTA treated HDFs are not clear, lncRNAs have been regarded as vital regulators of gene expression and have various biofunctions. There are a large number of evidences demonstrating that lncRNAs can regulate gene expression by forming RNAprotein, RNA-RNA, DNA-RNA, and protein-DNA interactions [18]. Although no direct relationship was found between the altered lncRNA and mRNA expressions, we are convinced that the clear changes of lncRNAs (RP11- 517C16.2-001, FR271872, LOC283352, RP11-401E9.3, FGFR3P, and XXbac-BPG16N22.5) in HDFs are in response to BoNTA and are related to the changes of protein-coding RNAs. We suspected that the decrease of NOS2 indicated the regulation of cell proliferation process [19,20]. The increase of FOS induced by BoNTA showed the regulation of proliferation, which is also involved in the cellular senescence process of HDFs [21][22][23][24]. BIRC5 has been reported to participate in modulation of diverse cellular processes such as proliferation, adhesion, apoptosis, migration and invasion during growth, development, repair, maintenance, and regression of a wide variety of mesenchymal tissues [25][26][27][28]. The downregulated BIRC5 induced by BoNTA indicated the decrease of apoptosis in HDFs. PLAC8 has been proved to regulate cell cycle and participate in the regulation of apoptosis and cell division [29][30][31].
GO database was applied to analyze the function of the differentially expressed genes. The results were divided into three sections: biological process (BP), cellular component (CC), and molecular function (MF). Besides the top ten of the number of differentially expressed genes and the significance analyzed by GO term demonstrated in Figure 4, the analysis showed that the differentially expressed genes were also involved in a variety of other biological functions, such as negative regulation of autophagy (GO: 0010507), positive regulation of collagen biosynthetic process (GO: 0032967), regulation of G1/S transition of mitotic cell cycle (GO: 2000045), and positive regulation of fibroblast proliferation (GO: 0048147). The analysis indicated that the lncRNAs can affect the function of HDFs by regulating the expression profiles of genes related to HDFs. KEGG analysis can provide some suggestive information about the potential relativity of the changed gene expression with the alterations of the pathways. The result of our study indicates several significant pathways related to a variety of functions, such as cell proliferation, cell cycle, apoptosis, and DNA replication. It had been reported that BoNTA can regulate the process of cell cycle and DNA replication by other researchers. For example, G. Karsenty et al. reported in their study that BoNTA obviously reduced LNCaP cell proliferation and increased apoptosis in a dose-dependent manner [32]. Park et al. found that BoNTA upregulates the expression of cell cycle related genes such as RhoA, Rac1, and Cdc42 in a dose-dependent manner [33]. Our prediction results accord with the functional analysis of BoNTA obtained from other investigations.
QRT-PCR was applied to verify the result of microarray analysis. We learned from other studies that BoNTA can regulate cell proliferation and collagen synthesis; the mechanisms may play important roles in skin rejuvenation effects of BoNTA [32][33][34]. The RNAs chosen for qRT-PCR confirmation meet at least one of two following criteria: 1: the RNAs being at least 2-fold differently expressed in BoNTA treated groups in comparison with the control group according to the data analysis; 2: the RNAs that had been proved in other studies to have a relationship with biological functions such as cell proliferation and collagen synthesis.     Figure 4: Bioinformatic analysis of the differentially expressed genes. The p value denotes the significance of GO terms enrichment in the differentially expressed genes. The lower the p value, the more significant the GO term (p value ≤ 0.05 is recommended). We can choose the target genes for further study based on the combination of the analysis provided by GO and the biologic significance. Significant pathways of downregulated genes −log 10 (p value) −log 10 (p value) Figure 5: Pathway analysis of the differentially expressed genes. Fisher's exact test was used to select the main pathway, and the significance threshold was defined with p value. The lower the p value, the more significant the pathway (p value ≤ 0.05 is recommended). We can get some information about the possibility between the differentially expressed genes and the change of cellular pathways.
consistent with the microarray data which confirmed the reliability of our microarray analysis. Next, we applied qRT-PCR to further investigate the expression changes of COL19A1 and FGFR3P after BoNTA treatment at different dosages and culture times. COL19A1, as one member of the fibril-associated collagens with interrupted triple helices (FACIT) group, is thought to act as a cross-bridge between extracellular matrix molecules (ECM) and is involved in the formation of the well-known striated fibrils [35]. Bioinformatic analysis reveals that FGFR3P which is located on chromosome 6 is defined as the pseudogene of Fibroblast Growth Factor Receptor 3 (FGFR3), noncoding RNA. FGFR3 is one member of FGFR family and had been reported to participate in the regulation of cell proliferation, apoptosis, migration, and angiogenesis in many cells including HDFs [36][37][38][39]. Although pseudogenes have been considered as the remnants of functional genes which have no coding ability over a long period of time, evergrowing number of studies have proved that some pseudogenes have diverse functions, including serving as miRNA decoys, functioning as antisense transcripts, encoding short peptides, and producing siRNAs or proteins [40][41][42][43]. Our results showed that downregulating effect of BoNTA on the expression of FGFR3P was time-related. Meanwhile, the dosage of BoNTA at a range from 2.5 U/10 6 cells to 7.5 U/10 6 cells increased the expression of FGFR3P and COL19A1 in HDFs as well. We speculated significant change of expression of these RNAs in a time-dependent manner after BoNTA treatment which may be of some clinical significance. Zhu et al. conducted one clinical experiment and found that topical application of BoNTA could enhance the rejuvenation effect of fractional CO 2 laser; they also found that best results of skin detection were in the last observation point at 3 months after treatment [34]. Zhu et al. also performed another clinical study and found that intradermal BoNTA injection showed its best rejuvenation results in the last observation point at 12 weeks after treatment as well [44]. The results of the previous clinical investigations showed that BoNTA had a persistent promoting effect of collagen synthesis with the extending of time. Taking these into consideration, we speculate that BoNTA have the antiaging effect not only by relaxing in muscle fiber, but also by promoting activity of HDFs directly. The underlying mechanisms still need further studies.
In conclusion, for the first time, our current study identified the changes of expression profiles of lncRNAs in BoNTA treated HDFs and found that BoNTA dynamically regulated the expression of COL19A1 and FGFR3P in HDFs, indicating the potential role of several lncRNAs in BoNTA treated HDFs. Therefore, further explorations are warranted to discover the mechanisms behind the dynamic changes BoNTA induced on COL19A1 and FGFR3P and the findings obtained in this study should lay foundations for further studies into the potential roles of these altered lncRNAs in BoNTA treated HDFs.