Introduction

The cutaneous manifestations associated with COVID-19 vary widely and can be broadly classified into two main categories: inflammatory reactions (eg, morbilliform (measles-like) rash, vascular rashes and vesicular rash) and vascular-origin lesions (eg, chilblain-like rashes, petechiae/purpura, and livedo racemosa-like patterns).^1,2 Vesicular eruptions appear early in the course of the disease (15% before other symptoms). The pseudo‐chilblain pattern frequently appears late in the evolution of the COVID‐19 disease (59% after other symptoms), while the rest tend to appear with other symptoms of COVID‐19.³ Despite numerous reports of diverse cutaneous manifestations linked to SARS-CoV-2 infection, dermatologists treating COVID-19 patients at the University of São Paulo Medical School Hospital—a major reference center and one of the largest university hospitals in Latin America confirmed that was not possible to attribute with certainty that any cutaneous lesions had a direct correlation to the SARS-CoV-2 infection in 86 patients confirmed with COVID-19.⁴ Instead, the dermatological presentations observed mirrored those commonly encountered prior to the pandemic, with fungal infections being the most frequent, followed by drug eruptions and viral infections. This underscores the need for further molecular investigations to unravel the immunological changes occurring in the skin during COVID-19. In same line, it has been verified a lack of association between chilblains and SARS-CoV-2 infection.⁵

SARS-CoV-2 enters host cells via the angiotensin-converting enzyme 2 (ACE2) receptor and TMPRSS2 protease, which are expressed in multiple tissues, including the lungs, kidneys, colon, heart, testicles, and skin.⁶ Despite this, studies have found no virally induced cytopathic alterations or intranuclear inclusions in the skin,⁷ nor detectable SARS-CoV-2 viral load in skin biopsy specimens.^8,9 Nevertheless, upon viral entry, after the latent period, depending on the variant, such as Delta, the latency period was 3.2 days.¹⁰ After this, a cascade of activation of innate immunity is rapidly activated, including the recruitment of dendritic cells, macrophages, and monocytes. These immune cells adhere to infected cells and release type I interferons (IFNs), which act as danger signals, along with pro-inflammatory and antiviral cytokines.

Understanding the molecular immunological events in the skin during SARS-CoV-2 infection requires a focus on innate immunity, which involves pattern recognition receptors (PRRs), Toll-like receptor (TLR) signaling, and cytokine-mediated antiviral responses. To determine the extent to which the skin is activated in response to SARS-CoV-2. If activation does occur, it is crucial to identify the controlling factors.

In this study, we conducted a transcriptomic analysis of skin samples from patients with acute COVID-19 as depicted in Figure 1A. Our findings reveal an integrated gene expression profile encompassing antiviral, antioxidant, and innate immune responses in the skin, shedding light on the complex host-pathogen interactions in this compartment.

Figure 1 (A) Schematic workflow of the study; (B) Unsupervised z-score heatmap of normalized gene expression counts, based on 150 down and 150 upregulated genes. The control group (CT) is represented in green, and the COVID-19 group (CV) in pink. Red shades indicate z-scores above the row average (higher expression), while blue shades indicate z-scores below the row average (lower expression).

Materials and Methods

Patients’ Enrollment

Skin samples used in this study were obtained from patients who succumbed to severe SARS-CoV-2 infection while hospitalized in the Intensive Care Unit of Hospital das Clínicas, São Paulo. Post-mortem biopsies were collected in the thigh region using a 5 mm punch, by the Death Verification Service. All patients were diagnosed with COVID-19 through nasopharyngeal detection of SARS-CoV-2 RNA using reverse transcriptase polymerase chain reaction (RT-PCR).

For RNA sequencing, we included skin samples from eight COVID-19 patients (4 men and 4 women), aged 32–70 years, who died between April and June 2020 (Supplementary Table 1). Additional skin samples from COVID-19 patients (n=10, 3 men and 7 women, aged 26–76 years) were used for real-time PCR analysis. At the time, vaccines against SARS-CoV-2 were unavailable. Post-mortem skin biopsies were collected within 20 hours of death from the same anatomical region for all patients. Post-mortem skin biopsies were taken up to 20 hours after death due to complications arising from a diagnosis of SARS-CoV-2. The exclusion criteria were cancer, HIV-1 infection, and others skin diseases.

As controls, skin biopsies were obtained from healthy donors (n=5) prior to the COVID-19 pandemic. These control samples were collected from Caucasian individuals of both sexes, aged 35–70 years. The autopsy was carried out with the authorization of the relatives, who signed an informed consent form approved by the National Research Ethics Committee (CONEP). The study was approved by the ethics committee of the HC-FMUSP (CAAE 30364720.0.0000.0068) and conducted following the Declaration of Helsinki.

Transcriptomic Analysis

Total RNA was extracted from autopsy-derived skin specimens using the RNeasy Plus Mini Kit (Qiagen, Hilden, Germany), following the manufacturer’s instructions. Tissue homogenization was performed with a TissueRuptor (Qiagen). RNA integrity was assessed using the TapeStation system (Agilent, Santa Clara, CA, USA).

Transcriptomic profiling was carried out using the Illumina TruSeq Stranded mRNA Library Prep Kit (Illumina, San Diego, CA, USA) at the FMUSP Large-Scale Sequencing Laboratory (SELA). Sequencing adapters were removed, and quality control of raw reads was performed using FastQC.¹¹ Reads were aligned to the reference human genome (GRCh38.p13, version 107) using the SUBREAD package¹² and quantified with FEATURECOUNTS.¹³

Differentially expressed genes (DEGs) were identified using the DESeq2 package from R/Bioconductor,¹⁴ with statistical significance determined at an adjusted p-value ≤ 0.05 (Benjamini-Hochberg correction). Functional enrichment analyses were conducted with the FGSEA package¹⁵ using KEGG¹⁶ and GO Biological Process¹⁷ databases. Protein-protein interaction networks were constructed using Cytoscape¹⁸ with data from the InnateDB database.¹⁹

Gene Expression by Real-Time PCR

Total RNA was quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific). Complementary DNA (cDNA) was synthesized using the iScript™ Reverse Transcription Kit (Bio-Rad, Hercules, CA, USA). Real-time PCR was performed with SYBR Green (Applied Biosystems, Waltham, MA, USA) and primers specific to the target genes (primer sequences listed in Supplementary Table 2). Amplification was carried out on a 7500 Real-Time PCR System (Applied Biosystems), and data analysis was conducted using 7500 software v2.0.6 (Applied Biosystems) based on the ΔΔCt method.²⁰

Statistical Analysis

Statistical comparisons between groups were performed using the Mann–Whitney test, with significance set at p ≤ 0.05. Detection of outlier samples, normalization, statistical tests and identification of DEGs were carried out using the package R DESEQ2 from the Bioconductor repository. The p-values were corrected at using the Benjamini-Hochberg method. Genes were considered differentially expressed for p-adj values ≤ 0.05 and absolute fold-change values, FC ≥ 1.5.

Results

Characteristics of the COVID-19 Patients and Healthy Controls

A group of 8 patients who died due to COVID-19 complications (CV) from April to June 2020, composed of 4 men and 4 women, aged between 32 to 76 years (mean = 55,25), was evaluated (Supplementary Table 1). All cases configure severity and death, and comorbidities were identified in almost all individuals. As a control group, we used skin samples from 5 healthy donors (HD), composed of 2 men and 3 women, aged between 29 and 70 years (mean = 62). HD skin samples were collected in 2019, before the COVID-19 pandemic. No comorbidities or infectious diseases were present in this group.

Transcriptomic Profile of Skin Samples

Skin samples were processed for RNA extraction, followed by library construction and sequencing. Computational analyses were performed to identify differentially expressed genes (DEGs), enriched molecular pathways, and gene interaction networks. Gene expression of the most important DEGs were validated using real-time PCR (Figure 1A).

Unsupervised z-score heatmap of normalized counts revealed a distinct transcriptomic profile between COVID-19 (CV) and healthy control (CT) skin samples (Figure 1B). This showed 150 upregulated and 150 downregulated genes. The more important genes related to innate immunity were labelled. A total list of genes can be found in the Supplementary Table 3. A volcano plot demonstrated with total genes, that 1188 upregulated and 803 downregulated genes in the CV group compared to the CT group (Figure 2). Principal component analysis (PCA) showed that the great majority of the samples grouped closely within their respective cohorts, except for one perhaps mild CV sample (Figure 2).

Figure 2 Principal component analysis (PCA) showing the explained variance projected onto the two largest orthogonal principal components; Volcano plot depicting the -log10(adjusted p-value) as a function of log2 fold change (log2FC). Upregulated genes are shown as red circles, downregulated genes as blue circles, and non-significant genes as grey circles. Vertical dashed lines denote fold change thresholds (absolute log2FC > 1.5), while the horizontal dashed line represents the adjusted p-value significance cutoff (p ≤ 0.05).

Enrichment molecular pathway analysis revealed upregulation of 25 pathways associated with inflammatory responses in the CV group compared to control group. These pathways included “regulation of inflammatory response”, “response to molecules of bacterial origin”, “leukocyte chemotaxis”, “signaling by interleukins”, “response to oxidative stress”, “matrix organization”, and the “JAK-STAT signaling pathway”, all of which are linked to innate immunity and inflammation (Figure 3A). Moreover, pathways analyzed included response to oxidative stress as “reactive oxygen species metabolic response” and “response to oxidative stress”. This finding was in line with downregulation of pathways related to controlling of inflammatory process, such “positive regulation of signaling”, “regulation of cell differentiation”, “negative regulation of response to stimulus”, “homeostatic process” (Figure 3B).

Figure 3 (A) Enrichment analysis of the MSIG molecular pathways database. The normalized enrichment score (NES) is displayed at the center of each circle, with shades of red indicating upregulated pathways. The circle size represents the number of genes in each pathway. (B) Over-representation analysis obtained with the downregulated genes against the MSIG molecular pathways database. The percentage of overlapping genes is shown at the center of each circle, with colors reflecting the degree of overlap. The circle size corresponds to the number of overlapping genes.

Over-representation analysis revealed downregulated molecular pathways in CV samples included those related to the skin matrix, such as “matrisome”, “collagen-containing extracellular matrix”, “epidermis development”, “regulation of cell differentiation”, and “connective tissue development” (Figure 3B).

The gene interaction network revealed several upregulated DEGs (pink nodes) associated with antiviral responses, including interferon-stimulated genes (ISGs) such as IFIT1, IFIT2, IFIT3, and interferon-regulatory factors IRF3 and IRF7. Additional antiviral genes included members of the S100 protein family, such as S100P, S100A7, S100A8, and S100A9. Cytokine and cytokine-related factors were also prominent, including STAT2, STAT3, SOCS3, NF-κB2, RELB, NFKBIA, MyD88, IL-6, CCL20, and CXCL2 (Figure 4). Some factors were downregulated (blue nodes), such as TNF, a multifunctional proinflammatory cytokine, and MAVS (mitochondrial antiviral signaling protein), which regulates IFN-β expression and contributes to antiviral innate immunity (Figure 4). Figure 5 shows the selected DEGs for validation in skin samples, which showed upregulation of S100P, GPX3, SOD2, SOCS3, SERPINE1, and CXCL2, alongside downregulation of TNF-α in the skin of COVID-19 patients compared to the CT group.

Figure 4 Network of Genes Involved in COVID-19 Skin Infection Pathways. The network visualization illustrates genes associated with COVID-19-related pathways in the skin. Each circle represents a gene, with the color intensity indicating the log2FC as shades of red represent upregulated genes (FC ≥ 1.5), and shades of blue represent downregulated genes (FC ≤ −1.5). Only genes with an adjusted p-value < 0.05 were included. Genes directly involved in COVID-19-specific pathways are highlighted with magenta borders. The size of each circle corresponds to the number of connections (degree) the gene has within the network, indicating its relative importance.

Figure 5 Selecting innate immune factors to validated genes in the skin of COVID-19 patients via qRT-PCR. Gene expression levels of S100P, S100A8, S100A9, MMP14, CXCL2, CCL20, TNF-α, SERPINE1, GPX3, SOD2, SOCS3, and ADAMTS1 were assessed in skin autopsies. Data are presented as medians. Statistical significance was determined using the Mann–Whitney U-test, with significance levels indicated as p < 0.05 (*), p < 0.01 (**), and p < 0.001 (***). Green dots represent healthy controls (CT), while pink dots represent COVID-19 patients (CV).

In addition, deconvolution analysis of immune cell types revealed an increased presence of macrophages and reduced dendritic cells in the skin tissue of COVID-19 patients, likely reflecting their migratory status (Supplementary Figure 1). Eosinophils were only present inside the dermal perivascular vessel, not in the tissue infiltration. There is mild lymphohistiocytic infiltration in skin tissue (see Supplementary Figure 2).

Discussion

The gene expression profile of the skin in COVID-19 patients revealed significant upregulation of innate immune factors including antiviral components, antioxidant defenses, and inflammatory markers. Although the skin is not a primary target organ for SARS-CoV-2, it engages protective mechanisms that may prevent severe cutaneous lesions. Network analysis highlighted the involvement of innate antiviral immune factors such as Mx1, APOBEC3G, IRF3, IRF7, IRF9, DDX58, IFIT3, and IFIT1. These genes are associated with antiviral responses triggered by interferon signaling, indicating a robust cutaneous antiviral response. Additionally, interferon-induced transmembrane (IFITM) proteins, known to be upregulated in multiple tissues such as the lung, gut, heart, and brain,²¹ were also expressed in the skin, further underscoring its protective role. Despite the absence of detectable viral load in most cases (present in only two samples), all skin samples exhibited upregulation of antiviral genes.

In terms of innate immunity, the differentially expressed genes found in COVID-19 skin samples included members of the S100 protein family, such as S100A8, S100A9, and S100A7. These proteins, particularly S100A8/A9, are known alarmins that mediate host pro-inflammatory responses and have been identified as biomarkers of COVID-19 severity.²² Consistent with findings in nasal swabs, the upregulation of S100 genes was associated with disease severity.²³ In the skin, S100A8/A9 acts as a ligand for TLR4 and RAGE, promoting NF-κB signaling and inducing pro-inflammatory cytokines such as IL-6, which was notably upregulated in COVID-19 skin samples. This contrasts with the downregulation of TNF-α, a cytokine often implicated in the cytokine storm of severe COVID-19 cases.²⁴ The reduced expression of TNF-α in the skin may represent a counter-regulatory mechanism to limit excessive inflammation, possibly mediated by the upregulation of SOCS3, a known negative feedback inhibitor of cytokine signaling, and also includes NFKBIA which encodes a protein known as IκBα (NF-κB alpha inhibitor). These findings align with prior evidence of cross-regulation between TNF-α and type I IFNs triggering,²⁵ where TNF-α has been shown to suppress type I IFN production and pDC activity, while IFN-β can inhibit TNF-α secretion. Given that TNF-α may suppress an important antiviral factor such as type I IFN, the cutaneous reduction in TNF-α is timely.

Over-representation analysis obtained with the downregulated genes related with impaired development of the extracellular epidermis and connective tissue and influence the dynamic repair of skin tissue in skin of severe cases of COVID-19. However, the upregulation of CXCL2 plays a role in the stages of haemostasis, inflammation, proliferation and remodeling in the skin. This finding suggests that, rather than acting as an eosinophil/neutrophil chemoattractant, CXCL2 acts more in skin remodeling. This is particularly relevant given that skin infiltration in samples from patients with SARS-CoV-2 was primarily composed of lymphocytes and histiocytes (macrophages).

The number of downregulated genes related to innate immunity in the skin was much smaller than the number of upregulated genes. Such as, TNF-α, mitochondrial antiviral signaling protein (MAVS) and WNT2. MAVS plays a key role in the antiviral response. MAVS is activated in response to viral infection by RIG-I-like receptors that recognize viral RNA, mediating the activation of NF-κB and interferon-regulatory factors, as well as the induction of interferons, in response to viral infection. It is required for interferon induction in dendritic cells (DCs).²⁶ This observation aligns with previous reports that SARS-CoV-2 nucleocapsid protein inhibits type I IFN responses by interacting with MAVS, thereby suppressing host antiviral defenses. However, it remains unclear whether the downregulation of MAVS in the skin reflects a specific viral evasion strategy. The Wnt/β-catenin signaling pathway inhibits the replication of SARS-CoV-2 and other pathogenic RNA viruses in vitro.²⁷ Furthermore, the Wnt/β-catenin signaling pathway strongly suppresses peroxisome biogenesis. These are intracellular organelles containing enzymes for beta-oxidation, and they are depleted from cells and tissues during SARS-CoV-2 infection. This depletion may be expected to impair the IFN response.²⁸

Despite the downregulation of MAVS, our findings showed that some first-line defense antioxidants, such as superoxide dismutase (SOD) and glutathione peroxidase (GPX), are upregulated in the skin. These genes play critical roles in mitigating oxidative stress, with SOD2 converting reactive oxygen species (ROS) to hydrogen peroxide and GPX3 metabolizing hydrogen peroxide and lipid hydroperoxides toxins.²⁹ This robust antioxidant response highlights the importance of oxidative stress mitigation in SARS-CoV-2 pathogenesis and suggests potential therapeutic targets.

Another notable upregulated DEG was SERPINE1, which encodes plasminogen activator inhibitor 1 (PAI-1). Upregulation of PAI-1 may counteract the thrombosis frequently observed in COVID-19 by modulating fibrinolysis. SERPINE 1 is required to stimulate keratinocyte migration during the repair of cutaneous injuries.³⁰ A polymorphism in the human SERPINE1 gene increases influenza A virus susceptibility in vitro.³¹ Furthermore, the ability of SERPINE1 to inhibit TMPRSS2 has been demonstrated, resulting in a reduction in SARS-CoV-2 entry and infection.³² We observed that protease inhibitors such as SERPINE 1, was observed in heatmap as DEGs, in the expressions in skin from COVID-19, have shown to be valuable component of the antiviral innate immunity.

Analysis of immune cell composition revealed an increased presence of macrophages in the skin of COVID-19 patients, along with a reduction in DCs. The decreased number of DCs in the skin may reflect their migration to the draining lymph nodes, where antigen-specific adaptive immunity occurs. Conversely, eosinophils were only present in the perivascular dermal vessels and not in the tissue infiltrate, which was found to be mild lymphohistiocytic infiltration of the skin.

Conclusion

SARS-CoV-2 triggers a significant innate immune response in the skin, characterised by the upregulation of antiviral, inflammatory and antioxidant pathways. Several genes related to innate immunity are upregulated, including members of the S100 family, IFITs, ISGs, pro-inflammatory cytokines, and chemokines, indicating a robust cutaneous antiviral response but also with inflammatory potential. In this line, a counter-regulatory response is induced, by genes such as SOCS3, which controls cytokine levels, and NFKBIA, which controls NF-κB activity. These genes may, partly, control the inflammatory status in the skin. The severe COVID-19 patient showed intense innate immunity enriched with antiviral factors. This could explain why they did not exhibit significant cutaneous lesions. Further investigation into the interaction of immune pathways in the skin could lead to the development of new therapeutic strategies for controlling the inflammation associated with the disease.