Lncrna Rn7sk and Cacna1g-As1 Related to the Scar Fibroblasts Proliferation in Burn Scars and Human Skin Cell

Lihua Wang and  Zhiming Xin

Published Date: 2019-08-26
Lihua Wang1 and Zhiming Xin2*

1Department of Dermatology, Jinan Central Hospital Affiliated to Shandong University, Jinan, 250013, China

2Department of Burn and Plastic Surgery, Shengli Oilfield Central Hospital, Dongying, 257034, China

*Corresponding Author:
Zhiming Xin
Department of Burn and Plastic Surgery
Shengli Oilfield Central Hospital, No. 31 Jinan
Road, Dongying City, Shandong Prov
P. R. China
E-mail: xinzhiming612@163.com

Received Date: June 07, 2019 Accepted Date: August 16, 2019 Published Date: August 26, 2019

Citation: Wang L, Xin Z (2019) Lncrna RN7SK and CACNA1G-AS1 Related to the Scar Fibroblasts Proliferation in Burn Scars and Human Skin Cell. Cell Mol Med Vol. 5 No. 1:1.

Visit for more related articles at Cellular & Molecular Medicine: Open access

Abstract

Background: Several mechanisms are thought to be essential for the development of burn scars, but there is still a challenge to study the lncRNA expression response to burn injury in the skin. The purpose of this study was to examine the changes in lncRNA expression in human skin after burn injury.

Methods: This was done by comparing pre-injury tissue from otherwise healthy adults, by virtue of a microarray gene expression profile. We first ranked gene sets (or gene signatures) that identify each class. A multiple-class classifier was built based on the gene sets. Then, gene networks associated with scars and burn injury was built in order to identify the hub genes to each class. Finally, the hub genes in the network were identified as important genes and tested in an In-vitro cell model.

Results: The set of genes considered significant for each of the classes is determined by a common threshold for the posterior probability. Burn wound has been assigned 6 genes as signatures, pre-injury healthy skin 5 genes, while scar only 2 genes. By gene networks, lncRNA RN7SK and CACNA1G-AS1 were shown with the higher discriminant power. In fibroblasts, overexpression of RN7SK and CACNA1G-AS1 increased migration and wound healing-rate by approximately 20% and 10% compared with the control group.

Keywords

Burning injury; Human dermal fibroblast; lncRNA expression; Classification model; In-vitro cell models

Introduction

The development of hypertrophic scars is a common occurrence and the most common complication after burn injury [1]. Progress has been made in developing gene expression patterns that share some molecular characteristics with human scars [2,3]. However, these studies do not focus on non-coding RNAs. Additionally, the cost and expertise required to manage these large scales of gene expression studies may preclude many investigators from incorporating these models into their experiments [4].

Following a current consensus of burn researchers and clinicians, several research priorities were proposed. Particularly, it is pointed out that gene expression plays an important role in early and serial biopsies of human skin [5,6]. Further, recommendations were also made that future studies in burn research were needed to include high throughput gene expression analysis by cutting-edge bioinformatics pipelines [7,8]. Long non-coding RNA, which has similar sequence characteristics with proteincoding genes, such as being usually the 50 cap and 30 poly (A) tail, spliced and most transcribed by RNA polymerase II [9], is defined as a kind of RNA molecules with length longer than 200 nucleotides [10]. Although lacking the protein-coding function, lncRNA, acknowledgedly, plays a significant role in regulation of transcription and translation of the coding RNA in the scar formation [11,12], through regulating a series of biological processes across genomic imprinting [13], maintenance of genome integrity, cell cycle control [14,15], development and differentiation [16].

The purpose of this study was to identify lncRNAs in pre-injury, early and serial samples of untreated human burn wounds. In order to identify lncRNA gene signatures of burn injury and scars [17], we set up an analysis based on the geNetClassier algorithm [18], which is designed to build transparent classifiers and the associated gene networks based on genome-wide expression data. This algorithm shows a robust performance applied to patient-based gene expression datasets that study disease subtypes or disease classes.

Here in our research, we first ranked gene sets (or gene signatures), which integrates several existing machine learning and statistical methods, identifies burn injury and pre-injury sample. The feature ranking was achieved based on a Parametric Empirical Bayes method (PEB). Double-nested internal crossvalidation (CV) [19,20] was used for the feature selection process and to estimate the generalization error of the classifier. The machine learning method implemented in the classifier was a multi-class Support Vector Machine (SVM) [21]. Then, the gene networks associated with each class were identified [22,23]. Finally, the hub genes in the network were identified as key gene signatures and tested in an In-vitro cell model.

Methods

Data collection

Microarray data with the accession number E-MTAB-1323 was chosen and obtained for this study. The microarray dataset E-MTAB-1323 deposited by Vincent Gabriel was based on the platform of the Illumina HumanWG-6 v3.0 Expression BeadChip. In our study, we show its performance for a dataset that includes 30 microarray samples from 5 healthy skin subjects, 15 burn wound subjects, 10 scar controls. This dataset contains an expression profile with mRNA extracted from burn scar biopsies of patients of the 3 major types (burn wound, pre-injury, and scar).

Data pre-processing

Gene expression profiles were preprocessed using a chip description file. This chip description file, with gene-based annotation files for the Illumina expression microarrays, allows mapping the expression directly to genes (Ensembl IDs ENSG) instead of probe sets. To translate these Ensembl gene IDs into Gene Symbols for easier reading, the optional argument geneLabels from geNetClassier can be used. This option allows extending the annotation and labeling of the genes by providing a table that contains the gene symbol and other characteristics of the genes in the expression set. This option can be used with any annotation as long as it is provided in the correct format. The raw downloaded data in the original chip files were preprocessed using the Robust Multi-array Average (RMA) algorithm in Bioconductor (version 3.6) in R language (version 3.4.3).

The process could realize the background correction, quantile normalization and probe summarization of the raw microarray data. The gene expression value was log2 transformed. Finally, we got 48,803 gene expression profile data. A total of 503lncRNAs were selected for further analysis.

Classification

Filtering data and calculating the ranking of the genes: A first version of the ranking is built by ordering the genes decreasingly by their posterior probability for each class. To resolve the ties, our model uses the expression difference between the mean for each gene in the given class and the mean in the closest class. The first step of classification algorithm is to determine a ranking of genes for each class based on the analysis of the expression signal. To create this ranking, we use a Parametric Empirical Bayes method [19,20,24]. This method implements an Expectation- Maximization (EM) algorithm for gene expression mixture models, which compares the patterns of differential expression across multiple conditions and provides a posterior probability. The posterior probability is calculated for each gene-class pair, and represents how much each gene differentiates a class from the other classes; being 1 the best value, and 0 the worst. In this way, the posterior probability allows finding the genes that show significant differential expression when comparing the samples of one class versus all the other samples (One-versus-Rest comparison). In addition, the genes with a posterior probability greater or equal to 0.95 are filtered out before proceeding into further steps.

The classifier’s generalization error: The genes are taken in order from the genes ranking of each class until any of the classes reaches gets to the maximum number of genes or until zero error is reached. The error for each of the classifiers and the number of genes used to construct them are saved.

Construction of the classier: Selection of a subset of genes to train the classier through 8-fold cross-validation. The selected genes are used to train the classier with the complete set of samples. The fastest execution would be training the classier exploring a reduced number of genes (maximum Genes Training cycle is 100).

Construction of the gene networks: A gene network is built for each one of the classes using the pairwise gene-to-gene correlations and interactions. The gene networks are built calculating the relations derived from gene to gene co-expression analysis (by default, Pearson correlation) and the interactions derived from gene mutual information analysis (using minet package) [22]. In order to identify the importance of the genes within the network, we calculate the degrees of nodes (genes) in the network [18].

Validation with In-vitro cell models

Adult human dermal fibroblasts (aHDFs): Adult human dermal fibroblasts (aHDFs) were purchased from Lonza Group, Ltd. (Walkersville, MD, USA) and maintained in fibroblast basal medium-2 (FBM-2) supplemented with growth kit (10 ml of fetal bovine serum, 0.5 ml of insulin, 0.5 ml of gentamicin sulfate amphotericin-B (GA-1000) and 0.5 ml of r-human fibroblast growth factor-B, Lonza). Cells were incubated at 37℃ in a 5% CO2 atmosphere. To build the RN7SK and CACNA1GAS1 overexpression aHDFs cell models, the full-length RN7SK and CACNA1G-AS1 cDNA were transfected into aHDFs cells to generate aHDFs/RN7SK and aHDFs/CACNA1G-AS1, respectively.

Cell migration assay (wound-healing assay): A wound closure seeding model was constructed using silicon culture inserts (Ibidi, LLC, Munchen, Germany) with two individual wells for cell seeding. Each insert was placed in a culture dish, and 8 × 103 cells of aHDF were plated in each well and grown to form a confluent and homogeneous layer. Twenty-four hours after cell seeding, the culture insert was removed and a cell-free area, the “wound” made by the culture insert, could be observed. The wound was approximately 500 nm wide. The cells were treated with 10 g/ml of mitomycin C in serum-free media without growth factors for 2 h to suppress proliferation. Healing of the wound by migrating cells after lncRNA cDNA transfection was observed over time by light microscopy (IX-70, Olympus) and analyzed using Image J software.

Results

Genes ranking

As a result of this process, even if a gene is found associated with several classes during the expression analysis, each gene can only be on the ranking of one class. When the analysis performed for two classes, the one vs rest approach only provides one list of genes. This list of genes, the genes that differentiate the classes, is labeled both classes. The sign of the expression value reflects in which class it is up or down. The genes ranking obtained for each class is used for the gene selection in the classification procedure and it is provided as Figure 1.

cellular-molecular-medicine-overlap

Figure 1: (A) Scheme representing the overlap between the sets of genes that each disease may affect. All the genes that affect each disease (ovals) and selects as significant, the genes that are unique (differentially expressed) to each group (B).

Significant genes

The set of genes considered significant for each of the classes is determined by a common threshold for the posterior probability. This common threshold provides a way to quantify the size of the gene signature assigned to each group (as always: compared to the other diseases in the study). In this way, the algorithm provides a framework to compare biological states, i.e. the biological or pathological conditions represented in the samples. The results (Figure 2) showed the differences in the size of the gene sets assigned to a disease: at lpThreshold 0.95 burn wound has been assigned 6 genes, while scar only 2 genes.

cellular-molecular-medicine-probabilities

Figure 2: Plot of the posterior probabilities of the genes of burn wound class, scar class, and pre-injury healthy skin, ordering the genes according to their rank and setting the lpThreshold at 0.95.

This observation showed the biological interpretation that small gene signatures might be an indication of relatively less systemic changes. The results may help us to unravel that differences based on the gene signatures. It indicated that differences in lncRNA gene expression in the scar is relatively small, compared to the burn wound in this study.

Gene selection procedure

The optimum number of genes for training the classier is selected by evaluating the classifiers trained with an increasing number of genes. This is done using several iterations of 8-fold crossvalidation. Each cross-validation iteration starts with the first ranked gene of each class: it trains an internal classier with these genes and evaluates its performance.

The final selection is done based on the genes selected in each of the iterations. For each class, the top-ranked genes are selected by taking the highest number of genes selected in the crossvalidation iterations. This allows identifying a stable number of genes while accounting for the differences in sampling. The results were shown in Figure 3.

cellular-molecular-medicine-iteration

Figure 3: (A): the number of genes selected in each iteration. The algorithm runs until exploring a maximum number of genes in any class (maxGeneTrain=100) or until zero error is reached (continueZeroError=FALSE). In each iteration, the minimum number of genes with minimum error is selected. (B): Plot of the gene-selection iterations. Each line represents an iteration and the error rates observed for each number of genes (starting at 5, one per class).

Classification model built

We finally built the classier associated to each class, and the genes ranking and further information about the selected genes. This classier was built based on total 17 of support vectors with parameters as follows: SVM-Type: C-classification; SVM-Kernel: linear; cost: 1; gamma: 0.083. In the dataset, the genes were labeled with the gene symbols provided by a gene-based probe mapping file. The best lncRNA with best discriminant power was shown in Figure 4.

cellular-molecular-medicine-discriminant

Figure 4: Plot of the discriminant power of gene RN7SK (A) presents the highest discriminant power in burn wound while CACNA1G-AS1 (B) in scar. Plot of the expression profiles for gene RN7SK (C) and CACNA1G-AS1 (D). A high discriminant power can help to identify gene markers.

Gene networks

Gene networks provide the available information about the genes, which was shown in Figure 5. The gene name is the node label. The gene expression is shown with the node color, meanwhile, the discriminant power determines its size.

cellular-molecular-medicine-wound

Figure 5: Gene network obtained for class burn wound and scar selecting the top 100 genes from the gene ranking, presenting all nodes. The network legend indicates the meaning of the shapes and colors given to the nodes and edges.

Experimental validation of hub genes by analysis of the gene effect of the migration rate of skin cells

In the cell viability assay, Figure 6 shows the effects of overexpression of RN7SK and CACNA1G-AS1 about the migration of normal human skin cells. Compared with the control group, overexpression of aHDFs/RN7SK and aHDFs/CACNA1G-AS1 cells migrated faster, and the most effective concentration was varied along with cell type. In fibroblasts, overexpression of RN7SK and CACNA1G-AS1 increased migration rate most effectively and improved wound healing by approximately 20% and 10% compared with the control group.

cellular-molecular-medicine-migration

Figure 6: Plot of three groups of cell migration rate. aHDFs/RN7SK and aHDFs/CACNA1G-AS1 are both experimental groups, the aHDF represents control group. Y-axis is cell migration rate.

Discussion

In this study, we used a dataset including pre-injury, early and serial samples of untreated human burn wounds. We aimed to mine the hub lncRNA in burn or scar from aesthetic scarification individuals. The creation of an aesthetic scar, which adopted as a form of body modification, may be done by applying a heated object to the skin, frostbite, electrocautery or full-thickness skin excision [8,25]. Aesthetic Scarification has a history in many cultures and in some cases may even be used as a traditional medical treatment without supporting evidence [7,26]. This kind of scars is different from patients with severe burning scars which have a risk of developing sepsis or other complications.

Although RNA expression is thought to be significant in regulating dermal collagen production in scars, less has been described regarding the function of lncRNA both in scars and aHDFs [27]. Therefore, we selected the lncRNA genes for our candidates. Because the geNetClassier algorithm relies on a pre-specified number of clusters of genes [18], genes were enriched by the SVM classification based on the expression of lncRNAs [20]. In addition, two new significant genes are described here. First is the rapid transient response of the gene RN7SK in burn wound. Considering the important role that inflammation and the Positive Transcriptional Elongation Factor b (P-TEFb) [28] or signaling proteins are thought to have in scar development, this finding implies that there may be an opportunity to intervene early after an injury to investigate potential scar modification treatments rather than later after the scar has been established. Second is the sustained regulation of lncRNA including CACNA1GAS1. Using the algorithm of geNetClassier, some consistencies were observed. Bioinformatics analysis suggested that lncRNA CACNA1G-AS1 may be crucial to keloid formation [27]. This observation consists of the fact that differentially expressed lncRNAs, such as CACNA1G-AS1, may play a key role in keloid formation in aHDFs.

There are clear limitations in this study. First, the sample is small. Given the sample size, it is impossible to know what if any effect each of these factors may have the result. However, given the exploratory nature of this study as well as the available resources with poorly described population, expanding the study was not possible. Secondly, using a small volume biopsy, cell and beadarray technology did not allow for additional histologic descriptors of the wound and resultant scar. However, this methodology has been previously reported in skin and scar investigations [11,29].

Conclusion

Overall, the response to injury is consistent between subjects and cell models. This cohort provides a unique insight into the skin transcriptional changes that occur after burn injury. Other studies have described transcription differences of lncRNA RN7SK and CACNA1G-AS1 between normal skin and scar or validated in cell models. However, there was considerable variation in sampling time that needs to be further investigated.

References

open access journals, open access scientific research publisher, open access publisher
Select your language of interest to view the total content in your interested language

Viewing options

Flyer image

Share This Article