Vesicular hand eczema transcriptome analysis provides insights into its pathophysiology.

Hand eczema is a common inflammatory skin condition of the hands whose pathogenesis is largely unknown. More insight and knowledge of the disease on a more fundamental level might lead to a better understanding of the biological processes involved, which could provide possible new treatment strategies. We aimed to profile the transcriptome of lesional palmar epidermal skin of patients suffering from vesicular hand eczema using RNA-sequencing. RNA-sequencing was performed to identify differentially expressed genes in lesional vs. non-lesional palmar epidermal skin from a group of patients with vesicular hand eczema compared to healthy controls. Comprehensive real-time quantitative PCR analyses and immunohistochemistry were used for validation of candidate genes and protein profiles for vesicular hand eczema. Overall, a significant and high expression of genes/proteins involved in keratinocyte host defense and inflammation was found in lesional skin. Furthermore, we detected several molecules, both up or downregulated in lesional skin, which are involved in epidermal differentiation. Immune signalling genes were found to be upregulated in lesional skin, albeit with relatively low expression levels. Non-lesional patient skin showed no significant differences compared to healthy control skin. Lesional vesicular hand eczema skin shows a distinct expression profile compared to non-lesional skin and healthy control skin. Notably, the overall results indicate a large overlap between vesicular hand eczema and earlier reported atopic dermatitis lesional transcriptome profiles, which suggests that treatments for atopic dermatitis could also be effective in (vesicular) hand eczema.

INTRODUCTION

Hand eczema (HE) is a common skin disease, with a 1‐year prevalence of up to 10% of the general population. Since HE has a big socioeconomic impact, , so far, studies have focused on the treatment of this disease. , , , However, hitherto limited research concentrated on the pathogenesis of HE. Publications are mainly restricted to epidemiological association studies and genetic investigations closely linked to atopic dermatitis (AD). For example, it was shown that mutations in the filaggrin gene (FLG) predict a persistent form of HE in patients with AD, , and others reported no association between FLG mutations and HE in individuals without AD. ,

Hand eczema is a multifactorial skin disorder in which both endogenous and exogenous factors seem to play a role. Based on aetiology, HE can be classified into allergic contact dermatitis, irritant contact dermatitis, atopic HE, unclassified HE, protein contact dermatitis (PCD) and a combination of these. Morphologic classification of HE includes pulpitis, recurrent vesicular HE (VHE), hyperkeratotic HE, chronic fissured HE and nummular HE. ,  VHE, which can be present without any known aetiological factors, is morphologically well characterized by frequent eruptions of vesicles on the palms, palmar or lateral aspects of the digits. , Its episodic nature grants the opportunity to follow the development of disease morphology.

In this study, we performed RNA‐sequencing (RNA‐seq) to determine the transcriptome of vesicular lesional epidermis and non‐lesional epidermis of VHE patients compared to epidermis from healthy controls. For a large set of either biologically relevant or highly up and downregulated genes according to RNA‐seq, we validated mRNA expression by real time quantitative PCR (RT‐qPCR). Additionally, for a subset of these candidate genes we studied protein level expression using immunohistochemistry (IHC). We here present the first transcriptome analysis of VHE lesion skin, which reveals a large overlap with the AD transcriptome based on previously reported studies. Non‐lesional VHE skin, however, is not significantly different from healthy controls as is the case with non‐lesional AD skin.

METHODS

Study design and participants

Skin biopsies were collected from 10 adult patients diagnosed with VHE, as defined by the Danish Contact Dermatitis Group and current guidelines. , Exclusion criteria were usage of topical corticosteroids on the hands 2 days before taking biopsies; systemic treatment in the last 2 weeks; ultraviolet radiation therapy in the last 4 weeks; active bacterial, fungal or viral infection of the hands; other skin diseases of the hands; proven contact sensitization with clinical relevance to the hands, in which exposure to allergens is not avoided. To compare lesional (L) and non‐lesional (NL) VHE skin samples with controls (C), biopsies from normal palmar skin were collected from 10 healthy adult subjects. Furthermore, patient and disease characteristics were collected. All procedures performed in this study were in accordance with the ethical standards of the institutional research committee (Medical Ethical Review Board of the University Medical Center Groningen, reference: METc 2018/018) and with the Declaration of Helsinki and its later amendments or comparable ethical standards. Informed consent was obtained from all participants included in this study, including usage of photographs of the hands.

Biopsies

Biopsy locations were marked and photographed in advance. Two punch biopsies of 5 mm were taken from patients with VHE. The first biopsy was taken from lesional skin with visible vesicles, preferably with minimal other signs of inflammation. The second biopsy was taken from NL skin of the hands, defined as palmar skin without any vesicles and none or minimal other signs (e.g. erythema, squamae, infiltration). Both biopsies were taken with a minimum distance of 1 cm from each other. From healthy controls (C), the biopsy was taken from the hypothenar region of the non‐dominant hand. Biopsies were cut in half; one half was formalin‐fixed and embedded in paraffin for histological analysis. The epidermis of the second half of the biopsies was separated from the dermis by dispase (Roche) treatment for 24 h at 4°C as previously described, and the epidermis was subsequently frozen at −80°C until further use (for RNA isolation).

RNA isolation, RNA‐sequencing and data analysis

RNA isolation was performed as described earlier. RNA quality control, Illumina TruSeq strand‐specific mRNA polyA library preparation, whole transcriptome sequencing, sample demultiplexing and basic data analysis was performed by BaseClear B.V. (Leiden, The Netherlands). Sequencing was performed on an Illumina NovaSeq system with paired‐end 150 cycles protocol and indexing. RNA‐sequencing reads were adapter trimmed and quality‐checked using Trim Galore version 0.6.2.  Next, RNA‐seq reads were mapped to the human reference genome using STAR aligner version 2.7.1.  The aligner identifies splice junctions between exons based on the mapping results. Assignment of mapped reads to genomic features (counting) was performed with featureCounts version 1.6.3.  The normalized counts per sample were calculated using the TPM (Transcripts Per Kilobase Million) method that takes into account the gene length and the sequencing library size.

Real‐time quantitative PCR

RT‐qPCR analysis was performed as described earlier. All primers were designed and used as described previously.  Target gene expression was normalized to the expression of the house keeping gene human acidic ribosomal phosphoprotein P0 (RPLP0). The ΔΔCt method was used to calculate relative mRNA expression levels.  The primers used in this study are listed in Table S1.

Morphological and immunohistochemical analysis

Tissue sections of 6 µm of hand skin biopsies were stained with hematoxylin and eosin (Sigma‐Aldrich, St. Louis, MO) or processed for IHC analysis. Sections were stained with antibodies using an indirect immunoperoxidase technique (Vectastain, Vector Laboratories, Burlingame, CA) as previously described. An overview of all used antibodies in this study is presented in Table S2.

Quantification of protein expression in IHC biopsy sections

Image acquisition of stained IHC biopsy sections was performed by a ZEISS Axiocam 105 colour Digital Camera and a 10× or 20× objective. The ZEISS Axiocam 105 colour is a compact 5 megapixel camera (2560 × 1920 pixels) for high‐resolution images with a 1/2.5” sensor. The images were stored in CZI format. The images were analysed with the cell image analysis software CellProfiler (Broad Institute). In CellProfiler different algorithms for image analysis are available as individual modules that can be modified and placed in sequential order to form a pipeline that can be used to identify and measure biological objects and features in images. Pipelines for Ki‐67, loricrin, filaggrin and pan‐LCE3 analysis were created (available on request). Data visualization and statistical analysis were performed using Instant Clue Software.

Statistical analysis

Differential gene expression analysis was performed with DESeq2, with a two‐factor design for paired samples (i.e. patient L vs. patient NL) and with a single‐factor design for all other comparisons (i.e. healthy control vs. patient). The p‐values for statistical differences between the contrasts were Benjamini‐Hochberg multiple‐testing corrected (FDR).

The fold‐change in gene expression is expressed as log2 ratio in average expression for a gene between the sample groups. Quantitative PCR ΔCt values were analysed using one‐way ANOVA Kruskal‐Wallis analysis followed by Dunn's post‐hoc analysis to test for significant differences between L, NL and C skin. For quantification of protein expression data a Kruskal‐Wallis test was conducted to examine the differences in cell counts and protein expression. For all tests, p < 0.05 was considered statistically significant.

RESULTS

Descriptive characteristics of the study volunteers

The study included 10 VHE patients and 10 healthy controls. All participants were Caucasian and the healthy controls were sex‐matched with the patient group, 7 females and 3 males. The mean age was 48.8 years (range 21–65 years) and 55.2 years (range 40–64 years) for the patient group and controls, respectively. Patients had a mean age of onset of 42.4 years (range 16–62 years) with a mean disease duration of 6.4 years (range 1–24 years). Patient characteristics are shown in Table 1, and information of the control group can be found in Table S3, respectively. All patients showed moderate to very severe vesicular eruptions on the palms and palmar or lateral sides of the fingers, combined with erythema.

TABLE 1

Patient characteristics

Cases	Age	Sex	Disease duration, years	Atopy a	Atopic comorbidities	Age of onset AD	Contact sensitization b	Severity (Photoguide)	Contributing aetiological factors
1	26	F	1	‐	‐	‐	‐	Moderate	‐
2	52	M	2	+	‐	‐	+	Very severe	‐
3	56	F	6	+	pAD, AR/AC, asthma c	17	+	Moderate	‐
4	21	M	5	+	‐	‐	+	Severe	‐
5	64	F	24	+	pAD, AR/AC, asthma c	50	+	Severe	ICD
6	29	M	1	+	AR/AC	‐	‐	Severe	‐
7	65	F	3	‐	‐	‐	+	Moderate	ICD
8	46	F	4	+	AR/AC	‐	+	Severe	‐
9	65	F	15	+	‐	‐	+	Very severe	ICD
10	64	F	3	‐	‐	‐	+	Very severe	‐

Abbreviations: AR/AC, allergic rhinitis/conjunctivitis, based on the ARIA criteria ; F, female; ICD, irritant contact dermatitis; M, male; pAD, previous atopic dermatitis, based on the UK criteria.

Atopy based on specific IgE inhalant allergens >0.99 kU/L.

≥1 positive patch test reaction to allergens from the European Baseline Series.  None of the patients had current exposure to these allergens.

Based on the GINA guideline.

RNA‐sequencing data analysis shows distinct gene expression profiles for lesional compared to non‐lesional skin and healthy control skin

Quantity and quality of epidermal RNA was checked by the sequence provider (data not shown), and 1 µg of each sample was used for library preparation for RNA‐seq. For a detailed overview of all normalized gene counts per study sample see File S1. RNA‐seq analysis yielded between 9.8 and 16.9 million readings obtained per sample (12.4 million on average). A large part of the reads (92.7%) was unique and could be mapped to the human genome (see RNA‐seq data statistics in File S2). Multivariate gene expression analysis by principal component analysis (PCA) shows a distinct expression profile for lesional VHE samples (Figure 1). This PCA reflects the distances between study samples (L, NL and C) based on their individual gene expression profiles for all genes measured in the total dataset, showing that L samples are separated from the NL and C samples. In total, 28 614 genes were analysed of which 52 were differentially expressed according to the predetermined stringent criteria (for a complete overview of differential expressed genes see the volcano plots at Figure S1). These criteria can be adjusted in the File S1 to alter the list of differentially expressed candidate genes. Of these 52 genes, 41 are upregulated, while 11 are downregulated. In contrast to lesional skin, a large overlap was observed between C and NL transcriptomes (Figure 2). Lesional skin of VHE patients showed significant upregulation of several genes involved in keratinocyte host defense and inflammation (e.g. LTF, LYZ, LCN2, LCE3A, PI3, and genes belonging to the S100 gene family: S100A7A, S100A7, S100A8 and S100A9), epidermal proliferation and differentiation (e.g. SPRR2A, SPRR2B, SPRR2D, SERPINA3, SERPINB3, SERPINB4, and keratin family members KRT6, KRT16 and KRT17), and immune signalling (MMP12, CHI3L2, CCL22, IL4R). A relative small number of genes was significantly downregulated in lesional skin, most notably LOR (epidermal differentiation), IL37 (immune signalling), and LCE1D/E (host defense). S100A7A was the most differentially expressed gene (fold‐changes of 13.98 and 11.14 for L vs. NL and C, respectively). S100A7, S100A8 and S100A9 were also represented at the top of this list of most upregulated genes in lesional skin (Figure 2). In addition, we found a set of upregulated genes in lesional skin with a relationship to HE and AD (e.g. HAS3, CDH3, TNC, DPP4, C1R, C1S, CD207, TMEM173, CD1A), and also a number of upregulated genes that have not yet been described to be expressed in skin or related to HE or AD (e.g. HEPHL1, BIRC3, PRSS53, LCP1, LAPTM5, TYMP, PARP9). In deeper downstream bioinformatics analysis of the transcriptomics data, Gene Set Enrichment Analysis (GSEA) based on GO terms hinted towards the involvement of adaptive immune system processes in lesional VHE skin, as a strong and significant enrichment of upregulated genes belonging to this GO:0002250 gene set was detected in comparison to matching non‐lesional VHE skin samples (normalized enrichment score of 3.0 for NL to L, Gene Ratio of 0.5, FDR p = 0.0025, see Figure S2, and File S3 for more details).

FIGURE 1

Multivariate gene expression analysis by PCA shows a distinct expression profile for lesional VHE samples. This representation of a principal component analysis (PCA) reflects the distances between study samples based on the expression (as normalized counts) of all genes for each sample (n = 28 614 genes). The black arrows point to the top 20 individual genes (names shown in italics) which best explain compositional differences between the samples. The genes in the left quadrants associate with C and NL samples (with in green as in correspondence with univariate analysis, see Figure 2), genes in the right quadrants associate with L samples (idem, in red). The PCA was made using Canoco 5.12  software using default settings of the analysis type ‘Unconstrained’, with log‐transformation set to 1 and centering of gene expression data. C, control (green circles); NL, non‐lesional (orange squares); L, lesional (red diamonds)

FIGURE 2

The most up/downregulated genes in VHE according to RNA‐seq. The top up and downregulated genes according to RNA‐seq when applying a threshold of threefold differences in gene expression between study samples (C, NL and L). Fold‐change is expressed as the log2 ratio in average expression for a gene between the study sample groups, and gene expression is expressed as the average normalized counts for a sample group. A gene is listed here only if it is significantly differentially expressed for at least one of the three study contrasts (according to FDR‐corrected p‐values as calculated by DESeq2), and if it meets the threshold of >500 averaged normalized counts between all sample groups. Note that fold‐change numbers in red are not significant. C, control; NL, non‐lesional; L, lesional

Real‐time quantitative PCR for validation of RNA‐seq data

To validate our RNA‐seq results, we selected a number of genes that were found to be upregulated in lesional skin for replication with RT‐qPCR. We also included genes that were reported to be up or downregulated in the existing HE literature, for example FLG, FLG‐2, HRNR, KLK5, KLK7, CST6, KRT10, CLDN1, MKI67, and a set of cytokeratins. , , In addition we added three immune signalling genes that are known to be involved in AD (CCL17, TSLP and IL32 ). For description of all genes see Table S1. The same RNA samples included in the RNA‐seq dataset were used for validation by RT‐qPCR. For 13 candidate genes (including the household gene RPLP0), RNA yield was sufficient to perform validation on the complete set of samples (n = 10), while for the remaining set of 25 genes we used less samples in the C and NL group (n = 6) as shown in Table S4A–B.

RT‐qPCR analyses confirmed that most of the investigated genes involved in epidermal differentiation were significantly upregulated (KRT6A, KRT14, KRT16, KRT17, SPRR2D, SERPINA3) in lesional VHE skin, in correspondence to the RNA‐seq data (Figure 3A). Four epidermal differentiation genes (LOR, KRT2, FLG, FLG‐2) were found to be significantly downregulated according to RT‐qPCR analysis (Figure 3A), as well as in our RNA‐seq analysis (fold‐changes are respectively −3.16, −2.97, −1.41, and −2.24) as shown in Figure 2 or calculated from File S1. We found no significant differential expression for MKI67, HRNR, KRT9, KRT10, CLDN1, and CST6, by RT‐qPCR analysis.

FIGURE 3

Gene expression of VHE skin. RT‐qPCR data represents fold differences in mRNA expression of NL and L skin compared to C (set to 1) of (A) epidermal proliferation and differentiation genes, (B) keratinocyte host defense genes and (C) immune signalling related genes. p‐values of one‐way ANOVA statistical analysis comparing all three contrasts are indicated with *p < 0.05, **p < 0.001, and ***p < 0.0001. Asterisks above L plots only represent significance of these contrasts compared to C (we found no significant differences between NL and C). Significance between NL and L are indicated by asterisks above the line. Data are represented as mean ± SD

The significant differential expression of keratinocyte host defense and inflammation genes in the RNA‐seq dataset (contrast: C vs. L) was also confirmed by RT‐qPCR analysis (Figure 3B). This set of genes is not only differentially expressed but is also the most abundant set of genes in lesional palmar skin as shown in the right columns of Figure 2 (gene expression) and in Table S4. The late cornified envelope genes (LCE3A, LCE1D and LCE1E), which were found to be significantly up or downregulated in the RNA‐seq dataset (with fold‐changes of 5.80, −3.38 and −2.91, respectively), could not be reproduced by RT‐qPCR.

The last group of genes that we validated concerns molecules involved in immune signalling: MMP12, CHI3L2, IL4R, and CCL22. RT‐qPCR analysis confirmed their significant upregulation in lesional skin (Figure 3C). IL32 and CCL17, two genes that were not included in the RNA‐seq data list of 52 differentially expressed genes (as they did not reach the threshold of >500 averaged normalized counts between the groups and, thus, are only very lowly expressed), were found to be significantly upregulated in lesional skin (Figure 3C). From File S1 fold‐changes of 5.17 and 6.99 were calculated for both genes. Finally, in line with the RNA‐seq results, IL37 was found to be downregulated in RT‐qPCR analysis, and the kallikreins KLK5 and KLK7 showed no significant differences according to both methods (both genes showed a slight downregulation in the RNA‐seq dataset with fold‐changes of −1.86 and −0.57).

Immunohistochemical analysis confirms transcriptome data

To confirm the gene expression data at the protein level, we performed IHC analyses of all samples in our cohort (n = 10, for each of each study group). An example of the clinical manifestation of lesional VHE skin and the location of the biopsy punch is shown in Figure 4A (the location of a biopsy in non‐lesional patient skin is shown in Figure S3). The clinically characteristic vesicles of lesional skin are clearly visible in the H&E staining (Figure 4B upper left). Large vesicles with infiltrate are present in the highly acanthotic epidermis. Epidermal thickness is significantly increased in L compared to C and NL skin, but no difference was observed between C and NL skin, as shown in Figure S4A. We also determined the targeted expression of specific proliferation, differentiation, host defense and immune signalling proteins (Figure 4B).

FIGURE 4

Epidermal morphology and protein expression of VHE lesional skin. (A) Clinical photograph with representative example of lesional skin, black circle indicates biopsy site. (B) H&E and immunohistochemical staining microscopic images of keratinocyte‐expressed structural proteins, antimicrobial proteins and immune cell markers (see Figure 3) in C and L biopsies. Bottom right black boxes show optical magnifications for each microscopic image. Images are representative for n, 10 individuals per group. C, control; L, lesional. Scale bar =100 µm

Proliferation of keratinocytes is similar for VHE lesional skin and healthy control skin (Figure S4B). The epidermal differentiation‐ and skin barrier‐associated proteins filaggrin, hornerin and loricrin show decreased expression levels in lesional skin of VHE patients (Figure 4B). This is also reflected by the disruption of the granular epidermal layer, which is likely affected by the vesicle formation within the epidermis. We further quantified IHC images for protein expression of filaggrin, loricrin and LCE3, which were all three strongly reduced in the L group (Figure S4C–E). Moreover, we found that based on these three characteristic keratinocyte proteins, NL skin is significantly different from L skin, but closely resembles C skin (Figure S4C–E). Cystatin M/E protein was absent in all palm biopsies. We also stained for three keratin proteins as some keratins were differentially expressed (up or downregulated) on mRNA level as described above (see also Figure 3 and Table S4). Total amount of keratin 10 protein might be elevated in L skin, due to epidermal acanthosis, but its expression profile in NL is comparable to C skin (data not shown). In contrast, lesional skin shows strong reduction of keratin 2 protein compared to C skin, not restricted to the vesicular area. This nicely reflects its gene expression pattern we found in the transcriptome analyses (fold‐change =−2.97) and qPCR analysis (see Figure 3). Keratin 16 is clearly present in lesional skin only. The H&E staining showed infiltrating cells in the dermis, the epidermis and in the vesicles of lesional skin (Figure 4). We additionally stained CD45, a leukocyte marker, to reveal infiltration, inflammation and migration of T‐cells from the dermis up to the epidermis in lesional skin. In inflammatory skin diseases like AD, the level of antimicrobial proteins (AMP) is elevated in the skin lesions. Therefore, we determined the expression of psoriasin, S100‐A8, SKALP and LCE3 (Figure 4). For all of these AMP, except for LCE3, we found that expression was absent in C and NL skin and only induced in the lesional skin samples. The amount of protein was especially high for psoriasin and S100‐A8, and notably, most suprabasal keratinocytes expressed these two proteins. Total LCE3 antimicrobial protein per biopsy is expected to be greater in lesional skin, as the epidermis is substantially thicker, however quantification reveals that the percentage of LCE3 protein is lower in L compared to NL and C epidermis (Figure S4E). In contrast to elevated mRNA expression for LCE3A (Figures 2 and 3), total LCE3 protein (LCE3A‐E) seems to be reduced in lesional skin.

DISCUSSION

Here we report the first RNA‐seq analysis of the VHE transcriptome. Overall, high expression of genes/proteins involved in keratinocyte host defense and inflammation was found in lesional VHE skin, with statistically significantly differential expression compared to healthy skin or internal controls for many of these candidates. Furthermore, molecules involved in epidermal differentiation were specifically up or downregulated, while immune signalling genes were upregulated in lesional skin, albeit with relatively low expression levels.

Previously, proteome analysis of lesional palmar skin of patients with chronic HE and healthy controls revealed that several skin barrier proteins, most notably filaggrin, filaggrin‐2 and hornerin, were downregulated in diseased skin compared to unaffected healthy skin, suggesting that skin barrier dysfunction plays a role in the pathogenesis of chronic HE. In another study, using a targeted approach on six candidate skin barrier genes in patients diagnosed with chronic HE, normalization (upregulation) was found for the mRNA expression levels of FLG, LOR and KRT10 after treatment with alitretinoin, which is indicated for treatment of severe chronic HE. Although both studies are limited in size and only focus on a small predefined set of genes and proteins, the expression levels of the studied skin barrier molecules are in agreement with our data. In contrast to both studies, we could not confirm upregulation of TSLP in lesional skin, nor downregulation of the cysteine protease inhibitor CST6 and the serine proteases KLK5 and KLK7.  Moreover, we observed complete absence of cystatin M/E protein expression in the hand palm biopsies, which is remarkable as this protein is strongly expressed in the stratum granulosum of keratinocytes at other skin locations. Proliferation of keratinocytes was similar between groups, as demonstrated by the MKI67/Ki‐67 gene and protein expression levels. This is in contrast to a recent study, in which epidermal hyperproliferation and an increased number of Ki‐67 positive cells was found in hyperkeratotic HE.  This might be explained by the studied HE subtype that is characterized by hyperkeratosis and desquamation. To identify chronic HE phenotypes, the mRNA expression of a set of genes involved in skin barrier function was recently studied using qPCR. Genes involved in filaggrin degradation, natural moisturizing factor synthesis, and structural barrier genes were downregulated suggesting that skin barrier dysfunction is an important parameter to discriminate chronic HE patients. However, a comparison with our data is difficult, since Tauber et al. analysed lesional dorsal skin, whereas we studied lesional palmar skin. Furthermore, this study did not compare lesional HE skin with non‐lesional and healthy control skin.

The most differentially upregulated gene in lesional VHE skin, S100A7A, has previously been identified in inflamed hyperplastic psoriatic skin and is being expressed throughout the epidermis of chronic AD. , Other S100‐family members (S100A7, S100A8 and S100A9), which we also found to be highly expressed in lesional VHE skin, are known to be involved in the regulation of a broad range of cellular processes, and show antimicrobial properties as well. It has been suggested that these S100 calcium binding proteins promote cytokine production, which is critical in the inflammatory response, and lead to decreased filaggrin and loricrin expression levels in the epidermis. , Enhanced antimicrobial and innate defense responses attenuates skin barrier function. Our RNA‐seq data also shows highly increased expression levels of host defense genes such as PI3, LCN2, LTF, and LYZ, which may contribute to an impaired skin barrier function in VHE.

The most downregulated gene in the VHE transcriptome is c5orf46, also known as skin and saliva secreted protein 1 (SSSP1), which is located on chromosome 5q32 in between the family of SPINK genes. Its function is still unknown, but it is hypothesized that c5orf46 interacts with TMBIM6, which modulates endoplasmic reticulum calcium homeostasis. Calcium homeostasis might play an important role in the pathogenesis for multiple skin disorders, including AD, since it is considered to be a central regulator for keratinocyte differentiation. ,

As we have successfully validated the reliability of our RNA‐seq data by RT‐PCR and IHC, we believe it's valid to discuss and evaluate genes from this dataset that have not been analysed by RT‐qPCR. One such example is FABP7, a lipid‐associated gene which is significantly downregulated in VHE. Decreased expression of FABP7 mRNA has recently also been reported for psoriasis and for a form of autosomal recessive congenital ichthyosis. , Fatty acid binding proteins (FABP) are involved in long‐chain fatty acid metabolism, by facilitating their delivery to intracellular sites. The upregulation of family member FABP5 in AD and psoriasis is thought to be a response to the abnormal differentiation and proliferation of keratinocytes. ,  Therefore, the decreased expression levels of FABP7 might suggest that an impaired lipid metabolism is a contributing factor in HE pathogenesis as well.

In addition, we found DPP4 to be upregulated in lesional VHE. DPP4 is an intrinsic membrane glycoprotein, which was found to be inducted by type 2 cytokines like IL‐4 and IL‐13 in type 2 inflammation disorders, including AD, asthma and chronic rhinosinusitis, and has been opted to be a biomarker in type 2 inflammation.

One gene that was strongly upregulated in the VHE transcriptome, SPINK6, may also play an important role in skin barrier function. The protein encoding this gene, serine peptidase inhibitor Kazal‐type 6, has been reported to control the protease activity of kallikreins (KLK5, KLK7 and KLK14) that regulate skin desquamation by cleaving corneodesmosomal components. An imbalance between enzyme activity of serine proteases and their inhibitors might be an contributing factor in the pathogenesis of HE. However, SPINK6 expression is reported to be decreased in AD lesional skin.

Theories on HE pathogenesis are mostly based on AD literature, which mainly focuses on skin barrier defects. The link between HE and AD is still largely unknown, but AD is certainly a risk factor for developing HE. Our data on the VHE skin transcriptome are largely comparable to published gene expression data of AD lesional skin (Table S5). Both are predominantly characterized by the upregulation of keratinocyte host defense and inflammation genes, immune signalling genes, and the downregulation of several skin barrier related genes. The similarities in gene expression may indicate that new treatments for AD might also be effective for HE, as current treatment options for severe HE are limited. Most of the novel therapies for AD are focusing on pathway‐targeted therapies, for example dupilumab, a monoclonal antibody that inhibits signalling of the hallmark AD cytokines IL‐4 and IL‐13. , In the RNA‐seq data set, IL4R was found to be highly upregulated in lesional HE skin vs. healthy control skin. Noteworthy, dupilumab has been found to be effective in HE patients with AD, and in patients with isolated hand eczema. Furthermore, JAK3 was also found to be upregulated in lesional HE skin vs. healthy control skin. However, this gene did not reach the threshold of >500 averaged normalized counts between the groups and, thus, is only very lowly expressed. This finding still may be promising for pan‐Janus Kinase (JAK) inhibitors, e.g. gusacitinib and the topic JAK‐inhibitor delgocitinib. Future research should focus on similar pathways in HE.

In AD lesional skin, KRT6, KRT16 and KRT17 are highly expressed, while KRT1 and KRT10 are downregulated, , most likely due to suppression by IL‐4 and IL‐13. In HE, we also found high expression levels of KRT6, KRT16 and KRT17, but we did not observe downregulation of KRT10 gene expression, which we confirmed by protein expression analysis. On the other hand, we observed a very significant downregulation of KRT2 gene and protein expression, which has not been reported in AD literature. Absence of KRT2 might lead to disturbance of the epidermal integrity of plantar skin, which suggests that this protein has an important role in HE pathogenesis.

It is well known that non‐lesional AD skin significantly differs from healthy skin, including delayed and intermittent expression of barrier genes and their proteins (e.g. LOR, FLG) and increased expression of multiple immune‐related genes (e.g. IL‐13, CCL17, CCL22, S100A7). In the NL skin transcriptome of VHE, no significant differences with healthy skin were found, which is striking as it has been hypothesized that a primary dysfunctional epidermal barrier is one of the key pathways for developing HE. , , , ,

One of the limitations of our study is the small sample size of ten HE patients vs. ten healthy controls. Secondly, it is important to note that we used stringent criteria for the RNA‐seq data as presented in Figure 2 (fold‐change ≤−3 or ≥3, >500 average normalized counts). We could not to compare our study to other HE studies in literature, since these reports included other subtypes or non‐specified subtypes of HE, and different techniques such as immunofluorescence staining.

In conclusion, we have identified a significantly different gene expression profile in palmar lesional skin of VHE patients, compared to non‐lesional skin and skin of healthy controls. The large overlap between the transcriptome of VHE and AD lesional skin suggests that treatments for AD might be also effective in HE, at least in VHE. Transcriptome analysis of HE subtypes other than VHE could provide insight into specific differences between HE subtypes, which may lead to a new classification of HE based on endo(pheno)types, as suggested for AD. ,

CONFLICT OF INTEREST

Dr. Schuttelaar received a research grant from Sanofi Genzyme; received consultancy fees from Sanofi Genzyme and Regeneron Pharmaceuticals; and is advisory board member for Sanofi Genzyme, Regeneron Pharmaceuticals, Pfizer, Abbvie, LEO Pharma, Lilly. The other authors have no conflicts of interest to disclose.

AUTHOR CONTRIBUTION

All authors have read and approved the final manuscript. Conceptualization: M.L.A.S., P.L.J.M.Z., A.N.V., H.N. (equal). Methodology: M.L.A.S. (lead), P.L.J.M.Z. (lead), A.N.V. (lead), H.N. (lead), T.H.A.E. (supporting). Software: T.H.A.E. (lead). Validation: all authors (equal). Analysis: H.N. (lead), P.L.J.M.Z. (supporting), T.H.A.E. (supporting), Investigation: all authors (equal). Resources: all authors (equal). Data curation: all authors (equal). Writing – original draft: A.N.V. (lead), H.N. (lead), M.L.A.S. (supporting), P.L.J.M.Z. (supporting). Writing – review & editing: M.L.A.S. (lead), P.L.J.M.Z. (lead), all other authors (supporting). Visualization: all authors (equal). Supervision: M.L.A.S, P.L.J.M.Z. (equal). Project administration: all authors (equal). Funding acquisition: A.N.V., M.L.A.S. (equal).

App S1. Supplemental information

Click here for additional data file.

File S1. RNA‐seq results

Click here for additional data file.

File S2. RNA‐seq alignment statistics

Click here for additional data file.

File S3. Supplementary table GSEA based on GO (BP)

Click here for additional data file.