Comparative Transcriptome Analysis of Acne vulgaris, Rosacea, and Hidradenitis Suppurativa Supports High Dose Dietary Zinc as a Therapeutic Agent

Identification of differentially expressed genes and enriched pathways through transcriptome analysis

We first conducted differential expression analysis on the eleven datasets (Fig. 1b and c; Supplementary Fig. 1b and c), and filtered DEGs using a threshold of adjusted p value [<[0.05 and |log2 Fold Change|[>[1. The examination of the Acne vulgaris dataset (GSE6475) revealed 96 upregulated genes and 32 downregulated genes in lesional samples compared to healthy controls, along with 103 upregulated genes and 9 downregulated genes in lesional samples compared to non-lesional controls. From PPR GSE65914, we identified 767 upregulated genes and 492 downregulated genes in lesional samples when compared with healthy controls. From PPR Dataset 2, 345 upregulated genes and 84 downregulated genes were identified in lesional samples compared to non-lesional controls. In HS GSE148027, our analysis revealed 1226 upregulated genes and 883 downregulated genes in lesional samples compared to healthy controls, as well as 1303 upregulated genes and 823 downregulated genes in lesional samples compared to non-lesional controls. In Acne vulgaris, the most prominently upregulated genes (Lesional versus Healthy) encompassed matrix metalloproteases (MMP1, MMP3, MMP12), members of the serpin family (SERPINA1, SERPINB3, SERPINB4), and genes from the epidermal differentiation complex (SPPR gene family, S100 gene family). Notably, PPR and HS also exhibited substantial upregulation in these genes. Among the most consistently downregulated genes across all three datasets were the motility-promoting protein TSPAN8 and the anti-inflammatory IL-37 in lesional samples compared to healthy controls. This pattern was mirrored in both the PPR and HS datasets, which also demonstrated marked downregulation of the Wnt signaling inhibitor WIF1, the antimicrobial DCD, and the KRT77 genes.

We then conducted gene set enrichment analysis (GSEA) on the three datasets to identify enriched Gene Ontology and KEGG pathways (Fig. 1d for lesional versus healthy controls, Supplementary Fig. 1d for lesional versus non-lesional controls). All three diseases showed significant upregulation of immune response pathways, such as innate and adaptive immune cell recruitment and activation, associated with regulatory/inhibitory immune responses such as regulatory T cell differentiation. Furthermore, all three diseases demonstrated heightened serine-type endopeptidase activities, leading to increased extracellular matrix disassembly. All three diseases exhibited downregulation in genes involved in the Wnt signaling, especially for genes that possess Wnt inhibitor functions (including WIF1 and DKK2), suggesting a dysregulation of the pathway in the diseases. Both Acne and HS showed downregulation of lipid-related pathways, such as lipid synthesis and modification. PPR showed upregulation in those pathways when comparing to healthy controls but downregulation when comparing to non-lesional controls.

Identification of disease-specific molecular features

Despite the shared clinical manifestations, Acne vulgaris, PPR, and HS represent distinct diseases with potentially differing etiologies. To gain deeper insights into their molecular distinctions, we conducted a comparative analysis of the DEGs across these three conditions, identifying DEGs that are either exclusive to each disease or commonly shared (Fig. 2a; Supplementary Fig.2a). Our approach began by comparing the DEGs lists sourced from the three primary datasets (Acne vulgaris GSE6475, PPR GSE65914, and HS GSE148207), and then validated these genes by cross-referencing them with the additional eight datasets listed in Supplementary Fig. 1a. (e.g., we verified that the DEGs unique to Acne vulgaris did not appear in the DEG lists from all PPR and HS datasets). In the case of PPR, there were 216 uniquely upregulated and 209 uniquely downregulated DEGs. Notably, HS exhibited the most extensive set of exclusive DEGs, encompassing 651 upregulated and 658 downregulated genes. This pattern of distinct DEGs remained consistent across all diseases when comparing lesional samples to non- lesional controls (Supplementary Fig.2a).

• Download figure
• Open in new tab

Figure 2. Identification of molecular features that are unique to each disease.

a) Strategy for identification of unique DEGs. Venn plots were generated by overlaying DEGs from the three diseases, considering upregulated and downregulated DEGs separately. DEGs that were unique to each disease were validated against additional datasets to ensure the uniqueness. The validated unique DEGs were quantified and visualized using a bar graph. b) A protein- protein interaction (PPI) network of DEGs from unique PPR upregulated genes that are enriched in the lipid metabolism pathway. c) A PPI network of DEGs from unique PPR downregulated genes that are enriched in the Homeobox gene family and the extracellular matrix pathway. d) A PPI network of unique HS upregulated genes that are enriched in extracellular matrix organization and collagen biosynthesis and modifying enzymes. e) A PPI network of unique HS upregulated genes that are enriched in B cell and T cell activation and regulation pathways. f) Upregulation of gene expression in B cell/plasma cell-related genes in HS comparing to PPR and Acne vulgaris.

We then performed STRING protein-protein interaction network analysis on the DEGs that were specific to PPR and HS (Fig. 2b-f and Supplementary Fig. 2b-e). When comparing lesional to the healthy controls, PPR displayed uniquely upregulated genes that were enriched in lipid metabolism pathways (Fig. 2b). This observation aligns with the GSEA comparison results depicted in Figure 1D, which also highlighted the differential upregulation of lipid metabolism-related pathways in PPR. Intriguingly, we also detected a distinctive downregulation of development related Homeobox family genes in PPR, including HOXA1, HOXA10, HOXB2, HOXB7, HOXC4, HOXC9, HOXC10, and their interacting partners PBX3 and ISL1 (Fig. 2c). Within the same downregulated gene set, we also noticed an enrichment in extracellular matrix proteins, such as COL1A2, COL3A1 and COL6A5, which contrasted with an upregulation in extracellular matrix organization pathways in HS (Fig. 2c and d; Supplementary Fig.2b). In the context of comparing lesional to non-lesional samples, we identified that PPR displayed unique downregulation of the late cornification envelope (LCE) family genes implicated in keratinization (Supplementary Fig. 2c).

HS exhibited a more pronounced activation of adaptive immune responses when compared to PPR and Acne vulgaris, particularly in B cell/plasma cell responses (Fig. 2e and Supplementary Fig. 2d). We then scrutinized the original DEG list for HS and discovered that among the top 50 most-upregulated genes (ranked by Log2 fold change), 24 were associated with B cell/plasma cell functions (Fig. 2f). These same genes demonstrated relatively lesser upregulation in PPR and displayed minimal or no upregulation in Acne vulgaris. Of interest, we conducted a comparison between non-lesional samples and healthy controls within the HS dataset. Strikingly, we identified an upregulation of B cell-related genes even in non-lesional samples, implying a systemic elevation of B cell response among HS patients (Supplementary Fig. 2e).

Identification of disease-shared molecular features

Acne vulgaris, PPR and HS are all inflammatory skin disorders, and understanding of the commonality of the molecular features of the three diseases may foster the development of broadly applicable therapeutics. We identified the DEGs that were shared by all three diseases as depicted in Figure 2a and validated them by selecting the DEGs that were also present in the DEG lists of the additional eight datasets. Ultimately, we confirmed the validation of 22 upregulated DEGs, including antimicrobial peptides (S100A7, S100A8, S100A9 and S100A12), chemokines (CXCL1, CXCL8 and CCL19), Matrix Metallopeptidase (MMP1 and MMP12), and protease inhibitor (PI3, SERPINA1 and SERPINB4) (Fig. 3a). We performed STRING network analysis on these 22 DEGs and the results revealed the genes were enriched in stress response pathways, chemotaxis, IL-17 signaling pathways, metal sequestration antimicrobial proteins, serine-type endopeptidase activity and serine protease inhibitor (Fig. 3b). No shared downregulated DEGs passed validation.

• Download figure
• Open in new tab

Figure 3. Identification of molecular features that are shared by the three diseases.

a) Visualization of the 22 validated shared DEGs. The dot size represents the fold change of gene expression, and the dot color scale represents the adjusted P value. b) A protein-protein interaction network of the 22 genes. The enriched pathways were highlighted in colors.

Identification and validation of high dose dietary zinc as a therapeutic agent

To identify potential therapeutic agents, we focused on the shared DEGs and map the gene-drug interactions using the Drugbank database (https://go.drugbank.com/). We observed that multiple genes (SERPINA1, S100A7, S100A8, S100A9 and KRT16) interacted with zinc and its related salts (zinc sulfate, zinc chloride and zinc acetate) (Fig. 4a). Zinc is an essential mineral that plays a crucial role in various bodily function, such as immune system regulation. Although zinc supplementation has been explored for the treatment of inflammatory skin disorders, the outcomes have yielded mixed results (Dhaliwal et al. 2020). An examination of the effects of zinc supplementation on inflammatory skin diseases divulged that out of 14 studies focusing on zinc supplementation for acne vulgaris patients, 10 exhibited a beneficial impact on disease severity. Similarly, for HS patients, 3 out of 3 studies demonstrated favorable outcomes from zinc supplementation. However, a study involving rosacea patients failed to demonstrate positive results in comparison to the placebo group. It’s worth noting that these studies have certain limitations, such as varied methodologies with different types of zinc salts and dosages (Dhaliwal et al. 2020). Furthermore, the mechanism underlying zinc’s potential efficacy remains elusive, underscoring the need for further research to comprehensively assess the effectiveness of zinc in treating skin inflammation, especially on the mechanisms underlying the beneficial effects at the molecular level.

• Download figure
• Open in new tab

Figure 4. Identification and validation of high dose dietary zinc as a therapeutic agent.

a) Gene-drug interaction analysis identified zinc and its related salts interacted with multiple validated shared DEGs. b) Schematic of the experimental design. c) Effect of dietary zinc on disease scores and bacteria load in an acne mouse model that was infected with C. acne HL043PA1 strain. d) Quantification of inflammatory cytokines from skin homogenates collected from infected mice under diets with different amount of zinc. All data were presented as medians from two independent experiment (n = 5 for each experiment).

We have previous developed an acne mouse model (Kolar et al. 2019). To directly query the role of zinc in acne, we maintained CD1 mice on a zinc-deficient diet (0 ppm, parts per million) and compared them with the age-matched mice subjected to diets either sufficient in zinc (30 ppm) or featuring a higher zinc concentration (300 ppm) over a 40-day period (Fig. 4b). The mice were challenged with C. acnes bacteria intradermally, and disease score and tissue cytokines production were surveyed on Day 3 (Fig. 4c and d). Independent of bacterial burden, mice fed on a Zn deficient diet showed significantly higher disease score and inflammation than the other two groups of mice, which showed comparable acne immunopathology. Among the mice fed with zinc, those receiving a high dose (300 ppm) demonstrated more substantial improvements than their counterparts who received a sufficient dose (30 ppm) of zinc.

We tested two C. acnes strains in this study, and notably, we found significant differences in immunopathology with HL043PA1 strain only, while disease score and inflammation were trending upwards in zinc deficient mice when acne was induced with less acnegenic strain HL043PA2 (Supplementary Fig. 4). Overall, our findings suggest that dietary zinc could play a significant role in mitigating acne immunopathology by controlling excessive inflammatory phenotype observed in Zn deficient mice.

Microarray and RNA-Seq data collection

Microarray and RNA-Seq studies published until July 2023, focusing on skin samples associated with Rosacea, Acne, and hidradenitis suppurativa, were sourced by querying the Gene Expression Omnibus (GEO) (Barrett et al. 2012) using the following key words: (Acne OR Rosacea OR hidradenitis suppurativa OR HS) AND (RNA-seq OR Expression profiling OR high throughput sequencing). The result was manually filtered that are only from expression profiling by Affymetrix array and Illumina high-throughput sequencing. Raw expression data and accompanying metadata were acquired using the GEOquery package v2.68.0 (Davis and Meltzer 2007). A complete list of studies included in the analysis can be found in Supplementary Figure 1a. Each Sample was manually labeled as “Healthy”, “Non-lesional” or “Lesional”, corresponding to healthy control, non-lesional sample from patients, and lesional sample from patients, respectively.

Gene expression quantification for microarray data and RNA-seq data

The microarray data were pre-processed through the Bioconductor affy package v1.78.0 in R to generate count matrix with Robust Multi-Array Average (RMA) expression measure (Davis and Meltzer 2007). The data were then annotated by hgu133plus2.db or corresponding annotation table by query the GPL accession number in GEO database (Carlson et al. 2016). The differential expression analysis was then performed by the limma package v3.56.2 in R. The up- and down-regulated differentially expressed genes (DEGs) were selected from each dataset using a threshold of log2 fold change of 1 and an adjusted p–value (Bonferroni correction) of 0.05 (Ritchie et al. 2015). For RNA-seq data, raw counts were initially annotated using Genome assembly GRCh38.p13. Differential expression analysis was subsequently conducted on the raw counts using the DESeq2 package v1.40.2, employing a negative binomial generalized linear model (Love, Huber, and Anders 2014). The same fold change and adjusted p-value thresholds applied to the microarray data analysis were used here as well. Notably, the outcome of differential expression analysis from rosacea dataset 2 was directly acquired from the Mendeley data website. (Shih et al. 2020). The Venn plots were created using tools from https://bioinformatics.psb.ugent.be/webtools/Venn/.

Gene set enrichment analysis

The analysis was conducted using the R package clusterProfiler v3.12.047. Gene sets and molecular signatures utilized in this analysis were obtained from the Molecular Signatures Database. Enrichment scores were normalized and calculated through 1000 gene set permutations. Subsequently, significant enrichment outcomes were filtered based on a P-value cutoff of 0.05.

Protein-protein interaction (PPI) network analysis

The STRING database (https://string-db.org) was employed to predict functional interactions among the differentially expressed genes. The network type was designated as a full string network. In the case of HS, a high confidence threshold (0.700) and high FDR stringency (1 percent) were applied. For the remaining cases, a medium confidence threshold (0.400) and medium FDR stringency (5 percent) were used.

Ethics statement

All mouse studies were approved under the guidelines of the University of California San Diego (UCSD) Institutional Animal Care and Use Committee (S18200). Outbred 6 weeks old female CD1 mice (The Charles River Laboratory) were housed in an animal facility at UCSD with a standard of care as per federal, state, local, and NIH guidelines.

C. acnes bacterial culture

This study involved the utilization of two acne-associated strains, namely HL043PA1 and HL043PA2. As established in our previous investigation [1], these strains are representative of the most acnegenic strains within the ribotype-type 4/5 panel utilized in our acne mouse model (Kolar et al. 2019). The clinical C. acnes strains sourced from frozen stock were cultivated anaerobically on blood agar plates using the BD BBL™ GasPak™ system for a duration of 96 hours at 37 °C. A single colony of C. acnes was anaerobically grown in 10 ml of Brain Heart Infusion (BHI) broth (Catalog no. #53286, Sigma-Aldrich, USA) for 3-4 days (OD=0.15-0.3), followed by washing once of the bacterial pellet with BHI media at 2300 × g for 5 minutes. The pellet was resuspended in BHI media to an OD of 0.5 for in vivo mouse infection.

Diets

Six-week-old mice were assigned randomly to three distinct experimental groups. These groups were characterized by their dietary intake: a zinc-deficient diet (egg white-based AIN-93G without added zinc; a diet with a 30 ppm concentration of zinc (egg white-based AIN-93G with added zinc); or a diet with 300 ppm concentration of Zn (egg white-based AIN-93G with 300 ppm added zinc). The diets were procured from Research Diets, Inc., and the diets were irradiated with 10 to 20 kGy gamma-irradiation. All the three groups of mice were provided with sterile food and water ad-libitum, and the mice were maintained on diets for a period of 40 days. The composition details of the diets are provided in Supplementary Table 2.

Acne mouse model

To model human acne disease, mice that were on different diets were intradermally (i.d.) infected with a C. acnes strain (2x10⁷CFU in 50μl volume of BHI media), followed by the topical application of synthetic sebum daily as described previously (Kolar et al. 2019). The i.d. infections were performed under vaporized Isoflurane (Fluriso, Vet One) anesthesia. On the third day post-infection, the disease score was evaluated, and euthanasia was carried out using CO_2. Skin lesions were aseptically excised and harvested in 400 µl sterile phosphate buffer saline (PBS, pH 7.4). The skin lesions were then homogenized and 25 μl was serially diluted (10-fold) in PBS to determine CFU on BHI agar plates. The BHI agar plates were incubated anaerobically at 37 °C for a period of 3-4 days. In addition, the homogenized skin lesions were centrifuged at maximum speed (13000 rpm) for 25 min, yielding a supernatant that was subsequently stored at -80°C for the purpose of cytokine analysis.

Disease scoring

The assessment of gross skin pathology involved scoring according to the following criteria: Erythematous change (none = 0, mild = 1, moderate = 2, and marked = 3); and papule characteristics (flat = 0, small = 1, large = 2, and extra-large = 3).

Determination of cytokines in skin lesions

IL-1β, IL-6, and TNF-α cytokine levels within stored skin homogenates at -80°C were quantified through a solid-phase sandwich ELISA, employing commercially available mouse cytokine ELISA kits (Biolegend, San Diego, CA, USA). For IL-1β and IL-6, skin homogenates (50 μl) were diluted at a 1:1 ratio with blocking buffer (1% BSA plus 1xPBS-Tween20), while undiluted skin homogenate (100 μl) was used for TNF-α. These were assayed alongside cytokine standards of known concentrations (provided with the kits) on each ELISA plate. The plates were developed, and the optical density (OD) was read at 450 nm with wavelength correction set at 570 nm, utilizing a multimode microplate reader (PerkinElmer, Waltham, MA, USA). The cytokine levels in the samples were determined using the standard curve derived from the OD of the cytokine standards.

Statistical analysis

GraphPad prism version 8 was used to analyze all data (GraphPad Software, San Diego, CA, www.graphpad.com). All mouse model data were presented as medians from two independent experiments. For comparisons involving multiple groups, a one-way ANOVA with Tukey’s post-hoc test was employed (HL043PA1 bacteria load and IL-1β). If normality assumptions were not met, non- parametric Kruskal-Wallis one-way ANOVA was utilized to analyze the data (HL043PA1 disease score, IL-6, and TNF-α).