Abstract
Li Yang1,*, Lechun Lyu2,*, Wenjuan Wu1,*, Dongyun Lei1,*, Ying Tu1, Dan Xu1, Jiaqi Feng1 and Li He1
1Department of Dermatology, The First Affiliated Hospital of Kunming Medical University, Kunming, China
2Technology Transfer Center, Kunming Medical University, Department of Physiology, Kunming Medical University, Kunming, China
*These authors contributed equally to this work
Correspondence to:
Li He, email: [email protected]
Keywords: sensitive skin; lncRNA; mRNA; RNA sequencing
Received: August 23, 2017 Accepted: November 14, 2017 Published: December 12, 2017
ABSTRACT
Sensitive skin (SS) is a condition of subjective cutaneous hyper-reactivity. The role of long non-coding RNAs (lncRNAs) in subjects with SS is unclear. Therefore, the aim of the present study was to provide a comprehensive profile of the mRNAs and lncRNAs in subjects with SS. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis presented the characteristics of associated protein-coding genes. In addition, a co-expression network of lncRNA and mRNA was constructed to identify potential underlying regulation targets; the results were verified by quantitative real-time PCR (qRT-PCR) and RNA-seq analyses in patients with SS and normal samples. Compared with the normal skin group, 266 novel lncRNAs and 6750 annotated lncRNAs were identified in the SS group. A total of 71 lncRNA transcripts and 2615 mRNA transcripts were differentially expressed (P < 0.05). The heat signature of the SS samples could be distinguished from the normal skin samples, whereas the majority of the genes that were present in enriched pathways were those that participated in focal adhesion, PI3K-Akt signaling, and cancer-related pathways. Five transcripts were selected for qRT-PCR analysis and the results were consistent with RNA-seq. The results suggested that LNC_000265 may play a role in the epidermal barrier structure of patient with SS. The data suggest novel genes and pathways that may be involved in the pathogenesis of SS and highlight potential targets that could be used for individualized treatment applications.
INTRODUCTION
Sensitive skin (SS) is a broad term used to describe a multitude of clinical findings that are attributed to different sensory perceptions, namely facial irritation, burning, stinging, tightness, tingling, pain, and pruritus [1]. The sensitive skin syndrome (SSS) is considered as a state of hyperactivity to specific environmental stimuli that is caused from a single and/or a number of underlying pathologies [2]. The main disadvantage encountered during the diagnosis of the disease is the lack of an objective-screening test [3]. The complex nature of the skin disease syndrome requires the use of a diagnostic algorithm and the need to test patients with multiple patch testing, prior to the establishment of a definite diagnosis, as it has been shown from the lack of association between different allergens in subjects with positive allergic reactions with regard to each allergen alone (SDS and/or lactic acid) [3, 4]. Nevertheless, several studies have suggested a link between SS and disruption of the epidermal barrier function, resulting in the perception of skin discomfort [5, 6]. Despite these promising findings, the molecular network that contributes to the development of SS has not been elucidated to date.
Long non-coding RNAs (lncRNAs) are a class of RNA sequences that are more than 200 nt in length and are involved in the regulation of translation process, although they do not possess protein coding potential [7]. A multitude of studies have shown that lncRNAs are involved in the regulation of developmental processes and in the progression of several human diseases [8–11], while their expression and localization varies among different cell types and subcellular compartments [12–17]. LncRNAs have been found crucial to genomic imprinting, dosage compensation, and pluripotency-regulation [18, 19]. The rapid progress of RNA sequencing (RNA-seq) promoted the exploration and research of non-coding RNAs, and novel lncRNAs have been identified by different pipelines using RNA-seq data [20, 21]. RNA-seq exhibits several advantages compared to the previously established methodologies that have been used for the evaluation of the complete set of transcripts in the cell, such as hybridization-based approaches and specialized microarrays. RNA-Seq has been successfully used to provide a ‘digital measurement’ of the gene expression difference between different tissues [22]. Although no reports on the contribution of lncRNA in the development of SS have been reported, it has been suggested that these RNA sequences play a significant role in skin homeostasis and related skin diseases [23, 24]. Several studies demonstrated the involvement of lncRNAs in the differentiation and maintenance of normal human keratinocytes and epidermal tissues [25, 26]. For example, lncRNA such as anti-differentiation non-coding RNA (ANCR) and terminal differentiation-induced non-coding RNA (TINCR) are vital for epidermal stability [27, 28]. Sonkoly et al. identified a novel lncRNA, namely psoriasis susceptibility-related RNA gene induced by stress (PRINS) that is involved in the susceptibility to psoriasis [29]. The authors suggested that PRINS may play an important role in psoriasis by evidence derived from psoriasis patients and in vitro cell culture experiments [29]. The use of bioinformatics methods has been adopted in the investigation of the genes involved in the development of atopic dermatitis and psoriasis. Notably, differentially expressed genes (DEGs) were associated with epidermis development and immune response in atopic dermatitis [30]. Similarly, enrichment analysis of psoriasis- correlated modules revealed that pathways involved in short chain fatty acid metabolism, olfactory signaling, and regulation of leukocyte-mediated cytotoxicity were the main pathways in which the DEGs were identified [31]. Of note is that more than 50% of the co-expressed genes in 18 psoriasis patients and 16 healthy controls were lncRNAs [31].
In view of the significant roles of lncRNA in the regulation and the differentiation of epidermal homeostasis and the disruption of the epidermal barrier function in SS, we hypothesized that lncRNAs may also take part in the pathogenesis of SS. The aim of the present study was to provide a more comprehensive and validated conclusion regarding the identification of differentially expressed lncRNAs in SS tissues.
RESULTS
Overview of RNA-Seq and mRNAs and lncRNAs identification
The RefSeq (Build 37.3) and the GENCODE v19 databases were selected as the annotation references for mRNA and lncRNA analyses, respectively. Tissues from three patients with SS and three normal skin tissues were used to construct the RNA-seq library. The reads were mapped using Cufflinks and this resulted in 233,945 assembled transcripts. The number of transcripts corresponding to more than two exons were selected in order to filter abundant of low-expression, low-confidence single exon transcripts. In addition, the transcripts with a length of >200 bp were selected and the known functional genes were removed. The transcripts that contained an exon area overlapping with the annotation database were screened by the Cuffcompare function. The lncRNAs with overlapping exon areas assembled in the database were included as annotation lncRNA database into subsequent analysis. Furthermore, the expression of each transcript was calculated by Cuffquant, and the transcripts with expected number of fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM) ≥0.5 were selected. This selection process yielded 183,814 filtered transcript isoforms. Finally, the protein-coding transcripts were excluded by Phylogenetic Coding Substitution Frequency (PhyloCSF), PfamScan (v1.3), Coding Potential Calculator (CPC), and Coding-Non-Coding-Index (CNCI) [21]. The transcripts identified by the four methods were deemed as confidently expressed lncRNAs. The final screen included a coding potential score lower than 0 (CPC < 0), CNCI and PLEK scores lower than 0 (<0) and Pfam and E-values lower than 0.001 (<0.001). The final results yielded 887 novel assembled lncRNAs that corresponded to 266 novel lncRNA transcripts (Figure 1A, 1B), including 236 (88.7%) lncRNAs and 30 (11.3%) antisense lncRNAs (Figure 1C).
Figure 1: Transcriptome analysis of lncRNA by RNA-seq in three skin samples of subjects with SS and three normal skin samples. (A) Venn diagram of screening results. The sum of the numbers in each large circle represents the total number of non-coding transcripts of the software, and the overlapping parts of the circle represent the non-coding transcripts common to the software. Fragments per kilobase of transcript sequence per millions base pairs sequenced (FPKM), phylogenetic coding substitution frequency (PhyloCSF), coding potential calculator (CPC) and coding-non-coding-index (CNCI). (B) The abscissa represents the screening step and the ordinate represents the corresponding steps after screening the number of transcripts. (C) Pie chart distribution of novel lncRNAs identified based on antisense and intronic forms.
A total of 6750 annotated lncRNAs, 2718 (40.27%) antisense lncRNAs, 2251 (33.35%) lncRNAs, and 542 (8.00%) intronic RNAs were identified (Supplementary Table 1). In the skin tissue of patient with SS, 71 lncRNA transcripts (33 up-regulated and 38 down-regulated) and 2615 mRNA transcripts (950 up-regulated and 1565 down-regulated) were differentially expressed compared with the normal skin (Figure 2A, 2B). The differentially expressed lncRNAs and mRNAs that were previously selected were further screened by heat maps between the SS patients groups and the normal groups using unsupervised hierarchical clustering analysis (Figure 3A, 3B). The analysis revealed that the SS samples could be distinguished from the normal samples as a different heat signature was evident in each case (Figure 3A, 3B). It is interesting to note that this discrimination could be achieved by both lncRNA and/or mRNA screening, whereas the heat signatures were very similar in both cases with regard to the up-regulated and down-regulated genes (Figure 3A, 3B). The top 20 differentially expressed mRNA and lncRNAs are listed in Table 1 and 2.
Figure 2: Volcano plots of DE transcripts. The difference of lncRNAs expression profiles (A) and mRNAs expression profiles (B) can be noted in the overall distribution of the transcripts. The red points in the plot represent upregulated transcripts and the green points represent downregulated transcripts. The filter threshold is p value < 0.05 by default.
Figure 3: Hierarchical heat maps of the differentially expressed profiles between sensitive and normal skin. lncRNA (A) and mRNA (B) were used for analysis of the gene expression data, for which the cluster analysis arranges samples into groups according to their expression levels. Each column represents a sample and each row represents a transcript. ‘Red’ indicates higher expression, and ‘blue’ indicates lower expression. Sensitive skin (SS), normal skin (N).
Table 1: Top 20 differentially expressed mRNAs with >1.5-fold change in 3 sensitive skin (SS) compared with normal skin tissue (N)
Upregulated mRNAs | Downregulated mRNAs | ||||
|---|---|---|---|---|---|
Transcript ID (Ensembl_Gene_ID) | P-value | Fold change (SS vs. N) | Transcript ID (Ensembl_Gene_ID) | P-value | Fold change (SS vs. N) |
ENST00000502213 (TLR1) | 1.92E-13 | 6.11601 | ENST00000557022 (ZFP36L1) | 1.96E-12 | –5.76884 |
ENST00000283684 (ATP8B1) | 2.94E-14 | 6.04138 | ENST00000616053 (TCF4) | 6.53E-12 | –5.72168 |
ENST00000620565 (UHRF1) | 3.42E-12 | 5.80807 | ENST00000502252 (CORIN) | 2.00E-10 | –5.27057 |
ENST00000510411 (HNRNPH1) | 2.18E-12 | 5.48796 | ENST00000535987 (FOS) | 8.45E-10 | –5.0481 |
ENST00000336665 (AGAP1) | 5.80E-09 | 4.84147 | ENST00000417268 (SCAF8) | 5.84E-10 | –4.97343 |
ENST00000621364 (NOMO2) | 7.76E-10 | 4.70092 | ENST00000602390 (COMMD3-BMI1) | 3.10E-09 | –4.9622 |
ENST00000347088 (YY1AP1) | 3.58E-08 | 4.67080 | ENST00000264266 (MFSD1) | 3.32E-10 | –4.94658 |
ENST00000541717 (MELK) | 4.46E-08 | 4.56854 | ENST00000442173 (DOCK9) | 9.59E-10 | –4.8981 |
ENST00000450331 (PNPLA6) | 1.31E-08 | 4.49646 | ENST00000613065 (ZNF254) | 9.90E-09 | –4.84424 |
ENST00000360661 (BAK1) | 6.78E-08 | 4.48298 | ENST00000514633 (HNRNPAB) | 6.47E-09 | –4.81581 |
ENST00000399202 (FAM214A) | 1.08E-07 | 4.42187 | ENST00000520492 (ZFPM2) | 8.54E-09 | –4.80934 |
ENST00000380989 (PITRM1) | 2.01E-07 | 4.36727 | ENST00000614805 (PLXNB2) | 9.12E-09 | –4.7907 |
ENST00000396499 (CCDC125) | 1.78E-07 | 4.34897 | ENST00000532891 (ARHGAP27) | 6.45E-09 | –4.78249 |
ENST00000506339 (HNRNPAB) | 2.36E-07 | 4.32278 | ENST00000503781 (PIEZO2) | 1.46E-08 | –4.69927 |
ENST00000523714 (ANXA6) | 2.51E-07 | 4.318563 | ENST00000383083 (PLSCR4) | 2.14E-08 | –4.66819 |
ENST00000368488 (CEP85L) | 3.63E-08 | 4.23751 | ENST00000283921 (HOXA10) | 1.68E-08 | –4.64941 |
ENST00000342788 (ERBB4) | 1.81E-07 | 4.22746 | ENST00000526927 (LTBP3) | 1.84E-07 | –4.39143 |
ENST00000577278 (AXIN2) | 2.27E-06 | 3.98945 | ENST00000325468 (GYLTL1B) | 7.09E-08 | –4.38332 |
ENST00000382422 (BAZ1A) | 3.47E-06 | 3.92167 | ENST00000251268 (MEGF8) | 2.39E-08 | –4.25831 |
ENST00000412006 (GABBR1) | 3.52E-06 | 3.90662 | ENST00000355973 (CCDC77) | 4.83E-07 | –4.24682 |
P-values < 0.05 were considered significant.
Table 2: Top 20 differentially expressed lncRNAs with >1.5-fold change in 3 sensitive skin (SS) compared with normal skin tissue (N)
Upregulated lncRNAs | Downregulated lncRNAs | ||||
|---|---|---|---|---|---|
Transcript ID (Gene_Symbol) | P-value | Fold change (SS vs. N) | Transcript ID (Gene_Symbol) | P-value | Fold change (SS vs. N) |
ENST00000624863.1 (AC003973.3) | 1.36E-05 | 3.10833 | LNC_000101 | 0.000429 | –3.0039 |
ENST00000585940.1 (CTD-2537I9.12) | 0.000452 | 2.87593 | ENST00000589456.1 (CTD-2537I9.12) | 0.000831 | –2.86602 |
ENST00000441722.5 (ZFAS1) | 0.000339 | 2.33685 | ENST00000457658.5 (TTTY15) | 0.001591 | –2.70333 |
ENST00000438107.1 (RP11-706O15.3) | 0.001474 | 2.20917 | ENST00000564363.1 (RP11-1100L3.8) | 0.000563 | –2.5488 |
ENST00000449766.5 (AC016735.2) | 0.008838 | 2.19405 | ENST00000564531.1 (RP11-1100L3.8) | 0.00051 | –2.53449 |
LNC_000156 | 0.005353 | 2.17926 | ENST00000528549.1 (TBX5-AS1) | 0.001077 | –2.51396 |
LNC_000208 | 0.013481 | 2.12202 | ENST00000567257.5 (LOXL1-AS1) | 0.006557 | –2.34528 |
LNC_000082 | 0.011171 | 2.09752 | LNC_000265 | 0.008368 | –2.26627 |
ENST00000607746.1 (LINC01215) | 0.005734 | 2.09382 | ENST00000566193.1 (RP11-424G14.1) | 0.011503 | –2.17116 |
LNC_000244 | 0.018443 | 2.04553 | ENST00000331787.2 (TTTY14) | 0.01933 | –2.0051 |
LNC_000144 | 0.014573 | 2.03612 | ENST00000414790.5 (H19) | 0.008537 | –1.99685 |
ENST00000424210.1 (SPAG5-AS1) | 0.021094 | 1.98796 | ENST00000503474.5 (HAND2-AS1) | 0.020475 | –1.9804 |
ENST00000381108.3 (RP11-706O15.3) | 0.021268 | 1.974543 | ENST00000439725.5 (H19) | 0.017798 | –1.95243 |
ENST00000562127.1 (RP11-445O3.3) | 0.023283 | 1.91734 | ENST00000507362.5 (RP11-319G6.1) | 0.016665 | –1.9427 |
ENST00000500989.2 (LINC00861) | 0.006481 | 1.89280 | ENST00000553510.1 (RP11-293M10.1) | 0.023008 | –1.92564 |
ENST00000559298.5 (IQCH-AS1) | 0.025628 | 1.84106 | LNC_000084 | 0.022983 | –1.89912 |
ENST00000414002.5 (SNHG5) | 0.010083 | 1.78812 | LNC_000008 | 0.027705 | –1.88524 |
ENST00000355358.1 (GATA3-AS1) | 0.013621 | 1.76931 | ENST00000567089.1 (RP11-265N6.2) | 0.028664 | –1.88521 |
LNC_000165 | 0.026313 | 1.76915 | ENST00000620503.1 (RP11-115H13.1) | 0.027862 | –1.88198 |
ENST00000417483.5 (RP11-557H15.3) | 0.025022 | 1.76723 | LNC_000087 | 0.026409 | –1.87013 |
Genomic features of lncRNAs
The characteristics of the lncRNAs and their comparison with protein-coding genes were further examined. The predicted novel lncRNA were greater in length (Figure 4A) and contained fewer exons (Figure 4B) compared with the mRNA transcripts. Furthermore, we identified that lncRNAs were shorter in ORF length (Figure 4C) and appeared less conserved than mRNA transcripts (Figure 4D). These findings were in agreement with previous studies [32]. The predicted lncRNA transcripts indicated different genomic characteristics compared with protein coding transcripts.
Figure 4: Genomic features of predicted lncRNAs. (A) Length distribution of coding transcripts (blue), known lncRNAs (purple), and new predicted lncRNAs (red). (B) Exon number distribution of coding transcripts and lncRNAs. (C) Orf length distribution of coding transcripts and lncRNAs. (D) Conservation score for coding transcripts and lncRNAs as determined by the phastCons software.
Alternative splicing events
We calculated the expression amounts of AS events. The results are shown in Figure 5 and Table 3. The results showed that AS events are involved in SS.
Figure 5: Expression amounts of alternative splicing (AS) events and differential analysis using rMATS (http://rnaseq-mats.sourceforge.net/index.html). SE: skipped exon; MXE: mutually exclusive exons; A5SS: alternative 5′ splice site; A3SS: alternative 3′ splice site; RI: retained intron.
Table 3: The type and quantitative statistic analysis of AS
EventType | NumEvents.JC | SigEvents.JC | ||
|---|---|---|---|---|
Total | Upregulated | Downregulated | ||
SE | 28594 | 27 | 11 | 16 |
MXE | 2902 | 3 | 3 | 0 |
A5SS | 2211 | 2 | 1 | 1 |
A3SS | 3402 | 2 | 1 | 1 |
RI | 3406 | 18 | 0 | 18 |
P-value < 0.05 were considered significant
GO and KEGG pathway analysis
To establish the function and connection of differentially expressed mRNAs, GO and KEGG pathway enrichment analysis were conducted. In the GO analysis, the apparently enriched terms are shown in Figure 6A and Supplementary Table 2. The gene networks were mainly associated with the following biological process terms: system development, multicellular organism development, anatomical structure development, and developmental process, whereas with regard to the cellular component terms the networks appeared to be active both in the intracellular part and the organelles (Figure 6A). The main molecular functions involved in the gene networks were protein binding and binding (Figure 6A).
Figure 6: GO enrichment analysis of differentially expressed protein-coding genes. (A) The ordinate represents the next level GO term of the three categories of GO. The abscissa represents the gene number under the term. In the directed acyclic graph (DAG) the downstream term corresponds to a subset of the upstream term. The depth of the color represents the degree of enrichment (B–D). The DAG of over-represented biological process (B), cellular component (C), and molecular function (D) terms in GO analysis between sensitive skin and normal skin.
Directed acyclic graphs (DAGs) were constructed in order to highlight the associations among the enriched GO terms (Figure 6B–6D), in which the downstream terms are subsets of the upstream terms. The GO analysis demonstrated that KRT27 and CLDN5 genes, which are associated with keratinocyte differentiation and epidermal development, were highly expressed in the SS GO term “system development”. In addition, IL-27RA and CCL18, which are key cytokine and chemokine genes involved in inflammation, were specifically expressed in the SS GO term “protein binding and binding” (Figure 6A).
KEGG pathway analysis deduces the systematic biological behavior by mapping the protein-coding DEGs to the KEGG reference pathway. The results indicated that the transcripts with the lowest q value (i.e. lowest false discovery rate FDR) were involved in transcriptional misregulation in cancer, cancer-related pathways, focal adhesion, and ECM receptor interactions (Figure 7). Among these, the highest number of genes involved was reported for the cancer-related pathways (Figure 7). A considerable high number of genes were enriched into the PI3K-Akt signaling pathway, although the q value was somewhat higher than that noted for the aforementioned biological processes (Figure 7).
Figure 7: Top 20 over-represented KEGG pathways of the common differentially expressed genes. The color coding (red, orange, green, light blue, purple) represents the q value i.e. the lowest false discovery rate (FDR), whereas the size of the dots represents the number of the genes involved in each pathway.
qRT-PCR validation
We randomly selected five transcripts (LNC_000265, ENST00000441722.5, ENST00000414790.5, ENST000 00624863.1, and ENST00000413119) for qRT-PCR in 10 SS and 10 normal skin samples in order to validate the data derived by RNA-seq. The results yielded a similar trend with regard to up- and/or down-regulation for each transcript compared with the RNA-Seq data. With regard to the up-regulated transcripts, RNA-seq appeared to be more sensitive compared with qRT-PCR (Figure 8). The data of the qRT-PCR assays were in agreement with the RNA-seq results and confirmed the reliability and the validity of the RNA-Seq analysis (Figure 8).
Figure 8: Relative expression analysis of candidate mRNA and lncRNAs in sensitive and normal skin tissues, as determined by qRT-PCR. (A) Validation of five candidate transcripts in 13 sensitive skin and normal skin tissues. The results were consistent. (B) The heights of the columns represent the log2 (fold changes). The trends of the qRT-PCR results were consistent with the RNA-seq data.
Predicted function of lncRNA
Numerous lncRNAs regulate their target protein-coding genes in a trans fashion. LNC_000265 (LNC_000265, ENST00000454741.5) was selected for further analysis. The trans-regulating target mRNAs of LNC_000265 that exhibited a PCC score higher than 0.98 were used to construct the network using the Cytoscape 3.2 software (Figure 9). The connections between lncRNAs and mRNAs were established through nodes and edges on the diagram. The complex network indicates that the same lncRNA can simultaneously control multiple protein-coding genes through antisense regulation, and the same protein-coding gene can further be regulated by multiple lncRNAs concomitantly.
Figure 9: Co-expression network for one annotated and one novel lncRNA. Two selected lncRNAs were co-expressed with 318 coding gene transcripts. The expression of LNC_000265 correlated with 153 coding gene transcripts (Pearson Correlation Coefficient >0.98).
DISCUSSION
The present study provided a comprehensive bioinformatics analysis regarding the expression of lncRNAs and mRNA transcripts in the skin tissue of subject with SS compared with normal skin. A total of 33 and 950 lncRNAs and mRNAs were up- regulated, whereas 38 and 1565 lncRNAs and mRNAs were down-regulated. Hierarchical clustering analysis revealed that the heat signature of the pathological tissue samples could be distinguished from the normal skin samples, whereas the majority of the genes that were present in enriched pathways were those that participated in focal adhesion, PI3K-Akt signaling, and cancer-related pathways. The data suggest that this analysis can yield hits in the discovery of key target genes involved in the development of SS.
The investigation of the key genes that are involved in SS has not been extensively reported to date possibly due to the lack of a definitive medical diagnostic consensus for this broad category of diseases [33]. Nevertheless, the identification of genes involved in skin diseases that are similar in pathology and can be medically defined has been previously reported. Diseases such as psoriasis, atopic dermatitis, and atopic eczema have been investigated in terms of bioinformatics and qPCR analyses [30, 34–39]. Although evidence regarding the aforementioned diseases from bioinformatics studies is limited, experimental studies have shown the contribution of certain genes in the development of skin pathophysiological diseases, namely genes involved in lipid homeostasis of the skin barrier (CERS4, [35]), genes involved in cell adhesion, migration, endocytosis and skin barrier (SDC-4, [37]), and genes that participate in the inflammatory response (IL-17, IL-20, IL-24, IL-31, and IL-33, [41, 43]. For example SDC-4 (syndecan-4) mRNA levels were increased in atopic dermatitis compared with normal skin samples [37] and the inflammatory mediators IL-17, IL-20, IL-24, IL-31, and IL-33 were also elevated in atopic dermatitis [34, 36]. The aforementioned experimental findings were confirmed by a bioinformatics study that identified 438 and 779 DEGs that were enriched in pathways involved in epidermis development and immune response, respectively [30], whereas a study that used human DNA microarrays revealed similar findings with regard to the participation of DEGs in immune response, lipid homeostasis, and epidermal differentiation in atopic eczema [40]. The present study identified 71 lncRNA transcripts (33 up-regulated and 38 down-regulated) and 2615 mRNA transcripts (950 up-regulated and 1565 down-regulated) that were differentially expressed in skin tissue of subject with SS compared with normal skin, which is in agreement with the studies by Ding et al. and Zhang et al. that demonstrated similar trends in the number of DEGs [30, 39]. In addition, GO analysis demonstrated that the KRT27 and CLDN5 genes, which are associated with keratinocyte differentiation and epidermal development, were highly expressed in the SS GO term “system development”, whereas IL-27RA and CCL18, which are key cytokine and chemokine genes involved in inflammation, were specifically expressed in the SS GO term “protein binding and binding”. These findings are in concordance with the studies of Ding et al., Zhang et al., and Nishiyama et al., in which the main pathways enriched were those corresponding to the inflammatory response and epidermis development [30, 39, 40].
With regard to the KEGG pathway analysis, the data indicated that the transcripts with the lowest q value were involved in transcriptional misregulation in cancer, cancer-related pathways, focal adhesion, and ECM receptor interactions. Importantly, the PI3K-Akt signaling pathway was also involved. The study by Zhang et al. identified a similar target gene, namely PI3KR1 (phosphoinositide-3-kinase regulatory subunit 1) in the chemokine signaling pathway that was the most enriched with regard to DEGs in atopic dermatitis [39]. Moreover, LAMA5, ITGB4, and other protein-coding genes associated with epidermal homeostasis were enriched in both pathways [41, 42]. Certain studies have suggested that the epidermal barrier function of SS is disrupted [43, 44]. Therefore, ECM-receptor interaction signaling and the PI3K-Akt signaling pathway may play an important role in the pathogenesis of SS. Furthermore, it is important to note that the pathways of cancer and inflammation are linked. For example, tobacco smoke acts as a tumor promoter by virtue of its high carcinogenic compound content can further trigger chronic inflammation [45]. Moreover, the majority of solid malignancies that appear in older individuals and cell senescence [46] are postulated to be tumor promoters that act through inflammatory mechanisms. In addition to its pro-tumorigenic effects, inflammation influences the host immune response to tumors and can be targeted to augment the response to chemotherapy [47]. Thus, it may be speculated that pathways involved in cancer may cross-talk with pathways involved in the inflammatory response, which provides an explanation for the results obtained by KEGG analysis in the present study.
Although the aforementioned studies have notably focused on the differential expression of mRNA transcripts, the contribution of lncRNAs in skin-associated diseases remains less studied. Ahn et al. demonstrated [31] that more than 50% of coexpressed genes, identified in three network modules that correlated with psoriasis and six network modules that correlated with biological treatment, were lncRNAs [31]. The current study identified the lncRNA LNC_000265 by RNA-seq and qRT-PCR analysis. This lncRNA was used to construct networks with genes that could be potentially correlated with regard to their transcripts. LNC_000265 and its predicted trans-regulated target CLDN5 were highly correlated. The protein Claudin5 encoded by CLDN5 plays a vital role in the epidermal barrier structure of tight junctions (TJs). TJs are intercellular contacts that seal the space between the individual cells of an epithelial sheet so as to collectively separate tissue compartments. In comparison with other junctions, TJs are composed of approximately 40 different proteins including claudins (Cldn) and zonula occludens (ZO), which are known to be important for barrier function in epithelial cells [48]. The data indicated that CLDN5 was expressed at low levels in skin tissue of subject with SS compared with normal skin (P < 0.05, log2FoldChange = –1.93758). Furthermore, LNC_000265 was closely related to its predicted trans-regulated target CLDN5 (PCC= 0.98). The results of the qRT-PCR assay demonstrated that LNC_000265 and CLDN5 were also expressed at low levels in skin tissues of subjects with SS. This gene cluster further included CLDN1, CLDN2, CLDN6, and CLDN12, and exhibited a strong correlation with TJs [49]. The findings suggest that LNC_000265 may play a role in the epidermal barrier structure of TJs by regulating the expression level of CLDN5.
Despite the aforementioned findings that expand the existing knowledge on the molecular gene network interactions that are observed in SS, the present study has several limitations. Although the data presented were consistent with previously published reports, the overall sample size was considerably small and only females were included, which could induce bias in the conclusions of the study. Additional bias is caused by the application of a gene set enrichment analysis that is based on the rational that gene networks are constructed based on network or pathway information. In contrast to gene set enrichment analysis, weighted gene coexpression network analysis (WGCNA) can be used in future studies to prioritize the main genes in a given network by the estimation of the connectivity of each gene as it has been noted by previous studies [31]. A WGCNA-based screen can achieve a higher validation rate compared with a differential expression analysis in biomarkers of glioblastoma [49]. Furthermore, the present study used sensitive skin tissue, in the absence of the isolation of specific cell types from the tissues. Future studies are needed to address the contribution of each gene network to a particular cell type, whereas the verification of the target genes can be further validated by the use of appropriate KO cell and animal models. Finally, no direct experimental evidence was provided to confirm the results of lncRNA function prediction. Our future research will involve the construction of appropriate cell models in order to address these limitations.
In conclusion, the present study provided a comprehensive analysis of the DEGs with regard to lncRNA and mRNA transcripts in skin tissue of subject with SS compared with normal skin by RNA-seq and bioinformatics analysis. This study indicated that lncRNAs contribute a critical role in the pathogenesis of SS, which may aid the identification of the pathogenesis of SS and the development of potential targets for personalized therapeutic strategies.
METHODS
Subjects and samples
A total of 26 patients who were scheduled to undergo facial nevus resection surgery at the dermatology department of the First Affiliated Hospital of Kunming Medical University were enrolled. The skin sensitivity was evaluated according to a questionnaire (Table 4) [50] and a lactic acid stinging test [51]. The subjects that responded positive to a minimum of 5 out of 7 questions and were associated with a stinging test score greater than or equal to 3 were classified as SS. The exclusion criteria were defined as the incidence of SS that resulted from acne, pityriasis rosea, contact dermatitis, eczema, and/or other skin diseases (Table 5). The skin tissue was obtained by trimming around the nevus. The tissue was mixed with RNA-later solution (QIAGEN, Germany) and stored at -80°C. The study protocol was approved by the Ethics Committee Board of the First Affiliated Hospital of Kunming Medical University. Informed written consent was obtained from each participant.
Table 4: Questionnaire for diagnosis of SS
1 | Would you say that your face/neck does not tolerate cold/hot weather or a cold/hot environment? |
2 | Would you say that your skin face/neck does not tolerate rapid temperature changes? |
3 | Have you already avoided the use of some cosmetic products that could, according to you, make your skin reactive? |
4 | Have you already had an adverse reaction on your face/neck to a cosmetic or hygiene product? |
5 | Would you say that your face/neck is reactive? |
6 | Have you already felt some itching, burning or tingling on your face/neck skin because of the wind or some cosmetics or hygiene products? |
7 | Is your face skin reactive to pollution, stress/emotions or menstrual cycle changes? |
Table 5: Summary of subjects characteristics
Characteristics | Sensitive skin (SS) | Normal skin (N) |
|---|---|---|
Total number | 13 | 13 |
Gender | Female | Female |
Age (y), mean ± SD | 38.5 ± 7.7 | 38.0 ± 7.6 |
Results of questionnaire, mean ± SD* | 6.3 ± 0.8 | 1.0 ± 0.7 |
Scores of lactic acid stinging test, mean ± SD* | 4.9 ± 0.9 | 1.6 ± 0.2 |
RNA sequencing | 3 | 3 |
Age (y), mean ± SD | 40.3 ± 6.7 | 39.3 ± 8.1 |
qRT-PCR | 10 | 10 |
Age (y), mean ± SD | 38 ± 8.2 | 37.6 ± 7.9 |
*P < 0.05 were considered significant by t-test
Library preparation and sequencing process
Total RNA was extracted from the pathological samples (n = 3) and the normal skin samples (n = 3). The extraction was carried out using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s instructions. The RNA quality was assessed using a Nano Photometer spectrophotometer (IMPLEN, CA, USA). The RNA concentration was measured using the Qubit RNA Assay Kit and the Qubit 2.0 Flurometer (Life Technologies, CA, USA). The RNA integrity was verified on an Agilent Bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA). The RNA library was constructed using a total amount of 3 μg of RNA per sample with an RNA integrity number (RIN) over 7.0. The Epicentre Ribo-zero rRNA Removal Kit (Epicentre, USA) was used to remove Ribosomal RNA (rRNA) according to the manufacturer’s protocol. Finally, the dUTP method was used to obtain strand-specific sequencing libraries by the NEB Next Ultra Directional RNA Library Prep Kit for Illumina (NEB, USA), according to the instructions provided by the manufacturer. The RNA-seq assay was carried out on an Illumina Hiseq 2000 platform and 100 bp paired-end reads were obtained. The preparation of the total transcriptome libraries was conducted by Novogene Bioinformatics Technology Cooperation (Beijing, China).
Reads mapping and transcriptome assembly
The filtered reads were compared with the human genome (hg19) using the TopHat2 software (v2.0.9) [52]. Based on the mapping results derived by TopHat2, the transcriptome was assembled and quantitatively analyzed by the Cufflinks software (v2.1.1) [53].
Differential expression analysis
Sequenced reads (raw data or raw reads) in the FASTQ format were first processed via in-house Perl scripts. Clean reads were prepared by deleting reads containing ploy-N, adapter, and low quality reads from raw data. GC content, Q20, and Q30 of the cleaned data were calculated. The false discovery rate (FDR) was calculated as described previously [54] and a FDR of 1% was used. Clean reads were mapped to the reference genome built with Bowtie (v2.0.6) [55]. The mapped reads were assembled by Cufflinks (v2.1.1) in a reference-based approach [53]. The DEGs were identified at a confidence interval of 95% for each experiment. The alternative splicing (AS) events were analyzed using rMATS (http://rnaseq-mats.sourceforge.net/index.html). In the case of a comparison group with a differential variable splice analysis, we first counted the types and quantities of variable splice events, and then calculated the amount of variable splice events for each category, and finally for each variable splice events for differential analysis.
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis
Differentially expressed transcripts were analyzed by GO enrichment and KEGG pathway in order to identify possible enrichment of genes with specific biological pathways (http://www.genome.jp/kegg/pathway.html). GO (http://www.geneontology.org/) is a database that provides vocabularies and classifications associated with the molecular and cellular structures and functions regarding the biological annotations of genes [56]. GO terms consist of three categories: biological processes (BP), cellular component (CC), and molecular function (MF). The KEGG (http://www.genome.jp/kegg/ or http://www.kegg.jp/) database includes information regarding the contribution of DEGs in known metabolic pathways and regulatory pathways, and facilitates the mapping of genes to KEGG pathways for systemic analysis of gene function. GO and KEGG pathway enrichment analysis were conducted for the upregulated and downregulated DEGs, respectively. KEGG analysis was used to see in which pathways the DEGs participated. If multiple transcripts are enriched in a specific pathway, the q value for the enrichment of the pathway was calculated, rather than the q value of each transcripts.
GO enrichment was performed with the GO Seq R package, with GO terms with P < 0.05 being considered as significantly enriched. The KOBAS (KO Based Annotation system) software (Beijing, China) was used to assess statistical enrichment in KEGG pathways. GO enrichment indicates the relationship between genes and GO terms. For each gene g and each GO term GOj, a score is generated, which is typically referred to as the gene ontology enrichment score. The enrichment score refers to the ratio of the number of genes located in the pathway entry in the differentially expressed genes to the total number of genes in the annotated gene located in the pathway entry (or “enrichment factor” = “the number of candidate genes”/”the number of background genes”). The greater the enrichment score, the greater the degree of enrichment. The q value is the value of the P value after multiple hypothesis test correction.
Prediction of the function of lncRNAs
The function of most LncRNAs in current databases is unclear. The biological functions of lncRNAs were analyzed through the position of target protein coding genes (cis) and expression correlation (trans). We set the co-location threshold to lncRNA upstream and downstream 100 kb, and analyzed the co-expression pairs by calculating Pearson correlation coefficients (PCC).
Quantitative real-time PCR (qRT-PCR)
Total RNA isolated from 10 SS patients and 10 normal facial tissues using Trizol reagent (Invitrogen, USA) was reverse transcribed to cDNA using All-in-One First- Strand cDNA Synthesis Kit (Gene Copoeia, USA), according to the manufacturer’s protocol. Two significantly upregulated (ENST00000441722.5 and ENST00000624863.1) and two downregulated lncRNAs (LNC-000265 and ENST00000414790.5) were selected to test the reliability of the analysis. In addition, one mRNA transcript (ENST00000413119) was selected to validate the result. The 2-ΔΔCt method was used with β-actin as an internal control in order to relatively quantify the detected transcripts. The primer sequences are listed in Table 6.
Table 6: Primer pairs used for qRT-PCR experiments
Gene ID | Primer sequence |
|---|---|
ENST00000413119 | Forward 5′- CCGAGTCGTACACTTTGCAC -3′ |
Reverse 5′- CCTTCCTGGACCACAACATC -3′ | |
ENST00000441722.5 | Forward 5′- CAGTCAGGAGAAAGAAGTGGAGG -3′ |
Reverse 5′- ACAAATAAGAGGGGACAGAGGTG -3′ | |
ENST00000624863.1 | Forward 5′- GGATGGGAGACAAGCATAGAAAAT -3′ |
Reverse 5′- TGTGAGGAGACCTGGTATAGAAAC -3′ | |
ENST00000414790.5 | Forward 5′- ACGAGTGTGCGTGAGTGTGAG -3′ |
Reverse 5′- ATTGATGATGAGTCCAGGGCT -3′ | |
LNC-000265 | Forward 5′- CCTTCCCTGATGTCTGATTTTTG -3′ |
Reverse 5′- GCCTCTTCTCCCATTTGTTTTTC -3′ | |
β-Actin | Forward 5′- CCAGGGCGTTATGGTAGGCA -3′ |
Reverse 5′- TTCCATATCGTCCCAGTTGGT -3′ |
ACKNOWLEDGMENTS AND FUNDING
The present study was supported by the Internal research institutions research project of Yunnan Province medical and health Institutions (2017NS005) and the National Natural Science Foundation of China (U1402223, 81360234 and 81460469). The Science and Technology Leading Talent Project of Yunnan province (grant no. 2017HA010). The funders had no role in the study design, data collection and analysis, decision to publish, and/or preparation of the manuscript.
CONFLICTS OF INTEREST
The authors declare that they have no competing interests.
REFERENCES
1. Misery L, Loser K, Stander S. Sensitive skin. J Eur Acad Dermatol Venereol. 2016; 30:2–8.
2. Berardesca E, Fluhr JW, Maibach H. (2006). What is sensitive skin? Dermatology: Clinical and Basic Science Series, Sensitive Skin syndrome. (New York: Taylor and Francis).
3. Lev-Tov H, Maibach HI. The sensitive skin syndrome. Indian journal of dermatology. 2012; 57:419–423.
4. Marriott M, Holmes J, Peters L, Cooper K, Rowson M, Basketter D. The complex problem of sensitive skin. Contact dermatitis. 2005; 53:93–99.
5. Roussaki-Schulze AV, Zafiriou E, Nikoulis D, Klimi E, Rallis E, Zintzaras E. Objective biophysical findings in patients with sensitive skin. Drugs Exp Clin Res. 2005; 31:17–24.
6. Saint-Martory C, Roguedas-Contios AM, Sibaud V, Degouy A, Schmitt AM, Misery L. Sensitive skin is not limited to the face. Br J Dermatol Syph. 2008; 158:130–133.
7. Spizzo R, Almeida MI, Colombatti A, Calin G. Long non-coding RNAs and cancer: a new frontier of translational research? Oncogene. 2012; 31:4577–4587.
8. Blackshaw S, Harpavat S, Trimarchi J, Cai L, Huang H, Kuo WP, Weber G, Lee K, Fraioli RE, Cho SH, Yung R, Asch E, Ohno-Machado L, et al. Genomic analysis of mouse retinal development. PLoS Biol. 2004; 2:E247.
9. Dinger ME, Amaral PP, Mercer TR, Pang KC, Bruce SJ, Gardiner BB, Askarian-Amiri ME, Ru K, Solda G, Simons C, Sunkin SM, Crowe ML, Grimmond SM, et al. Long noncoding RNAs in mouse embryonic stem cell pluripotency and differentiation. Genet Genome. 2008; 18:1433–1445.
10. Rinn JL, Kertesz M, Wang JK, Squazzo SL, Xu X, Brugmann SA, Goodnough LH, Helms JA, Farnham PJ, Segal E, Chang H. Functional demarcation of active and silent chromatin domains in human HOX loci by noncoding RNAs. Cell. 2007; 129:1311–1323.
11. Szymanski M, Barciszewska MZ, Erdmann VA, Barciszewski J. A new frontier for molecular medicine: noncoding RNAs. Biochim Biophys Acta. 2005; 1756:65–75.
12. Clemson CM, Hutchinson JN, Sara SA, Ensminger AW, Fox AH, Chess A, Lawrence J. An architectural role for a nuclear noncoding RNA: NEAT1 RNA is essential for the structure of paraspeckles. Mol Cell. 2009; 33:717–726.
13. Hutchinson JN, Ensminger AW, Clemson CM, Lynch CR, Lawrence JB, Chess A. A screen for nuclear transcripts identifies two linked noncoding RNAs associated with SC35 splicing domains. BMC Genomics. 2007; 8:39.
14. Mercer TR, Dinger ME, Sunkin SM, Mehler MF, Mattick J. Specific expression of long noncoding RNAs in the mouse brain. Proc Natl Acad Sci U S A. 2008; 105:716–721.
15. Ravasi T, Suzuki H, Pang KC, Katayama S, Furuno M, Okunishi R, Fukuda S, Ru K, Frith MC, Gongora MM, Grimmond SM, Hume DA, Hayashizaki Y, Mattick J. Experimental validation of the regulated expression of large numbers of non-coding RNAs from the mouse genome. Genet Genome. 2006; 16:11–19.
16. Sasaki YT, Ideue T, Sano M, Mituyama T, Hirose T. MENepsilon/beta noncoding RNAs are essential for structural integrity of nuclear paraspeckles. Proc Natl Acad Sci U S A. 2009; 106:2525–2530.
17. Sone M, Hayashi T, Tarui H, Agata K, Takeichi M, Nakagawa S. The mRNA-like noncoding RNA Gomafu constitutes a novel nuclear domain in a subset of neurons. J Cell Sci. 2007; 120:2498–2506.
18. Kung JT, Colognori D, Lee J. Long noncoding RNAs: past, present, and future. Genetics. 2013; 193:651–669.
19. Bonasio R, Shiekhattar R. Regulation of transcription by long noncoding RNAs. Annu Rev Genet. 2014; 48:433–455.
20. Weikard R, Hadlich F, Kuehn C. Identification of novel transcripts and noncoding RNAs in bovine skin by deep next generation sequencing. BMC Genomics. 2013; 14:789.
21. Cabili MN, Trapnell C, Goff L, Koziol M, Tazon-Vega B, Regev A, Rinn J. Integrative annotation of human large intergenic noncoding RNAs reveals global properties and specific subclasses. Genes Dev. 2011; 25:1915–1927.
22. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nature reviews Genetics. 2009; 10:57–63.
23. He Y, Meng XM, Huang C, Wu BM, Zhang L, Lv XW, Li J. Long noncoding RNAs: Novel insights into hepatocelluar carcinoma. Cancer Lett. 2014; 344:20–27.
24. Martens-Uzunova ES, Bottcher R, Croce CM, Jenster G, Visakorpi T, Calin G. Long noncoding RNA in prostate, bladder, and kidney cancer. Eur Urol. 2014; 65:1140–1151.
25. Orom UA, Derrien T, Beringer M, Gumireddy K, Gardini A, Bussotti G, Lai F, Zytnicki M, Notredame C, Huang Q, Guigo R, Shiekhattar R. Long noncoding RNAs with enhancer-like function in human cells. Cell. 2010; 143:46–58.
26. Mazar J, Sinha S, Dinger ME, Mattick JS, Perera R. Protein-coding and non-coding gene expression analysis in differentiating human keratinocytes using a three-dimensional epidermal equivalent. Mol Genet Genomics. 2010; 284:1–9.
27. Hombach S, Kretz M. The non-coding skin: exploring the roles of long non-coding RNAs in epidermal homeostasis and disease. Bioessays. 2013; 35:1093–1100.
28. Kretz M, Webster DE, Flockhart RJ, Lee CS, Zehnder A, Lopez-Pajares V, Qu K, Zheng GX, Chow J, Kim GE, Rinn JL, Chang HY, Siprashvili Z, Khavari P. Suppression of progenitor differentiation requires the long noncoding RNA ANCR. Genes and development. 2012; 26:338–343.
29. Sonkoly E, Bata-Csorgo Z, Pivarcsi A, Polyanka H, Kenderessy-Szabo A, Molnar G, Szentpali K, Bari L, Megyeri K, Mandi Y, Dobozy A, Kemeny L, Szell M. Identification and characterization of a novel, psoriasis susceptibility-related noncoding RNA gene, PRINS. J Biol Chem. 2005; 280:24159–24167.
30. Ding Y, Shao X, Li X, Zhai Y, Zhang Y, Wang S, Fang H. Identification of candidate genes in atopic dermatitis based on bioinformatic methods. Int J Dermatol. 2016; 55:791–800.
31. Ahn R, Gupta R, Lai K, Chopra N, Arron ST, Liao W. Network analysis of psoriasis reveals biological pathways and roles for coding and long non-coding RNAs. BMC Genomics. 2016; 17:841.
32. Qu Z, Adelson DL. Bovine ncRNAs are abundant, primarily intergenic, conserved and associated with regulatory genes. PLoS One. 2012; 7:e42638.
33. Richters R, Falcone D, Uzunbajakava N, Verkruysse W, van Erp P, van de Kerkhof P. What is sensitive skin? A systematic literature review of objective measurements. Skin Pharmacol Physiol. 2015; 28:75–83.
34. Cornelissen C, Marquardt Y, Czaja K, Wenzel J, Frank J, Luscher-Firzlaff J, Luscher B, Baron J. IL-31 regulates differentiation and filaggrin expression in human organotypic skin models. J Allergy Clin Immunol Pract. 2012; 129:426–433. 433 e421–428.
35. Ito S, Ishikawa J, Naoe A, Yoshida H, Hachiya A, Fujimura T, Kitahara T, Takema Y. Ceramide synthase 4 is highly expressed in involved skin of patients with atopic dermatitis. J Eur Acad Dermatol Venereol. 2017; 31:135–141.
36. Ma L, Xue HB, Wang F, Shu CM, Zhang J. MicroRNA-155 may be involved in the pathogenesis of atopic dermatitis by modulating the differentiation and function of T helper type 17 (Th17) cells. Am J Clin Exp Immunol. 2015; 181:142–149.
37. Nakao M, Sugaya M, Takahashi N, Otobe S, Nakajima R, Oka T, Kabasawa M, Suga H, Morimura S, Miyagaki T, Fujita H, Asano Y, Sato S. Increased syndecan-4 expression in sera and skin of patients with atopic dermatitis. Arch Dermatol Res. 2016; 308:655–660.
38. Trzeciak M, Sakowicz-Burkiewicz M, Wesserling M, Glen J, Dobaczewska D, Bandurski T, Nowicki R, Pawelczyk T. Altered Expression of Genes Encoding Cornulin and Repetin in Atopic Dermatitis. Int Arch Allergy Immunol. 2017; 172:11–19.
39. Zhang ZK, Yang Y, Bai SR, Zhang GZ, Liu TH, Zhou Z, Wang CM, Tang LJ, Wang J, He S. Screening for key genes associated with atopic dermatitis with DNA microarrays. Mol Med Rep. 2014; 9:1049–1055.
40. Nishiyama T, Amano S, Tsunenaga M, Kadoya K, Takeda A, Adachi E, Burgeson R. The importance of laminin 5 in the dermal-epidermal basement membrane. J Dermatol Sci. 2000; 24:S51–59.
41. Oldak M, Smola-Hess S, Maksym R. Integrin beta4, keratinocytes and papillomavirus infection. Int J Mol Med. 2006; 17:195–202.
42. Raj N, Voegeli R, Rawlings AV, Doppler S, Imfeld D, Munday MR, Lane M. A fundamental investigation into aspects of the physiology and biochemistry of the stratum corneum in subjects with sensitive skin. Int J Cosmet Sci. 2017; 39:2–10.
43. Pinto P, Rosado C, Parreirao C, Rodrigues L. Is there any barrier impairment in sensitive skin?: a quantitative analysis of sensitive skin by mathematical modeling of transepidermal water loss desorption curves. Skin Res Technol. 2011; 17:181–185.
44. Grivennikov SI, Greten FR, Karin M. Immunity, inflammation, and cancer. Cell. 2010; 140:883–899.
45. Rodier F, Coppe JP, Patil CK, Hoeijmakers WA, Munoz DP, Raza SR, Freund A, Campeau E, Davalos AR, Campisi J. Persistent DNA damage signalling triggers senescence-associated inflammatory cytokine secretion. Nat Cell Biol. 2009; 11:973–979.
46. Zitvogel L, Apetoh L, Ghiringhelli F, Kroemer G. Immunological aspects of cancer chemotherapy. Nature reviews Immunology. 2008; 8:59–73.
47. Brandner JM. Importance of Tight Junctions in Relation to Skin Barrier Function. Curr Probl Dermatol. 2016; 49:27–37.
48. Furuse M, Hata M, Furuse K, Yoshida Y, Haratake A, Sugitani Y, Noda T, Kubo A, Tsukita S. Claudin-based tight junctions are crucial for the mammalian epidermal barrier: a lesson from claudin-1-deficient mice. J Cell Biol. 2002; 156:1099–1111.
49. Horvath S, Zhang B, Carlson M, Lu KV, Zhu S, Felciano RM, Laurance MF, Zhao W, Qi S, Chen Z, Lee Y, Scheck AC, Liau LM, et al. Analysis of oncogenic signaling networks in glioblastoma identifies ASPM as a molecular target. Proc Natl Acad Sci U S A. 2006; 103:17402–17407.
50. Buhe V, Vie K, Guere C, Natalizio A, Lheritier C, Le Gall-Ianotto C, Huet F, Talagas M, Lebonvallet N, Marcorelles P, Carre JL, Misery L. Pathophysiological Study of Sensitive Skin. Acta Derm Venereol. 2016; 96:314–318.
51. Frosch PJ, Kligman AM. A method of appraising the stinging capacity of topically applied substances. J Soc Cosmet Chem. 1977; 28:197–209.
52. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, Salzberg S. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 2013; 14:R36.
53. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010; 28:511–515.
54. Lambrou GI, Zaravinos A, Adamaki M, Spandidos DA, Tzortzatou-Stathopoulou F, Vlachopoulos S. Pathway simulations in common oncogenic drivers of leukemic and rhabdomyosarcoma cells: a systems biology approach. Int J Oncol. 2012; 40:1365–1390.
55. Langmead B, Trapnell C, Pop M, Salzberg S. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009; 10:R25.
56. Liu F, Ji F, Ji Y, Jiang Y, Sun X, Lu Y, Zhang L, Han Y, Liu X. In-depth analysis of the critical genes and pathways in colorectal cancer. Int J Mol Med. 2015; 36:923–930.
