Abstract

Introduction:

Systemic Lupus Erythematosus (SLE) exhibits a pronounced sex bias, affecting females approximately nine times more frequently than males; however, males tend to experience a more severe clinical course yet the molecular basis for these differences remains unclear.

Methods:

Leveraging epigenomic, transcriptomic, and proteomic data from the whole blood of 720 SLE patients (679 females, 41males) and 84 healthy controls (77 females, 7 males), we conducted comprehensive multi-omic analyses to identify sex-specific molecular features of this disease. Specifically, differential expression analysis for each modality was conducted using a factorial design to identify differences between disease and healthy controls (SLE–HC) for each sex, as well as the interaction effects between send and disease ([Male SLE – Male HC] – [Female SLE – Female HC]). Benjamini & Hochberg false discovery rate (FDR) was used for multiple test correction.

Results:

The strongest signal differentiating males and females with SLE was the aberrant expression of the long non-coding RNA, XIST, in males. This XIST expression in males with SLE was bimodal, with 54% of males having elevated XIST expression, and correlated with disease severity. Males with SLE also exhibited significant hypermethylation of the X chromosome and transcriptional silencing of X-linked genes – hallmarks of X-chromosome inactivation (XCI), a process typically restricted to females.

Conclusion:

These results suggest that X-chromosome silencing by XIST may contribute to SLE disease in males.

1 Introduction

Systemic lupus erythematosus (SLE) is a systemic autoimmune disease with a striking female-bias (9:1) (1, 2), but a more severe clinical course in males (3, 4). Although the mechanisms driving this bias remains unclear, X-chromosome dosage and X-linked immune genes have been strongly implicated in SLE pathogenesis. For example, males with Klinefelter’s syndrome (XXY) have a 14-fold increased prevalence of SLE compared to XY males, while XO females have lower risk of developing SLE (2). Genome-wide association studies have identified specific X-linked genes such as TASL (CXorf21) (2, 5) and IRAK1 (6) associated with increased SLE risk, while gain-of-function mutations in the X-linked gene TLR7 (7, 8) have been linked to overactive toll-like receptor signaling in SLE. Another X-linked gene XIST, which equalizes X-gene dosage between sexes by modulating the epigenetic silencing of one of the two X chromosomes in women (9), has been associated with SLE both for its failure to completely inactivate immune genes which partially escape X-inactivation, like TLR 7 (10), and for its immunostimulatory capacity as a long diuridine-rich RNA (11). Despite these genetic and epigenetic insights, the biological mechanisms underlying the striking sex disparity in SLE incidence and severity remains incompletely understood.

These gaps in understanding SLE heterogeneity inhibit the development of new treatments. Therefore, broader and deeper approaches to understanding the molecular biology of this complex disease are needed. Transcriptomic analyses in SLE have consistently revealed a common hallmark of disease pathogenesis: the interferon (IFN) signature (12–16). Today, methylomic and proteomic approaches have been employed to search for biomarkers, predict response to therapy, and increase our understanding of disease pathogenesis (17–23). Specifically, combined methylation and expression data on IFN-responsive genes has been shown to be more sensitive than individual components as a biomarker, while complement components and cytokines have shown utility for diagnosis and subsetting in proteomics studies (17–23). In this study, we apply multi-omic analyses – including epigenetic, transcriptomic, and proteomic profiling – on 720 SLE patients (679 females, 41 males) and 84 healthy controls (77 females, 7 males) to uncover sex-specific molecular signatures.

We present details of the sexually dimorphic signatures of SLE in patients enrolled in the BRAVE-I (NCT03616912) (24) and BRAVE-II (NCT03616964) (25) clinical trials. Notably, SLE males exhibit a bimodal distribution of XIST expression (22 XIST-high, 19 XIST-low), accompanied by widespread X-chromosome hypermethylation and transcriptional silencing of X-linked genes. These features of X-chromosome inactivation (XCI) were present irrespective of XIST expression levels. We further explore molecular and clinical differences between XIST-high and XIST-low males, showing that XIST-high males may have more severe disease. Taken together, our data highlights sex-based differences in SLE manifestation, guiding further research on the role of XIST and X-inactivation in SLE.

2 Materials and methods2.1 SLE patient cohort

Our cohort consisted of the subset of patients from the BRAVE-I (NCT03616912) (24) and BRAVE-II (NCT03616964) (25) clinical trials who had given consent for all 3 omics testing (epigenomic, transcriptomic, and proteomic). This included 312 patients (293 females, 19 males) from BRAVE-I and 408 patients (386 females, 22 males) from BRAVE-II, all aged 18 years or older, with active SLE receiving stable background therapy. As previously described, these patients had a clinical diagnosis of SLE at least 24 weeks before screening; met at least 4 of 11 revised American College of Rheumatology (ACR) criteria for classification of SLE; were positive for at least one of antinuclear antibody, anti-dsDNA, or anti-Smith; had a total Systemic Lupus Erythematosus Disease Activity Index 2000 (SLEDAI-2K) score of at least 6 at screening and a clinical SLEDAI-2K score of at least 4 at baseline; and had at least one British Isles Lupus Assessment Group (BILAG) A score or two BILAG B scores at screening. The study also included 84 age-and sex-matched healthy controls – 77 females and 7 males. Serum samples (for proteomics) and whole blood samples (for transcriptomics and epigenomics) were collected at the same time in the same blood draw for each participant.

2.2 Proteomics

Baseline SLE serum samples were collected according to the clinical trial protocols. Blood was collected using a tube containing a clotting agent (BD SST 5ml tube). Circulating proteins were assayed using the Olink Explore 3072 panel (https://olink.com/products/olink-explore-3072-384), a proteomics platform that combines an antibody-based immunoassay with a proximity oligonucleotide extension assay and signal detection with next-generation sequencing on the NovaSeq6000 instrument (Ilumina Inc.). All serum samples were run as one batch without bridging samples, and each sample was assayed in singular. A set of internal controls (incubation, extension and amplification) was used to assess overall assay quality. Internal plate controls were used to normalize the data and to calculate the limit of detection (LOD) for each assay. Protein expression levels are given as normalized protein expression (NPX) values, which are arbitrary units used on a log2 scale that were normalized to the plate controls. For further quality control, we applied the following stringent steps to account for the observed technical variation in a subset of Olink assays that resulted in dropouts: 1) only data flagged as “PASS” according to the manufacturer’s standards passing technical QC were included; 2) only analytes with NPX values above the LOD in more than 75% of samples were included. This resulted in a total of 2,434 unique protein assays after QC that were used to explore potential novel protein markers and biological pathways in patients with SLE. Four SLE samples did not meet QC standards due to overall low expression patterns (greater than 4 standard deviations) across all assays and were excluded from the analysis.

2.3 RNA and DNA extraction

Whole blood RNA and DNA was extracted from PAXgene blood RNA tubes. Briefly, frozen tubes were thawed at room temperature and centrifuged according to the manufacturer’s protocol. The resulting pellet was then washed with nuclease free water, resuspended in sterile PBS, and split in half. RNA was extracted from one half of the suspension using the QIAsymphony PAXgene Blood RNA Kit (Qiagen). The other half of the suspension was used for DNA extraction using the Mag-Bind Blood & Tissue DNA HDQ Kit (Omega Bio-Tek) on the Kingfisher Flex (Thermo) instrument. RNA concentration and integrity were assessed using the Ribogreen Assay (Thermo) and a Fragment Analyzer (Agilent), respectively. The concentration and integrity for DNA was evaluated by Picogreen Assay and Agarose Gel electrophoresis.

2.4 RNAseq

Polyadenylated RNA was enriched using NEBNext Magnetic Oligo(dT)25 Beads and prepared into RNA-seq libraries using the NEBNext mRNA Library Prep Reagent Set for Illumina (New England BioLabs, Ipswich, MA, USA). Unique dual-index barcodes were incorporated as specified per manufacturer’s protocol. Library concentration was measured using the Quant-iT PicoGreen dsDNA Reagent (Thermo), and library size distribution and quality were checked using a DNA chip on Caliper Gx Touch (Revity). Library quantification was further refined by qPCR using the KAPA Library Quantification Kit (Roche). Paired-end sequencing (2 x 100 bp, 100 million reads per sample) was carried out on an Illumina NovaSeq 6000 platform (Illumina) following manufacturer’s recommended protocols.

Data was processed using nf-core/rnaseq v3.14.0 of the nf-core collection of workflows (26), utilizing reproducible software environments from the Bioconda (27) and Biocontainers (28) projects. The pipeline was executed with Nextflow v23.10.1 (29). Resulting count data was pre-processed using the edgeR package (30) in R, which included filtering out lowly-expressed genes using edgeR filterByExpr (min.count = 10) and performing quantile-normalized voom transformation.

2.5 Enzymatic methylation

Whole genome methylation libraries were constructed using 200ng gDNA input with NEBNext® Enzymatic Methyl-seq Kit as per manufacturer instruction with minor modifications for automation. Library concentration was measured using the Quant-iT™ PicoGreen® dsDNA Reagent (Thermo), and library size distribution and quality were checked using a DNA chip on Caliper Gx Touch (Revity). Library quantification was further refined by qPCR using the KAPA Library Quantification Kit (Roche). Paired-end sequencing (2 × 150 bp, approximately 30X coverage per sample) was carried out on an Illumina NovaSeq 6000 Plus platform (Illumina) following the manufacturer’s recommended protocols.

Data analysis was performed using nf-core/methylseq version 2.3.0 using the Bismark aligner with the GRCh38 human reference genome (GCA_000001405.15_GRCh38_no_alt fasta file with the hs38d1 decoy assembly). Most samples were deemed high quality with sufficient genome coverage > 20X. One sample was rejected as the sequencing data did not show an expected reduction in GC content observed in all other samples after sequencing, still averaging the pre-converted GC percentage of ~ 40%. The expected reduced GC content likely reflected the failure to enzymatically convert free cysteines (31, 32). Decoy sequences were padded with NN to enable extraction of relevant methylation contexts (CG, CHH vs CHG). EMseq specific trimming was set with the flag em_seq = true. Methylation coverage files from Bismark were imported and individual samples smoothed in R using bsseq BSmooth as the sample size was too large to jointly smooth across all samples. A gtf file was provided to summarize the smoothed methylation signals using the bsseq getMeth command resulting in a summarized methylation signal across the regions provided. Gene promotors were calculated with the GenomicRanges promoters command using default parameters for 2000 bp upstream and 200 bp downstream of TSS. Individual samples were aggregated into two methylation matrices representing gene and promotor methylation separately.

2.6 Chromosome read depth analysis

Mean coverage per chromosome was calculated for each patient using samtools (33), and average read depths were compared across sexes. As expected, females exhibited approximately twice the read depth for the X chromosome relative to males, while autosomal coverage was comparable between sexes (Supplementary Figure 1). Based on this approach, XXY Klinefelter’s males would be expected to show elevated X-chromosome coverage, comparable to females. Two outliers were identified: one male with abnormally high X chromosome read depth and one female with unusually low X-chromosome read depth. Both samples were excluded from all analyses.

2.7 RNA extraction & qRT-PCR

Total RNA was extracted and purified from whole blood collected in PAXgene tubes following manufacturer’s instructions. RNA concentration was determined by NanoDrop Spectrophotometer. First-strand cDNA was synthesized using SuperScript VILO cDNA Synthesis kit (Invitrogen) according to the manufacturer’s protocol. Then qRT-PCR was performed on QuantStudio 7 Flex Real-Time PCR System (Applied Biosystems) using TaqMan Fast Advanced Master Mix (Invitrogen) at 60 °C annealing temperature. Taqman assays used (XIST Hs01079824_m1, TSIX Hs03299334_s1, GAPDH Hs99999905_m1) were all obtained from Invitrogen. XIST and TSIX gene expression were normalized to GAPDH housekeeping gene expression in all samples analyzed in qRT-PCR.

2.8 Differential expression, hypergeometric, and gene set enrichment analysis

Differential expression analysis for each modality was conducted using limma (34) in R. A factorial design was applied to identify differences between disease and healthy controls (SLE-HC) for each sex, as well as the interaction effects between sex and disease ([Male SLE – Male HC] – [Female SLE – Female HC]). Benjamini & Hochberg false discovery rate (FDR) was used for multiple test correction, and hits with FDR ≤ 0.1 and |Fold Change (FC)| > 1.5 were considered significant. Genes identified as significantly differentially expressed in both male and female SLE patients from the RNAseq analysis were subjected to hypergeometric overrepresentation analysis using the enrichPathway function from the clusterProfiler (35) R package, which uses pathway terms from Reactome (36). For the other modalities (Olink and EMseq), genes were ranked by log fold change (LFC), and gene set enrichment analysis (GSEA) was performed using the gsePathway function from clusterProfiler. Significantly enriched pathways were defined by an adjusted p-value (FDR) ≤ 0.1.

2.9 Classification of male SLE patients based on XIST expression

To evaluate the potential bimodal nature of XIST expression in individuals with SLE, Hartigans’ dip test for unimodality was applied to both male and female XIST expression using the diptest package in R. Males with SLE were then systematically classified into XIST-high and XIST-low groups using kernel density estimates of qPCR-based XIST expression using the density function in R with the “SJ” parameter which implements the methods of Sheather & Jones (37) to select the bandwidth using pilot estimation of derivatives. The global minimum (trough) of the distribution was used as a threshold to define the two groups (Supplementary Figure 2A). Three male patients were missing qPCR-based XIST quantification. For these individuals, classification was instead based on XIST expression levels obtained from RNAseq data. Comparative visualization of the XIST-groupings using both qPCR and RNAseq XIST expression data revealed clear separation between XIST-high and XIST-low groups (Supplementary Figure 2B), supporting the robustness of the classification.

2.10 Molecular comparison of XIST-high and XIST-low males

Differential expression analysis was conducted using limma (34) in R to compare: i) XIST-high vs. healthy control males, ii) XIST-low vs. healthy control males, and iii) XIST-high vs. XIST-low males across all modalities – EMseq, RNAseq, and Olink. Multiple test correction was applied using Benjamini & Hochberg FDR, with significant hits defined by an adjusted p-value (FDR) ≤ 0.1 and |FC| > 1.5.

2.11 X-chromosome silencing

We constructed a custom gene set that maps each gene to their respective chromosome. Genes in the pseudoautosomal regions (PAR) were relabeled as “PAR” rather than annotated to the X- and Y- chromosomes. Using this custom gene set, we performed GSEA using the GSEA function in the clusterProfiler (35) R package comparing XIST-high, X-low, and females with SLE to their respective healthy controls across all modalities – EMseq, RNAseq, and Olink (Supplementary Figure 3). This approach enabled the identification of enrichments patterns at the chromosome level, rather that within predefined biological pathways. To further assess whether genes that were significantly hypermethylated (EMseq) or downregulated (RNAseq) in these subgroups were enriched in X-linked genes, we conducted hypergeometric enrichment using this custom gene set and the enricher function in the clusterProfiler (35) R package.

2.12 Clinical characteristics

Pairwise wilcoxon rank sum test of 95 numerical baseline clinical metrics were compared across groups: 1) XIST-high males vs females with SLE, 2) XIST-low males vs females with SLE, and 3) XIST-low vs XIST-high males with SLE. Correlations between XIST expression and the 95 clinical metrics was assessed using spearman correlation in R. Pairwise chi-square test were run to compare 31 categorical clinical metrics across the same 3 groups. In all cases, Benjamini & Hochberg FDR was used to correct for multiple test correction.

2.13 Correlation of XIST expression with type I IFN signature

Type I IFN signature was calculated by aggregating the expression of 6 type I IFN genes: MX1, OAS3, LY6E, USP18, IFI44, DDX60 (38). Spearman correlation was used to assess the association between XIST expression and the type I IFN signature in males and females separately.

2.14 Deconvolution of immune cell proportions from bulk RNA-seq

In-silico deconvolution of bulk RNA-seq data to estimate immune cell composition was performed using the granulator package in R (39). Raw gene counts were normalized to TPM. We evaluated 4 reference profiles (ABIS_S0-S3) (40) in combination with 7 deconvolution algorithms: dtangle (41), non-negative least squares regression model (nnls), ordinary least squares (ols), quadratic programming without constraints (qprog), quadratic programming non-negative and sun-to-one constraints (qprogwc), robust linear regression (rls), and support vector regression (svr). Notably, dtangle failed to generate cell proportions estimates for 2 of the 4 references profiles. This resulted in 26 unique deconvolution sets, which were subsequently benchmarked against immune cell proportions measured via flow cytometry (Supplementary Figures 4A-G). We selected the reference-method pair that achieved high concordance (Pearson correlation) with B-cell and NK-cell proportions from flow – ABIS_S0 with dtangle, resulting in the estimated composition of 17 immune cell types: naïve B-cells (B.Naive), memory B-cells (B.Memory), plasmablasts (Plasmablasts), memory helper T-cells (T.CD4.Memory), naïve helper T-cells (T.CD4.Naive), memory cytotoxic T-cells (T.CD8.memory), naïve cytotoxic T-cells (T.CD8.Naive). non-Vδ2 γδ T-cells (T.gd.non.Vd2), Vδ2 γδ T-cells (T.gd.Vd2), mucosal-associated invariant T-cells (MAIT),classical monocytes (Monocytes.C), non-classical & intermediate monocytes (Monocytes.NC.I), natural killer cells (NK), low-density neutrophils (Neutrophils.LD), low-density basophils (Basophils.LD),myeloid dendritic cells (mDCs), plasmacytoid dendritic cells (pDCs). Pairwise Wilcoxon rank sum test on the estimated cell types were compared across groups: 1) SLE females vs. female healthy controls, 2) XIST-high SLE males vs. male healthy controls, 3) XIST-low SLE males vs. male healthy controls, and 4) XIST-high vs. XIST-low SLE males. Additionally, Spearman correlations were computed between XIST expression levels with inferred cell type proportions.

3 Results3.1 Patient cohort

The study cohort comprised 720 patients with active SLE, including 312 participants (293 females, 19 males) from the BRAVE-I trial and 408 participants (386 females, 22 males) from the BRAVE-II trial. All individuals were aged 18 years or older and receiving stable background therapy at the time of enrollment. Patients had a confirmed clinical diagnosis of SLE and were seropositive for at least one autoantibody—antinuclear antibody, anti-dsDNA, or anti-Smith at screening. Disease activity thresholds included a total SLEDAI-2K score ≥6 at screening, a clinical SLEDAI-2K score ≥4 at baseline, and either one BILAG A score or two BILAG B scores. Summary clinical statistics for each patient cohort are shown in Table 1. Additionally, 84 age- and sex-matched healthy controls were included for comparative analyses –. All participants underwent baseline sample collection for subsequent multi-omic profiling.

Clinical variableBRAVE-IBRAVE-IIF (n=293)M (n=19)F (n=386)M (n=22)Mean age, years41.8436.6342.2438.50Mean time since onset of SLE, years9.287.658.456.08SLEDAI-2K score9.939.6810.1010.09CLASI score5.675.796.985.41SLICC score0.660.580.570.32PGA score60.1663.7459.4857.32MedicationAntimalarial2551831919Azathioprine492612Corticosteroid2111431820Hydroxychloroquine2301727216Immunosuppressants1731121011Methotrexate746845Non-steroidal anti-inflammatory drug903953RaceAmerican Indian or Alaska Native161211Asian23312310Black or African American544281Multiple1060White1991120710Native Hawaiian or Other Pacific Islander0010

Baseline demographics and clinical characteristics.

3.2 Proteomic, transcriptomic, and methylomic differences in SLE by sex

To elucidate the molecular differences between males and females with SLE, we analyzed baseline methylome (EMseq), transcriptome (RNAseq), and proteome (Olink) data collected from whole blood from the cohort of 720 SLE patients and 84 healthy controls. Leveraging a factorial design, we performed 3 differential expression analyses for each modality: (i) SLE vs. healthy controls in females, (ii) SLE vs. healthy controls in males, and (iii) interaction between sex and disease (Figure 1A). The interaction analysis accounts for baseline sex differences and isolates molecular features that differ in disease response between males and females. We first focused on the sex-stratified comparisons of SLE versus healthy controls (i and ii) to characterize disease-associated molecular signatures within each sex.

Multi-omics profiling of SLE reveals sexual dimorphism at all levels of genetic regulation. (A) Overview of study design. A total of 806 individuals—720 patients with SLE and 84 healthy controls—were enrolled from the BRAVE-I and BRAVE-II clinical trials. Baseline multi-omics profiling was performed using Olink proteomics, RNA sequencing (RNA-seq), and enzymatic methyl-sequencing (EM-seq). Sex-dependent molecular differences were assessed using a 2×2 factorial interaction model across each modality. (B-E) Venn diagrams depicting the number of significantly altered molecular features associated with disease (FDR ≤ 0.1 and |fold change| > 1.5) across four omics layers: promoter methylation (B), gene methylation (C), RNAseq (D), and Olink (E). Circle sizes are proportional to the number of significant features identified in each modality. For each Venn diagram, the left (pink) circle represents female-specific changes, and the right (blue) circle represents male-specific changes. Within each circle, upregulated features are shown in darker shades at the top, while downregulated features are shown in lighter shades at the bottom. (F) Volcano plot displaying all molecular features with significant sex:disease interaction effects across modalities. The x-axis represents the magnitude of the interaction effect, while the y-axis indicates statistical significance. Data points are color-coded by modality: promoter methylation (light blue), gene methylation (dark blue), Olink proteins (green), and RNA-seq transcripts (maroon).

3.2.1 Epigenetic profiling

Males with SLE exhibited significantly more alterations in the methylome compared to females. Of the 36,175 promoter and 34,981 gene methylations quantified using EMseq, we observed 104 and 110 significant changes (|FC|>1.5, FDR<=0.1) in promoter and gene methylation relative to healthy controls, respectively, in males, compared to just 15 and 25 in females (Figures 1B, C). Most of these changes represented hypermethylation of DNA, associated with epigenetic silencing. Only a single promoter (IFI44L) and 2 genes (IFI44L and IFI44) exhibited shared methylation changes in both males and females, all of which were hypomethylated, suggesting epigenetic activation. Among these shared loci, males demonstrated a greater magnitude of change compared to females. Most epigenetic changes in males with SLE were to promoters and genes on the X chromosome (84/104 promoters and 80/110 genes). This was not the case for females, in which only 1/15 promoter changes and 3/25 gene changes were X-linked. All statistically significant changes in methylation for both sexes can be found in Supplementary Tables 1 and 2A, B.

GSEA analysis of the methylomes (Supplementary Table 1 and 2C, D) show males with SLE exhibit reduced methylation, indicative of epigenetic activation, in pathways associated with mitophagy, autophagy, protein ubiquitination, and TNFR1 signaling compared to male healthy controls. Males with SLE also exhibited increased methylation, indicative of epigenetic silencing, in olfactory pathways dominated by olfactory receptor (OR) family genes, which are thought to be involved in modulating inflammatory macrophage activity (42, 43). Females on the other hand exhibited no significant epigenetic enrichments across any pathways, consistent with the relatively limited epigenetic changes observed in females compared to males.

3.2.2 Transcriptomic profiling

Males and females with SLE exhibited comparable numbers of differentially expressed transcripts relative to healthy controls, with 2,735 in females and 2,014 in males, out of 22,443 transcripts profiled by RNAseq (Supplementary Tables 3A, B). Of these differentially expressed transcripts, 1,080 were shared between males and females (Figure 1D). Hypergeometric enrichment analysis on this set of genes commonly dysregulated in both males and females with SLE revealed that the genes upregulated in both sexes were enriched for pathways related to interferon signaling (e.g. interferon alpha/beta signaling, interferon gamma signaling, ISG15 antiviral mechanism, and regulation of INFA/INFB signaling), complement cascade (e.g. initial triggering of complement, regulation of complement cascade, and creation of C4 and C2 activators), and cell cycle/mitosis (e.g. cell cycle checkpoints, DNA replication pre-initiation, mitotic G1 phase and G1/S transition, and separation of sister chromatids). Genes commonly downregulated in both sexes showed no significant enriched pathways. A comprehensive list of enriched pathways for commonly dysregulated genes is provided in Supplementary Tables 3C, D.

3.2.3 Proteomic profiling

Among the 2,434 proteins measured in the Olink panel, 259 showed significant changes in females (|FC|>1.5, FDR<=0.1), compared to just 92 in males (Figure 1E). We observed 33 proteins commonly upregulated in both sexes, including several well-characterized cytokines and chemokines. Only 3 proteins—KIT, CD63, and RNASE3—were consistently downregulated in both males and females. GSEA of the proteomes revealed shared upregulation of several immune-related pathways in both sexes, such as signaling by interleukins, cytokine signaling in immune system, and interleukin-10 signaling. In females, additional positive enrichment was observed in chemokine and interleukein-12 signaling; and negative enrichment of pathways involved in Rho GTPase signaling, tyrosine kinase signaling, phagocytosis, and death receptor signaling. The full lists of proteins with altered abundances in each sex is given in Supplementary Tables 4A, B, and all GSEA pathways are given in Supplementary Tables 4C, D.

3.3 XIST expression emerges as the dominant feature in sex-specific SLE

To investigate sex-specific molecular responses in SLE, we examined the interaction between sex and disease in our differential expression analyses. This approach identifies molecular features whose association with SLE differs between males and females, accounting for the natural differences in sex. From the EMSeq data, we identified 73 promoters and 67 genes exhibiting sex-dependent disease-associated methylation changes. Only three transcripts – XIST, ZNF630, and NEURL1-AS1 – and 11 proteins showed significant sex-disease interaction effects (Figure 1F). Complete lists of promoter methylation, gene methylation, transcripts, and proteins with significant sex-disease interaction effects are provided in Supplementary Tables 5A–D. Notably, expression of the long non-coding RNA XIST emerged as the strongest and most statistically significant sex:disease interaction across all modalities. While XIST is a well characterized female-specific transcript with naturally variable expression between males and females, baseline sex differences have been accounted for in this model, and its persistent emergence therefore indicates an unexpected disease-specific divergence in molecular response between the sexes. In short, XIST is upregulated in SLE, but only in males.

3.4 Bimodal expression of XIST in males with SLE

Focusing on the strongest molecular signal differentiating males and females with SLE – XIST expression – we observed that females with SLE showed no significant changes in XIST expression relative to healthy female controls (LFC = 0.12, FDR = 0.44), while males with SLE exhibited a dramatic upregulation compared to healthy male controls (LFC = 4.71, FDR = 3.9e-18; Figure 2A). Further visualization of XIST expression in males revealed a bimodal distribution. To confirm this, we applied Hartigans’ dip test for unimodality to XIST expression in males with SLE. This yielded a p-value of 0.0027, indicating a statistically significant deviation from unimodality, consistent with a bimodal distribution of XIST expression in males. Analysis of XIST expression in females with SLE resulted in a dip test p-value of 0.99, suggesting XIST expression is unimodal in women.

XIST expression and X-linked molecular changes in males with SLE. (A) Violin plots showing XIST expression as measured by RNA sequencing in females (left) and males (right) with SLE vs. healthy controls. (B) Scatter plot showing the correlation between XIST expression measured via RNAseq (x-axis) and qPCR (y-axis). (C) Boxplots of XIST expression measured via qPCR. Values here have been normalized based on expression in sex-matched healthy controls. (D-G) Faceted barplots illustrating the proportion of X-linked changes among all significantly hyper- and hypo-methylated promoters (D, E) and significantly up- and down-regulated genes and proteins (F, G). (H) Volcano plot showing X-chromosome enrichment across the proteome (green), trancriptome (maroon), gene methylome (dark blue), and promoter methylome (light blue); the dashed grey line indicates the significance threshold (FDR ≤ 0.1). Comparisons include SLE females vs. female healthy controls (square), SLE XIST-high males vs. male healthy controls (triangles), SLE XIST-low males vs. male healthy controls (astericks), and SLE XIST-high males vs. XIST-low males (diamonds). Negtive normalized enrichment scores (NES) indicate enrichned downregulation of X-linked gene changes, reflecting X-linked silencing in the proteome and transcriptome and activation in the epigenome, and vice versa. Scores for all other chromosomes can be found inSupplementary Table 6. (I) Cnetplot of differentially expressed X-linked genes in females, XIST-high males, and XIST-low males compared to healthy controls. Genes commonly dysregulated across all groups are highlighted and shown in the accompanying heatmap. Note: XIST was excluded from the plot to enhance visualization due to its dominant signal.

To rule out any assay-related artifacts that could be underlying the unexpected XIST expression in males, we validated these findings via RT-PCR for XIST using all males in the cohort (41 SLE, 7 healthy controls) as well as a small set of females (4 SLE, 5 healthy controls) for comparison. The RT-PCR results closely mirrored the RNAseq data (r = 0.99, p < 2.2E-16) (Figure 2B), confirming the highly unusual expression of XIST in males with SLE. Importantly, RT-PCR confirmed the bimodal nature of XIST expression among males with SLE (Figure 2C), with a dip test p-value of 0.0008. We subsequently classified males with SLE into XIST-high (n=22) and XIST-low (n=19) phenotypes using kernel density estimates (Table 2).

Clinical variableXIST-high (n=22)XIST-low (n=19)Mean age, years38.2336.95Mean time since onset of SLE, years4.729.22SLEDAI-2K score10.689.00CLASI score5.685.47SLICC score0.410.47PGA score58.7362.1119MedicationAntimalarial1918Azathioprine40Corticosteroid2014Hydroxychloroquine1716Immunosuppressants1210Methotrexate47Non-steroidal anti-inflammatory drug42Mean BILAG scoresConstitutional1.551.42Hematological2.272.21