Keratinocyte-Immune Cell Crosstalk in a STAT1-Mediated Pathway: Novel Insights Into Rosacea Pathogenesis

Rosacea is a common chronic inflammatory condition that mainly affects the central face. However, the molecular background of the normal central face and the transcriptional profiling and immune cell composition of rosacea lesions remain largely unknown. Here, we performed whole-skin and epidermal RNA-seq of central facial skin from healthy individuals, lesions and matched normal skin from rosacea patients. From whole-skin RNA-seq, the site-specific gene signatures for central facial skin were mainly enriched in epithelial cell differentiation, with upregulation of the activator protein-1 (AP1) transcription factor (TF). We identified the common upregulated inflammatory signatures and diminished keratinization signature for rosacea lesions. Gene ontology, pathway, TF enrichment and immunohistochemistry results suggested that STAT1 was the potential core of the critical TF networks connecting the epithelial–immune crosstalk in rosacea lesions. Epidermal RNA-seq and immunohistochemistry analysis further validated the epithelial-derived STAT1 signature in rosacea lesions. The epidermal STAT1/IRF1 signature was observed across ETR, PPR, and PhR subtypes. Immune cell composition revealed that macrophages were common in all 3 subtypes. Finally, we described subtype-specific gene signatures and immune cell composition correlated with phenotypes. These findings reveal the specific epithelial differentiation in normal central facial skin, and epithelial–immune crosstalk in lesions providing insight into an initial keratinocyte pattern in the pathogenesis of rosacea.


INTRODUCTION
Rosacea is a common disorder of the facial skin that affects about 5.5% of individuals in the general population (1). It is divided into four subtypes: erythematotelangiectatic (ETR), papulopustular (PPR), phymatous (PHR) and ocular. Convex parts of the central face are the region of the skin most frequently affected by rosacea. In fact, human skin shows remarkable diversity in its structure and function across anatomic sites (2). The embryologic development of the facial dermis derives from the neural crest, so facial skin has a distinct developmental origin different from other anatomic sites (3). Recently, a new concept of the immunological anatomy of skin highlighted the epithelialimmune interaction of skin in different anatomic sites (4). However, the molecular characteristics of the central facial skin and their potential relation with the pathogenesis of rosacea remain unknown.
Although the pathogenesis of rosacea is not fully understood, it is defined as a chronic, inflammatory skin condition. Abnormalities in the innate immune response system and neurovascular dysregulation are considered primarily implicated in the pathophysiology of rosacea (5). Immune cells involved in the innate immune response, such as macrophages, mast cells, and neutrophils, were described in rosacea lesions (6). However, the immune cell composition as well as the relative change in abundance of specific cell types in rosacea lesions are still largely undefined.
Dynamic epithelial-immune crosstalk fine-tunes epithelial homeostasis (7), which suggests that epithelial cells (keratinocytes) might be the primary trigger of the innate immune response at the onset and development of rosacea. However, evidence of keratinocyte-immune cell crosstalk in the rosacea lesion and what initially triggers this crosstalk remain poor unknown. Dynamic rosacea appears to feature a consistent inflammatory continuum, and a ''developmental march'' among different subtypes was assumed (8). Regardless, little is known about the distinct inflammatory properties of different clinical subtypes and their relationship with clinical manifestations.
To address these gaps in our understanding, the transcriptional profiling and the immune cell composition in rosacea, we used RNA-sequencing (RNA-seq) of central facial skin from healthy individuals and a large number of central facial skin lesions and matched normal skin surrounding the auricle from patients with different rosacea subtypes.

Sample Collection
Skin biopsies (within 0.5 cm 2 ) were obtained from the central facial skin of female healthy individuals undergoing plastic surgery (n=19 for whole skin RNA extraction; n=5 for epidermis RNA extraction; epidermis was isolated from whole skin with dispase) and lesions of the central face and corresponding normal skin surrounding the auricle from female patients with ETR (n=15 for whole skin; n=6 for epidermis), PPR (n=22 for whole skin; n=8 for epidermis) and PHR (n=9 for whole skin; n=4 for epidermis) (aged 20-60 years) from the Department of Dermatology in Xiangya Hospital, Central South University. The study was approved by the ethics committee of Xiangya Hospital, Central South University, and written informed consent was obtained from all participants. The experiments conformed to the principles set out in the WMA Declaration of Helsinki and the Department of Health and Human Services Belmont Report.

Library Construction and Sequencing
Whole skin or epidermis RNA was isolated by using TRIzol Reagent (Thermo Fisher Scientific, USA). A total of 3 µg RNA per sample was used as input material for RNA sample preparations. Sequencing libraries were generated by using the NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer's recommendations. The Illumina PE150 libraries were sequenced on a Hiseq 4000 platform according to the manufacturer's recommendations (Illumina) at the Novogene Bioinformatics Institute, Beijing. Finally, 150-bp paired-end reads were generated.

Quantification of Gene Expression
Raw data (raw reads) in the FASTQ format were first processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapters and ploy-N and low-quality reads from raw data. Raw reads were aligned to the human reference genome (hg19). Sequencing data were checked for sequencing quality by using FASTQC. FeatureCounts v1.5.0-p3 was used to count the read numbers mapped to each gene. Gene expression was estimated by using fragments per kilobase transcript sequence per million base pairs sequenced (FPKM) values.

Assessment of Immune Cell Types
Relative levels of distinct immune cell types in rosacea lesions were quantified by using the xCell algorithm (13) and CIBERSORT algorithm (https://cibersort.stanford.edu) (14).

Statistical Analysis
Statistical analysis was conducted with SPSS 20.0 and GraphPad 8.0. Data are presented as the mean ± SEM. We evaluated the data for normal distribution and similar variance between groups. Statistical significance (*P < 0.05, **P < 0.01) was assessed by 2-tailed unpaired Student's t-test for comparisons between 2 groups and 1-way analysis of variance (ANOVA) with relevant post hoc tests for multiple comparisons. We performed the 2-tailed Mann-Whitney U test for statistical analysis when the data were not normally distributed or exhibited unequal variances between the two groups.

Normal Central Facial Skin Presents Anatomy-Specific Gene Signatures Involved in Epithelial Cell Differentiation
In an initial analysis, we used transcription profiling of central facial skin from healthy individuals (healthy skin [HS], n=17) compared to normal skin surrounding the auricle (normal skin [NS], n=45) ( Figure 1A). We identified 984 differentially expressed genes (DEGs) (fold change in expression > 2 and false discovery rate [FDR] < 0.05): 799 upregulated and 185 downregulated in central facial skin ( Figure 1B). Hierarchical clustering analysis distinguished HS from NS, thus indicating a different gene expression pattern between the two anatomic sites ( Figure 1C).
To examine the transcriptional mechanisms underlying the DEGs, we analysed transcription factors (TFs). A total of 74 TFs were identified, including 21 downregulated and 53 upregulated, in HS. The top clusters and networks involved the biological processes epithelial cell differentiation (Log10(q)= -6.81), cell fate commitment (Log10(q)= -6.74) and regionalization (Log10 (q)= -6.74) (Figure 2A). The MCODE algorithm subclustered the PPI network into 2 subclusters: the activator protein-1 (AP1) pathway (Log10(p)= -11.5) and cellular response to hormone stimulus (Log10(p)= -9.2) ( Figures 2B, C). The co-expressed TFs had the most significantly enriched functional annotations ( Figure 2D). We developed a scheme to classify TFs in terms of their biological function: downregulated TFs were mostly clustered in neural development and upregulated TFs in nerve growth and the AP1 pathway ( Figure 2E). Furthermore, we identified AP1 TFs that were covered in the ChEA3 database, including JUN, JUND, JUNB, FOS, FOSB, ATF3 and MAFF ( Figure 2G). By immunohistochemistry, we further confirmed that AP1 pathway was significantly upregulated in the central facial skin (Supplementary Figure S1). Thus, gene signatures involved in metabolism of lipids and epithelial cell differentiation (cornification) were characteristics of central facial skin, and AP1 TFs may be required for maintaining homeostasis in HS.
To accurately and comprehensively estimate the cell types of HS, we used absolute deconvolution with the xCell algorithm. As compared with skin surrounding the auricle, in HS, the number of basophils and CD4 + and CD8 + T cells was significantly increased ( Figure 2F), which indicates an inflammationsensitive background.
We reveal a common gene signature derived from upregulated genes, enriched in biological processes including inflammation, and downregulated genes, involved in PPAR signalling and keratinization.
Furthermore, the ChEA3 database covered 9 of 11 TFs identified, including STAT1, lymphocyte differentiation genes (IKZF1, TBX21 and IRF1), immune response genes (TFEC and BATF2) and humoral immune response genes (IKZF3, SPIB and IRF8) (Figures 4D, E). Given that the core transcription regulatory circuitry (CRC) consists of a small group of self-regulated TFs and their interconnected regulatory loops, dbCoRC analysis showed a STAT1-containing core regulatory circuitry, including the interactions STAT1-interferon regulatory factor 1 (IRF1) and STAT1-IRF8 ( Figure 4F). We manually analysed the correlation between TFs and constructed the regulatory network ( Figure 4G). Consistently, by immunohistochemistry we verified the increased nuclear localization of phosphor-STAT1, IRF1 and IRF8 in both epidermis and dermal infiltration of rosacea lesions ( Figure 4H).

Epidermal RNA-Seq Data Further Revealed an Epithelial-Derived STAT1 Signature in Rosacea Lesions
To confirm whether the STAT1 signature derived from epidermis contributed to rosacea lesions, we used epidermal RNA-seq and identified a gene set overlapping between LS (n=18) versus HS (n=5) and LS (n=18) versus NS (n=18), with 264 genes showing elevated expression in both comparisons. For the upregulated genes, the enriched GO terms included defense response to other organism (Log10(q)= -23.12) and response to bacterium (Log10(q)= -17.02). Among the 264 upregulated genes were 12 overlapped TFs in rosacea lesions, including STAT1, FOXE1 and ERG ( Figure 5A). These data confirmed that the epithelial-derived STAT1 signature contributed to rosacea lesions.

Epidermal IFNg/STAT1/IRF-1 Shared Across ETR, PPR and PHR Subtypes
To identify the common epidermal gene signatures across rosacea subtypes, we compared ETR (n= 6), PPR (n= 8), and PHR (n=4) subtypes with normal central facial skin from healthy individuals (HS, n=5). Enrichment analysis of upregulated genes revealed an IFN-g-mediated signaling pathway across ETR, PPR and PHR subtypes ( Figures 5B-D). We further observed the common 65 upregulated genes shared across the 3 subtypes: 4 TFs were identified, namely STAT1, IRF1, MECOM and SP110 ( Figure 5E). By immunohistochemistry we also demonstrated that the nuclear localization of phosphor-STAT1 and IRF1 was increased in epidermal cells of rosacea lesions ( Figure 4H). These results demonstrated the STAT1/IRF1 signature across the ETR, PPR and PHR subtypes.

Subtype-Specific Gene Signature
To elucidate the gene signature differences among ETR, PPR and PHR subtypes, we compared one subtype with the other two subtypes. The heatmap ( Figure 6A) shows that each subtype has a unique expression pattern. To determine whether these subtypes are specific to biological processes, we examined upregulated gene enrichment. As compared with PPR and PHR subtypes, ETR showed strong enrichment for biological processes related to muscle contraction (Log10(q)= -4.63) and cornification (Log10 (q)= -1.53). As compared with ETR and PHR subtypes ( Figure 6B), PPR showed a strong enrichment of both innate (Log10(q)= -11.72) and adaptive immune responses (Log10(q)= -7.55) ( Figure 6C). We also noted a strong enrichment of formation of cornified envelope (Log10(q)= -8.99) in the PPR subtype. Of note, PHR was mainly related to metabolism of lipids (Log10(q)= -18.32), monocarboxylic acid metabolic process (Log10(q)= -15.49) and organ acid catabolic process (Log10(q)= -12.10) ( Figure 6D).

Immune Cell Composition for ETR, PPR and PHR Subtypes
Using the CIBERSORT algorithm, we first investigated the common and different immune infiltration features among ETR, PPR and PHR subtypes ( Figure 6E). Resting dendritic cells and macrophages were common to all 3 subtypes. In the ETR subtype, the most abundant immune cells were resting dendritic cells (18.3%), regulatory T cells (12.6%) and resting mast cells (11.9%), whereas in the PPR subtype, the most abundant immune cells were M0 macrophages (13.9%), resting dendritic cells (12.7%) and regulatory T cells (11.9%). The top 3 abundant immune cell types in the PHR subtype were resting dendritic cells (18.8%), resting CD4 + memory T cells and regulatory T cells (11.1%). As compared with PPR and PHR  subtypes, ETR showed a higher proportion of resting mast cells, whereas as compared with ETR and PHR subtypes, PPR showed abundant activated memory CD4 + and CD8 + T cells by xCell and CIBERSORT algorithms ( Figure 6F).

DISCUSSION
Previous whole-skin transcriptome analysis revealed that Th1/ Th17 polarized inflammation and macrophage infiltration were hallmarks in ETR, PPR and PHR subtypes (6), and separated signals of immune and keratinization gene expression in PPR (17). In addition, prominent permeability barrier alterations had also been descried in PPR via whole-skin RNA sequencing analysis (18). In the present study, through both whole-skin and epidermal transcriptome analysis, we characterized the atlas of transcriptomic profiles indicating the site-specific gene signatures and epidermal homeostasis for central facial skin. We further identified common inflammatory signatures, the TF network in whole-skin and epidermis, and epidermal STAT1/ IRF1 signature and epithelial-immune crosstalk for rosacea lesions ( Figure 5F). We also compared the molecular characteristics and immune cell composition for different rosacea subtypes based on their transcriptome profiles.

Normal Central Facial Skin Exhibited Anatomy-Specific Molecular Characteristics Involved in Metabolism of Lipids, Steroid Biosynthetic Process and Epithelial Cell Differentiation
We compared transcription profiling of central facial skin from healthy individuals and normal skin surrounding the auricle. Upregulated gene signatures were mainly associated with metabolism of lipids, steroid biosynthetic process and epithelial cell differentiation. We identified enrichment of the epithelial-cell differentiation pathway and AP1 TFs (FOS, FOSB, JUN and JUN) in central facial skin. AP1 TFs are key regulators of epidermal keratinocyte survival and differentiation (19). The epidermis contains proliferating keratinocytes that differentiate to form the corneal layer. Corneocytes together with intercellular lipids isolate the body from the external environment and protect it against infection. Thus, persistent activation of AP1 pathways may be necessary to maintain keratinocyte differentiation and the structural barrier of central facial skin, and its dysfunction might lead to rosacea. Our results were consistent with the hypothesis that AP1 TFs (FOS, FOSB, JUN and JUND) were downregulated in rosacea lesions. This situation may reflect intrinsic defects accelerating keratinocyte differentiation capacity, which may contribute to compromised barrier function in rosacea lesions.
To summarize, the DEGs associated with keratinocyte differentiation are required for maintenance of the local homeostasis of central facial skin and might explain the regional difference in rosacea susceptibility.

STAT1 Mediated Inflammatory Signatures for Rosacea Lesions
Although the clinical manifestations of rosacea are heterogeneous, they are all related to the presence of characteristic facial or ocular inflammation involving both the vascular and tissue stroma (20). We attempted to present a common link between the triggers and the immune response of rosacea lesions and used a comprehensive comparison among central facial skin, skin surrounding the auricle and rosacea lesions to identify shared gene signatures for rosacea lesions ( Figure 5F).
The common upregulated gene signatures of rosacea lesions were mainly associated with innate immune responses including the defense response to other organisms and response to IFN-g, which promotes macrophage activation, mediates antiviral and antibacterial immunity, enhances antigen presentation, and orchestrates activation of the innate immune system. Consistently, we identified 8 of 12 TFs (STAT1, IRF1, IRF8, IKZF1, IKZF3, TBX21, EOMES and BATF2) involved in IFN-g signaling. We also confirmed the interaction of the 8 TFs by co-expression and functional enrichment. Finally, core TF loop analysis revealed that STAT1 is the core TF of the IFN-g/STAT1 pathway.
Familial rosacea associated with a gain-of-function mutation (C324R) in STAT1 suggests that STAT1 could lead to chronic facial inflammation as the cause of rosacea (21). Innate immunesystem perturbations are known to contribute to the pathogenesis of rosacea (20). The IFNg/STAT1/IRF-1 signaling pathway was previously found involved in the processes of keratinocyte differentiation and the inflammatory interplay between innate immune responses and keratinocytes (22). The macrophage IRF8/IRF1 regulome was found required for protection against infection and is associated with chronic inflammation (23). IKZF1 regulates dendritic cell development and cytokine signaling pathways and CD4 + Th cell differentiation (24). TBX21 and EOMES are the only T-box proteins expressed in the immune system and control key checkpoints of natural killer cell maturation (25). We also identified the up-regulated transcription factors enriched in the IL12 signaling pathway, and up-regulated expression of IL12RB2 in rosacea lesions. The expression of IL12RB2 is up-regulated by IFN-g in Th1 cells and plays a role in Th1 cell differentiation via leading to the formation of high-affinity IL12 binding sites and reconstitution of IL12 dependent signaling (26). Overall, distinct transcriptional profiles are linked to inflammation at lineagespecific transcriptional regulators in rosacea lesions. Taken together, we suggest that IFNg/STAT1 is the core of the critical gene regulatory network connecting the epithelial-immune crosstalk in rosacea development, and targeting IFNg/STAT1 might be a potential therapy for the treatment of rosacea.
Epidermal-Derived IFNg/STAT1/IRF1 Signature Contributes to Rosacea Lesions As immune sentinels, keratinocytes produce innate immune mediators and act as non-professional antigen-presenting cells (27). Exposure to acute solar ultraviolet (SUV) leads to upregulation of IFNg and downstream STAT1/IRF1 signaling in the epidermis (22). These findings suggest that the IFNg/ STAT1/IRF1 signature is derived from epidermis in rosacea. According to our epidermal RNA-seq data, the upregulated genes were mainly related to the IFNg signaling pathway for rosacea lesions. Furthermore, STAT1 and IRF1 were identified across all rosacea subtypes. Our findings, together with previous reports, support that epidermis is pivotal for immune cell activation, which in turn leads to amplification of local innate immune responses. Similar "transcription circuits" including STAT1/IRF1 were found in psoriatic epidermis (28). Our findings provide a new model that diverse inflammatory processes in rosacea may be driven largely by just a small number of hubs within the epidermal transcription circuitry, such as STAT1/IRF1.
Given the importance of epithelial-immune crosstalk, there is value in measuring the expression from immune cells in rosacea lesions. Deconvolution is an accurate method for detecting a signal that is specific for a cell type and is expressed in a consistent pattern among the cell types of a heterogeneous sample (14). Resting dendritic cells and macrophages were common in ETR, PPR and PHR subtypes, which highlight the participation of innate immunity in rosacea pathogenesis. Furthermore, in terms of the effect of STAT1/IRF1 on M1 polarization, our previous study found that ADAMDEC1 plays a pro-inflammatory role in rosacea by modulating the M1 polarization of macrophages (29). Macrophages and keratinocytes influence the expression of inflammatory modulators that manifest clinically in the formation of pustules, papules, and sensations of heat and pain (30). The dialogue between keratinocytes and immune cells is fine-tuned to facilitate coordinated responses to preserve homeostasis and mount a host defense (7). Our data suggest that over-active keratinocyte-macrophage crosstalk might be an argument for activation of the innate immune system for rosacea development.

Subtype-Specific Gene Signatures Correlated With Clinical Phenotypes
Various subtypes can be manifested in one patient, thus leading to overlapping or combined clinical manifestations. To gain exact molecular insight into rosacea subtypes, we collected typical erythema, papules and phymatous lesions from ETR, PPR and PHR subtypes and compared the rosacea subtypes with healthy controls. We found that cornification and formation of the cornified envelope in ETR and PPR indicated epidermal dysfunction (skin barrier) as an early phase of rosacea. Previous study of the PPR subtype found it characterized by a profoundly diminished skin barrier (18); our results further support the prominent permeability barrier alterations in ETR at the molecular level, which highlight the diminished skin barrier as a potential initiator of rosacea.
We further explored the subtype-specific characteristics by comparing one subtype with the other two. The ETR subtype featured striated muscle contraction, which is consistent with the clinical symptoms of frequent episodes of transient facial erythema and/or persistent erythema. The immunoregulatory interaction was predominantly identified in the PPR subtype, which indicates an inflammatory infiltrate leading to papules and pustules. Compared with the PPR subtype, the down-regulated genes enriched in the GO terms of keratinization, neutrophil chemotaxis, monocyte chemotaxis and extracellular matrix disassembly in PHR subtype. The down-regulated genes involved in monocyte chemotaxis including CCL4 and CCL5, indicated the macrophage activation from M1 to M2 type (31), which plays critical roles in tissue remodeling and repair. Our data is consistent with the previous study (32), which indicated non-phagocytic macrophages favors a role in limiting fibrosis. PHR also exhibited upregulation of many genes involved in lipid metabolism, thus indicating sebaceous gland hyperplasia with a strong connection between chronic inflammatory processes and skin fibrosis development. Our work provides a detailed characterization of transcriptomic diversity of rosacea subtypes and additional evidence that this can be linked to phenotypes, from epidermis (ETR and PPR) to sebaceous gland (PHR), vascular reaction (ETR) and cell infiltration (PPR) of acute inflammation to hyperplasia (PHR) of chronic inflammation.
Among the immune cell types, the transcriptome composition of mast cells was particularly high in the ETR subtype, which suggests that mast cells participate in the initial phase of rosacea evolution. Mast cell proteases recruit immune cells to amplify the inflammatory response (cathelicidin LL-37 activation) and cause vasodilation (33). For inflammatory outbreaks in the PPR, CD4 + and CD8 + T cells and B cells were predominant in PPR, which supports the cell-mediated immune response. As compared with ETR and PPR, PHR shows absence of distinct immune cell composition. Here our novel immune-cell composition analysis provided new insights into the abundance of immune cell types of rosacea subtypes from patients with early erythema to fibrotic changes.
In conclusion, we extensively and comprehensively investigated the characteristics of central facial skin and rosacea lesions in a stepwise manner. We potentially revealed the communication between the skin barrier and innate immune system, keratinocyte-immune crosstalk, and provided insights into the initial keratinocyte pattern for the pathogenesis of rosacea. Finally, our findings suggested STAT1 might be a potential therapeutic target for rosacea.

DATA AVAILABILITY STATEMENT
All data needed to evaluate the conclusions in the study are provided in the manuscript and/or the Supplementary Materials, further inquiries can be directed to the corresponding author/s. Whole-skin and epidermal sequencing data from rosacea patients have been deposited in the genome sequence archive under accession number HRA000379 and HRA000809 (http:// bigd.big.ac.cn/gsa-human/browse).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the ethics committee of Xiangya Hospital, Central South University. The patients/participants provided their written informed consent to participate in this study.