Transcriptional changes in organoculture of full-thickness human skin following topical application of all-trans retinoic acid.

In this study, we developed an organoculture of human skin to investigate the effect of topical applied all-trans retinoic acid using a gene array approach. We could by using this approach confirm previous studies on genes activated by RA in keratinocyte monocultures and also provide new insights on genes that are relevant to RA-activation in human skin. The results in the present study show this model represent a valuable pre-clinical model for studying the effects of retinoids in skin.

Introduction

One of the most studied classes of skin-targeted compounds is retinoids. All-trans retinoic acid (RA) also referred to as tretinoin, and its synthetic or natural derivatives (retinoids) affect epidermal growth and differentiation 1. Synthetic retinoids are today used widely in the treatment of psoriasis and other disorders of keratinization 2. Also, acne has been widely treated with various forms of natural and synthetic retinoids for more than 40 years 3.

Later, retinoids has been used in the treatment for hyperpigmentation disorders such as post-inflammatory hyperpigmentation and melasma 4.

Separate from the effect of retinoids in skin disease, Kligman and Willis were the first to introduce retinoids for use as antiphotoaging agents 5. After its application, the authors noticed improvement in skin depigmentation and rejuvenation 6. Today, retinoids are extensively used for this indication 7,8.

Further studies have shown that RA's clinical effects include improvement in wrinkles, surface roughness, mottled pigmentation and skin appearance as a whole when used on photo-damaged skin 9–11. However, topical use of RA has also shown to develop erythema, scaling and burning/pruritus 12. Therefore, there is a constant thrive within the cosmetic industry to develop retinoid-like anti-ageing ingredients with significantly reduced side effects.

At a cellular level, retinoids are known to modulate the proliferation and differentiation of epidermal keratinocytes 13,14,14,15 by binding to nuclear retinoic acid receptors (RARs) and 9-cis-retinoic acid receptors (RXRs), which results in upregulation or downregulation of the transcription of target genes 16.

RA is the biologically active form of vitamin A (retinol). Endogenous production of RA by epidermal keratinocytes involves the uptake of preformed vitamin A from the surroundings followed by a series of metabolic activation steps involving retinol dehydrogenases (RDH1, RDH4, RDH10, RDH12 and DHRS9) and retinal dehydrogenases (RALDH1, RALDH2 and RALDH3) 17–19. DHRS9 specifically is encoding an enzyme that mediates conversion of retinol into RA.

Whilst endogenously produced or exogenously applied, RA is transported through the cytoplasm by specific intracellular retinoic-binding proteins (CRABP2) 8.

Once synthesized, the cellular levels of RA are controlled by several cytochrome P450-dependent enzymes (CYPs) – CYP26 A1, B1 and C1 – which metabolize RA into 4-hydroxy-RA, 4-oxo-RA and 18-OH-RA 20,21. CYP26B1 is known to be higher expressed than both CYP26A1 and CYP26C1 in human keratinocytes 19,22.

Over the last quarter century, 532 genes have been validated to be regulated by RA 23. In some cases, this control is direct, driven by a liganded RAR:RXR heterodimer bound to a DNA response element; in others, it is indirect, reflecting the actions of intermediate transcription factors, non-classical associations of receptors with other proteins or even more distant mechanisms.

Cell cultures have been studied extensively to characterize and prophesy the effect of retinoids on epidermal differentiation and growth 13,15,24–27.

In addition, a selection of genes using RT-PCR has been studied in full-thickness skin treated with retinoids 28,29.

Today, recent gene array techniques allow the characterization of the mRNA expression status of a large number of genes in cells or tissues after retinoid treatments with more than 170 studies conducted.

Some studies using gene arrays have also been carried out on reconstructed epidermis and human epidermis in vivo 30.

However, to our knowledge, gene array analysis following the treatment of RA on organoculture of full-thickness skin has not been investigated earlier.

The aim of the present study was to determine whether it is possible to use an organoculture of full-thickness skin as a preclinical model to evaluate modulation of gene expression profile following topical application of a commercially available retinoid cream.

Materials and methods

Ex vivo study – skin explant model

Skin from surgical waste material from breast reductions was collected and immediately put in DMEM at 8°C. All subcutaneous fat was removed with a scalpel to ensure full diffusion nutrification of the tissue. 8-mm punch biopsies were taken from the skin and referred to as full-thickness skin explants. Thereafter, the explants were put on Millipore Millicell Culture Plate Inserts (12 mm Ø) (Millipore Corporation, Stockholm, Sweden). Inserts containing skin explants were put in 6-well plates (1 insert/well), and 1 mL of supplemented keratinocyte medium (M154) (Life Technologies, Stockholm, Sweden) was added to each well to allow survival of the explants.

Three donors were used for this study. They were all healthy Caucasian females of 23, 41 and 60 years of age. The collection of skin samples was approved by the local Ethics Committee of Stockholm. All participants gave written informed consent for the skin sample collections, and the studies followed the guidelines of the World Medical Association's Declaration of Helsinki (2000) concerning biomedical research involving human subjects.

Exhibiting an enhanced functional barrier function compared with reconstructed skin 31 (ref), the use of viable skin explant allows topical application of finished product to study cosmetic effects on a molecular level. In this study, 5 mg cm−2 of cream containing RA (tretinoin) (0.05%) (Aberela®, Jannsen-Cilag) was applied topically on each explant, using positive displacement pipette. 5 mg cm−2 is a typical finite dose applied on skin penetration test and considered as a non-occlusive dose 32,33.

The study was vehicle controlled with a cream containing all ingredients except for RA, referred to as placebo. The placebo cream was applied in the same conditions as the test cream.

The explants were incubated with the RA cream for 24 h at 37°C in 5% CO2 humidified air. In the end of the incubation period, two Ø 3 mm biopsies were taken from each explants following RNA extraction.

Skin explant viability

Alamar blue® (Invitrogen, Stockholm, Sweden) was used as reagent for measuring skin viability. Ø 3 mm biopsies were used for corresponding Alamar blue tests according to manufacturer's protocol. Briefly, Ø 3 mm biopsy was taken from each organoculture and put in 200 μL of Alamar blue dye and media at 1 : 10 ratio for 15 h at 37°C. Absorbance was measured at wavelengths of 570 nm and 600 nm to calculate Alamar Blue reduction, which is a measure of reduction reactants of cellular metabolism.

RNA isolation

Total RNAs were extracted from skin samples using QIAzol Lysis Reagent (Qiagen). Isolation begins with immersion of skin explants into QIAzol Lysis Reagent (Qiagen) followed by mechanical homogenization of the explants with a rotator–stator homogenizer (TissueRuptor, Qiagen) on ice for 20 s. RNA isolation was performed using the RNeasy mini kit (Qiagen), and DNase I digestion was performed following the manufacturer's instruction.

Microarray expression analysis

RNA concentration was measured with ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE), and RNA quality was evaluated using the Agilent 2100 Bioanalyzer system (Agilent Technologies Inc, Palo Alto, CA).

The 250 nanograms of total RNA from each sample were used to prepare biotinylated fragmented cRNA according to the GeneChip® 3′ IVT Express Kit Manual (Affymetrix Inc., Santa Clara, CA). Affymetrix GeneChip® expression arrays (Human Genome U133 Plus 2.0) were hybridized for 16 h in a 45°C incubator, rotated at 60 rpm. The arrays were washed and stained using the Fluidics Station 450 and finally scanned using the GeneChip® Scanner 3000 7G (Affymetrix).

Microarray data analysis

Subsequent analysis of the gene expression data was carried out in the freely available statistical computing language R (http://www.r-project.org) using packages available from the Bioconductor project (www.bioconductor.org) 34. The raw data were normalized using the robust multi-array average (RMA) method 35,36. To search for the differentially expressed genes between the X samples and the Y samples group, an empirical Bayes moderated t-test was then applied using the ‘limma’ package 37. To address the problem with multiple testing, the P-values were adjusted using the method of Benjamini and Hochberg 38.

Regulation of individual probe sets in gene ontology categories using gene-set enrichment analysis

The enrichment analysis 39 was performed to determine the probability that gene ontology groups are specifically overexpressed in RA-treated skin as compared to vehicle-treated skin. For this analysis, we used the public DAVID Bioinformatics database; http://david.abcc.ncifcrf.gov/ 40. The enrichment score of each gene was determined by calculating the relative expression in the RA-treated skin compared with placebo-treated skin.

Genes that were significantly upregulated or downregulated more than twofold were included in this analysis to evaluate enriched functional gene ontology groups.

To determine the overrepresented categories, a GO category-specific list for each category of interest was compiled containing those probe sets identified as differentially expressed. A GO enrichment score of >1 was set as a limit to be discussed as a finding.

As a further step, genes that were more than 10-fold differentially expressed were analysed manually to investigate which genes were affected most predominantly by RA.

Real-time reverse-transcriptase polymerase chain reaction (RT-PCR)

To validate the mRNA microarray results, three significantly upregulated genes were selected and their expression was analysed with real-time quantitative polymerase chain reaction (RT-qPCR) (Bio-Rad, Stockholm, Sweden). Skin explants from facial or breast skin from nine patients with a median age of 37 years were used for this conformational study. It is noteworthy that these donors differ from the donors used for the micro-array study. Total RNA from these skin explants treated with RA, as described above, were prepared with RNeasy mini kit (Qiagen) and treated with DNase I (Qiagen) to remove contaminating genomic DNA. Reverse transcription was performed with iScript cDNA synthesis kit (Bio-Rad). SsoFast EvaGreen supermix (Bio-Rad) was used in the RT-qPCR together with specific primers (from Qiagen) for the genes of interest, that is, DHRS9, CYP26B1 and HBEGF. Hot-start two-step PCR, consisting of 95°C for 3 min, followed by 40 cycles of 95°C for 10 s, 57°C for 30 s was performed using CFX96 Real-Time PCR Detection System (Bio-Rad). Relative quantification of gene expression in the different explants was calculated using the delta-delta-CT method using the TATA box-binding protein (TBP) as a reference gene.

Results

Alamar blue assay

Explants from the three donors tested for the gene array were also tested for viability using Alamar blue assay. The viability of explants was investigated after topical application of RA for 24 h. The results showed that RA was not affecting the viability of the explants during a 24-h application (Fig.1).

Figure 1

RA does not affect the viability of skin explants. Explants were exposed for 24 h to topical RA treatment (Aberela®) at a concentration of 0.05%. Data are expressed as the mean ± SEM of three independent experiments on three explants from different donors. Two-way ANOVA test was used for statistical analysis.

Identification of RA-responsive genes by 24 h stimulation ex vivo (Step 1)

We compared the transcriptional profiles of RA-treated and placebo-treated skin explants (n = 3) using Affymetrix oligonucleotide microarrays. The microarrays simultaneously measure the levels of 22 000 mRNAs. Among them were 93 significantly regulated by RA, within the limits of our selection criteria (less than twofold differentially expressed, P < 0.05). Of these 93 genes, 60 genes were upregulated and 33 downregulated.

To understand the biological meaning behind such a large list of genes, functional annotation clustering analysis was performed using DAVID Bioinformatics resources. The clusters with the smallest P-value (P < 0.05 adjusted using the method of Benjamini and Hochberg 38) which were retrieved after 24-h stimulation with RA are presented in Fig.2. When we summarized the functional categories of the genes upregulated by RA, we found that more than 30% of them fall into one of the following ontological categories: clusters (lipid, hormone and retinoid metabolism) and in development (organ, tissue and epithelial development) (Table1).

Figure 2

Heat map of hierarchical clustering of 93 genes in RA-treated explants compared with placebo. Hierarchical clustering was performed using the genesis software 67, with default settings (e.g. Euclidian distance, average linkage) of genes differentially expressed more than less than twofold, that is, 93 genes. Green indicates reduced expression, black indicates the unaltered expression, and red indicates increased expression in RA-treated as compared to placebo-treated explants (n = 3). The colour scale bar is shown at the top of each figure.Cluster analyses of altered genes involved in metabolism genes and development genes. Hierarchal clustering analysis was performed in both the gene (row) and experiment (column) dimension. Yellow arrows indicate genes enriched in gene annotation clustering (David Bioinformatics) which resulted in cluster 1 representing genes involved in metabolism (Enrichment score 3.6). 5 genes of 10 genes in the metabolism cluster were defining genes involved in metabolism of RA. Blue arrows indicate genes enriched in cluster 2 which were defined as involved in development (Enrichment score 2.0) by David Bioinformatics.

Table 1

Summary of the corresponding genes to the clusters within Fig.2

Gene name	Gene symbol	Gene ID	Fold change
Cluster 1a: Metabolism
Cytochrome P450, family 26, subfamily A, polypeptide 1	CYP26A1	206424_at	22.61
Dehydrogenase/reductase (SDR family) member 9	DHRS9	224009_×_at	11.96
Cytochrome P450, family 26, subfamily B, polypeptide 1	CYP26B1	219825_at	5.26
Acyl-CoA synthetase long-chain family member 3	ACSL3	201662_s_at	3.66
Proprotein convertase subtilisin/kexin type 5	PCSK5	213652_at	2.69
Nuclear receptor subfamily 3, group C, member 1	NR3C1	201866_s_at	2.36
Mitogen-activated protein kinase 14	MAPK14	211561_×_at	−2.13
Retinol dehydrogenase 12 (all-trans/9-cis/11-cis)	RDH12	242998_at	−2.42
Diacylglycerol O-acyltransferase homologue 2 (mouse)	DGAT2	226064_s_at	−4.58
Arachidonate 12-lipoxygenase	ALOX12	207206_s_at	−4.73
Serine carboxypeptidase 1	SCPEP1	218217_at	−6.41
Cluster 1 b; Retinoic acid metabolism
Cytochrome P450, family 26, subfamily A, polypeptide 1	CYP26A1	206424_at	22.61
Dehydrogenase/reductase (SDR family) member 9	DHRS9	224009_×_at	11.96
Cytochrome P450, family 26, subfamily B, polypeptide 1	CYP26B1	219825_at	5.26
Retinol dehydrogenase 12 (all-trans/9-cis/11-cis)	RDH12	242998_at	−2.42
Serine carboxypeptidase 1	SCPEP1	218217_at	−6.41
Cluster 2: Development
Keratin 4	KRT4	213240_s_at	26.97
Cytochrome P450, family 26, subfamily A, polypeptide 1	CYP26A1	206424_at	22.61
Dehydrogenase/reductase (SDR family) member 9	DHRS9	224009_×_at	11.96
Epithelial membrane protein 1	EMP1	213895_at	7.31
F-box protein 32	FBXO32	241762_at	6.39
Paired-like homeodomain 1	PITX1	209587_at	4.54
Sciellin	SCEL	1554921_a_at	4.06
Nuclear-receptor-interacting protein 1	NRIP1	202600_s_at	3.70
Heparin-binding EGF-like growth factor	HBEGF	38037_at	10.8
Ephrin-B2	EFNB2	202669_s_at	3.18
Proprotein convertase subtilisin/kexin type 5	PCSK5	213652_at	2.69
Nuclear receptor subfamily 3, group C, member 1 (glucocorticoid receptor)	NR3C1	201865_×_at	2.11
Sprouty homologue 1, antagonist of FGF signalling (Drosophila)	SPRY1	212558_at	−2.12
Mitogen-activated protein kinase 14	MAPK14	211561_×_at	−2.13
Dystroglycan 1 (dystrophin-associated glycoprotein 1)	DAG1	205417_s_at	−2.28
Crystallin, alpha B	CRYAB	209283_at	−2.66
Prolactin receptor	PRLR	227629_at	−3.28
Nebulin	NEB	205054_at	−4.53
Chromosome 1 open-reading frame 68	C1orf68	217087_at	−5.92

Table includes gene name, gene symbol, gene ID and fold change compared with placebo. Explain more the redundancy.

Cluster 1: Metabolic process (including RA metabolic genes) (enrichment score 3.6)

As shown by the functional annotation clustering analysis, metabolic process was the top cluster identified in our study (enrichment score 3.6) highlighting RA metabolism pathways with totally 11 genes included (Table1). Two main enzymes were identified, involved in the metabolism of RA metabolism, that is, cytochrome P450 26 A1 (CYP26A1), which was highly upregulated and to a lesser extend, and CYP26B1.

Also, the SDR family member dehydrogenase reductase 9 (DHRS9), acyl-CoA synthetase long-chain family member 3 (ACSL3), proprotein convertase subtilisin/kexin type 5 (PCSK5), and nuclear receptor subfamily 3 (NR3C1) was found to be upregulated in our study.

Genes that shown to be downregulated in this cluster were the inflammatory gene arachidonate 12-lipoxygenase (ALOX12), diacylglycerol O-acyltransferase (DGAT2), mitogen-activated protein kinase 14 (MAPK14), serine carboxypeptidase 1 (SCPEP1) and retinol dehydrogenase 12 (RDH12).

Cluster 2: Development genes (enrichment score 2.0)

The second cluster was identified as the development genes including tissue development and epithelial development (Table1 and Fig.2). In total, 19 genes were identified in this cluster. Keratin 4 (KRT4) was found to be highly upregulated. Within the same cluster, heparin-binding EGF-like growth factor (HBEGF), CYP26A1, DHRS9, epithelial membrane protein 1 (EMP1), F-box protein 32 (FBX032), paired-like homeodomain 1 (PITX1), sciellin (SCEL), nuclear-receptor-interacting protein 1 (NRIP1), ephrin-B2 (EFNB2), proprotein convertase subtilisin/kexin type 5 (PCSK5) and NR3C1 were identified as genes induced by RA.

Downregulated genes within the development cluster were as follows: sprouty homologue 1 (SPRY1), MAPK14, dystroglycan (DAG1), crystallin, alpha B (CRYAB), prolactin receptor (PRLR), nebulin (NEB), and chromosome 1 open-reading frame 68 (C1orf68), CYP26A1, DHRS9, NR3C1 and MAPK14 were identified in both cluster 1 and 2.

Validation of microarray data with qRT-PCR

To confirm independently the microarray results obtained, we performed quantitative RT-PCR analysis of three representative genes, that is, CYP26B1, DHRS9 and HBEGF (Fig.3). These three genes were chosen based on following criteria; they were significantly upregulated by RA and included in the gene ontology groups that were detected via functional annotation clustering using DAVID Bioinformatics resources. These genes are also known to be regulated by RA and expressed in the two top 2 clusters of genes in our explant study: vitamin A metabolism (CYP26B1, DHRS9) and development (HBEGF). Results showed upregulation of CYP26B1 (1.52-fold), HBEGF (3.22-fold) and DHRS9 (3.53-fold) by topical application of RA.

Figure 3

Validation of microarray findings by qPCR analysis. qPCR confirmed upregulation of CYP26B1 (1.52-fold), HBEGF (3.22-fold) and DHRS9 (3.53-fold) by RA in explants from different donors. Wilcoxon signed-rank test used for statistical analysis. Data are presented as relative concentration of gene expression (mean + SD, n = 9).

Step 2: Analysis of genes affected more than 10-fold by RA

Upon analysing the gene array manually for genes affected more than 10-fold, RA was shown to affect 8 genes more than 10-fold (Table2).

Table 2

The 8 genes that were stimulated more than 10-fold by topical application of RA

Gene title	Gene symbol	Probe Set ID	Fold change
Flavin-containing monooxygenase 2 (non-functional)	FMO2	228268_at	50.17
Sterile alpha motif domain containing 9	SAMD9	228531_at	40.02
Carcinoembryonic antigen-related cell adhesion molecule 6	CEACAM6	211657_at	33.64
Keratin 4	KRT4	213240_s_at	26.97
Cytochrome P450, family 26, subfamily A	CYP26A1	206424_at	22.61
Chromosome 10 open-reading frame 54	C10orf54	225372_at	14.42
Dehydrogenase/reductase (SDR family) member 9	DHRS9	224009_x_at	11.96
Lipocalin 2	LCN2	212531_at	10.73

P < 0.05.

Flavin-containing monooxygenase (FMO) was in our study stimulated 50-fold by topical application of RA. Next, sterile alpha motif domain containing 9 (SAMD9) was upregulated 40-fold in our study.

Furthermore, carcinoembryonic antigen-related cell adhesion molecule 6 (CEACAM6/cd66c) was induced 34-fold.

Keratin 4 (KRT4) was also highly upregulated (27-fold) and was linked to the cluster of development genes by analysis using DAVID Bioinformatics. Same accounts for cytochrome P450, family 26, subfamily A (CYP26A) which was already identified in the metabolism cluster.

Chromosome 10 open-reading frame 54 (c10orf54) was upregulated 14-fold. This gene encodes for platelet receptor Gi24. As already described, dehydrogenase/reductase (SDR family) member 9 (DHRS9) was highly upregulated (12-fold). Finally, lipocalin 2 (LCN2) was upregulated 11-fold.

Discussion

Several previous reports describe the genetic targets of retinoids using gene array techniques 41. However, the genetic profile of RA has not earlier been studied in explants of full-thickness skin. In this study, we developed an ex vivo skin explant organoculture model to study the complete gene expression following RA treatment, adopting cDNA micro-array analysis.

Concomitantly, we described an approach for recovering RNA of full-thickness skin taken from the explant model skin. This procedure permits a direct quantitative and qualitative assessment of topical treatment with dermatological active ingredients, such as RA.

We found that 93 genes were significantly more than two-fold regulated by RA in three organocultures of human skin. The marked alterations of the expressions of 93 genes were found to be involved in mainly two clusters, that is, lipid and retinoid metabolism and in development, that is, organ, tissue and epithelial development. Comparing to previous studies in monocultures of keratinocytes, Lee et al. showed that these cells responded to RA very fast with 315 genes affected already at 1 h. Please note that these authors used a threshold of 1.62-fold change to show effect of RA 41.

Lipid and retinoid metabolism

As shown by the functional annotation clustering analysis, metabolic process was the top cluster identified in our study, highlighting lipid and RA metabolism pathways. Two main enzymes were identified with regard to retinoic acid metabolism, that is, CYP26A1, which was highly upregulated (22.6-fold) and to a lesser extend, and CYP26B1 (5.3). The importance of these enzymes with regard to retinoid metabolism in the human skin has been confirmed by other authors 22,42. For example, Lee et al. showed induction of CYP26B1 (12.8-fold) by RA in epidermal keratinocytes after 4-h stimulation 41. In contrast to our study, these authors have showed that CYP26B1 is the most abundant RA metabolizing CYP26 enzyme in unstimulated keratinocytes at the mRNA level 42,22,43. On the other hand, they suggested that CYP26A1 may play a role during situations of excess of RA in epidermis which is the situation in our study.

The expression of DHRS9 was also increased in our study (11.96-fold) and is in agreement with Lee et al.'s study on epidermal keratinocytes where DHRS9 was induced 3.5-fold after 24 h 41.

Other genes affected involved in metabolism were ACSL3 which was upregulated in our study which is concomitant with reports on the effect of RA on rat hepatoma cells 44. However, to our knowledge, such an effect of RA has not been reported in human skin to date.

The lipoxygenase, ALOX12, involved in fatty acid oxidation, was reduced after RA treatment. These results are compatible with Lee et al.'s study in human keratinocytes where ALOX12B was reduced after RA treatment 41.

Moreover, we found DGAT2 to be reduced (4.58-fold). Same accounts for the expression of RDH12, which was also reduced (2.42-fold) in response to exogenous RA in our study. To our knowledge, RDH12 has not been described in human skin previously.

In addition, MAP kinase 14, which is essential for embryonic development and regulates different cellular functions such as proliferation, differentiation, cell death, adhesion, migration, as well as the response to stress, and many metabolic pathways were found to be downregulated in our study. A similar enzyme, MAPK13, was downregulated in human keratinocytes by RA treatment. The inhibitory effect of RA on MAPK is documented by other authors and has been hypothesized to have a role in the inhibitory action of RA on UV-induced skin damage 45.

Regulation of development genes by RA

Retinoid signalling is known to be important in the development of the eye 46, lungs 47, heart 48 and skin 49,50. Adult skin continues to grow throughout the life of humans, and therefore, it is not surprising that genes involved in epidermal development are stimulated after RA application.

In our study, development genes were identified as the second cluster of importance by gene ontology analysis. In this cluster, KRT4 was found to be highly upregulated which is also well known in human keratinocytes 41. It has even been suggested as a future marker for retinoid bioactivity in epidermis in vivo 41. Within the same cluster, HBEGF was increased which also is consistent with literature data on the proliferative effect of RA in human keratinocytes 41 and in reconstituted human epidermis 30.

HBEGF was confirmed to be regulated by RA by qPCR in skin explants (Fig.3). The discrepancy in fold change (10.8 in the microarray study and 3.22 by qPCR) could be due to the difference in donors as well as skin site (facial and breast). In fact, Cyp26a1, Cyp26b1, Rbp1, Crabp1, Hbegf and Krt4 are all known marker for retinoid activity (Gericke et al i.e. ref 65).

Interestingly, EMP1 was found to be highly upregulated (7.31-fold), which is consistent with reports in F9 teratocarcinoma stem cells 51. EMP1 has been prominently found in the gastrointestinal tract, skin, lung and brain but not in liver. The proteins of this family are likely to serve similar functions possibly related to cell proliferation and differentiation in a variety of cell types 52.

Another gene in the development cluster, sciellin (SCEL), a precursor of the cornified envelope 53, was highly upregulated. The specific function of this gene in human skin is not known to date, but deletion of this gene in mice did not affect normal barrier function 54.

NRIP1 was also highlighted in the gene ontology analysis for development of genes, and interestingly, this biomarker has not been found in skin to date. However, it has been described in context of regulation of mesenchymal stem cell differentiation by retinoids as NRIP1 is induced by RA in stem cells 55.

The upregulation of EfnB2 seen by retinoic acid in our study could implicate the regulation angiogenesis and blood vessel development as this has been reported in the literature 56. However, the role of EfnB2 has not been investigated to date in human skin.

It is noteworthy that NR3C1 also known as glucocorticoid receptor is upregulated in human skin by topical application of RA in our study. To our knowledge, this mechanism has not been reported in human skin. However, a study in cultured chick embryonic skin showed that the addition of 2 μM retinol resulted in a threefold increase in specific dexamethasone binding to the glucocorticoid receptor 57. Whether the same mechanism applies to our model remains to be elucidated.

Prolactin receptor (PRLR) was profoundly downregulated which is in line with a study in patients with acne where isotretinoin caused suppression of several pituitary hormones such as prolactin in patients' sera 58.

Upon analysing the gene array manually for genes affected more than 10-fold by RA without RA was found to upregulate 8 genes more than 10-fold. One of these was FMO, which was upregulated 50.17-fold after 24-h stimulation with RA. FMO has been found to be present in human skin 59. It is known to metabolize xenobiotics, but to our knowledge, not known to metabolize retinoids.

SAMD9 was upregulated 40.02-fold by topical application of RA for 24 h. The function of SAMD9 in human skin is to date unknown, but it has been suggested to have a role in regulation of cell death in malignant glioma 60.

Furthermore, CEACAM6/cd66c was induced 34.64-fold. CEACAM6 is found on epithelial and other cells. It has to our knowledge not been described to be expressed in human skin. Other genes from the same family, for example CEACAM1, have shown to be important for the wound healing process in human skin 61.

KRT4 was also highly upregulated (26.97-fold) and was linked to the cluster of development genes by analysis using DAVID Bioinformatics. Same accounts for cytochrome P450, family 26, subfamily A (CYP26A), which was already identified in the metabolism cluster. C10orf54 was upregulated 14.42-fold. This gene encodes for platelet receptor Gi24. The function of this receptor in human skin is to our knowledge unknown. As already described, DHRS9 was highly upregulated (11.96-fold). LCN2 was upregulated 10.73-fold. The induction of LCN by RA has also been confirmed by other authors. The confirming studies were investigating the effect of oral application of isotretinoin in patients with acne 62,63 but also of RA in monocultures of keratinocytes 41.

In conclusion, we have in the current study shown that a full-thickness organoculture of human skin is an efficient non-invasive model to study the effects of topically applied compounds at the gene expression level.

We have in this study confirmed several genes affected by RA, which has been previously found by other authors in monocultures of keratinocytes 41, reconstructed epidermis 30 and in skin biopsies after oral application of isotretinoin 64, for example, CYP26A1/B1, DHRS9, ALOX12, KRT4, HBEGF and LCN2. In fact, KRT4, HBEGF and CYP26A1/B1 are considered to be retinoid target genes 65 (Gericke et al., 2013) where KRT 4 is a marker for keratinocyte differentiation 66 (Wanner et al. 1996), CYP26 for retinoid metabolism and HBEGF for increased keratinocyte proliferation and epidermal thickness 8.

We have also disclosed several genes that have never been identified in human skin before, for example SCEL, EfnB2, SAMD9 and CEACAM6. Consequently, the roles of these genes in skin have yet to be discovered.

Overall, we believe that this organoculture model possessing a physiological complexity such as a fully functioning barrier function compared with cell cultures and reconstructed skin models could serve as a model to evaluate the effect of topical application of cosmetic ingredients. The model has been shown in this study to be suitable for assessing short-term stimulation with RA. Therefore, this model can be used as a first screen to identify more effective therapeutic strategies by testing functional connections between drugs and genes using databases such as connectivity map (Cmap, Broad Institute). However, further evaluation will be needed including a larger sample group to confirm this hypothesis.