T and B cells recognize splicing-derived neutrophil neoantigens in lupus

Patients with SLE

We obtained freshly drawn blood from adult (n = 34) patients with SLE, healthy individuals (n = 20), and patients with rheumatoid arthritis (n = 15) to serve as disease controls. Study participants were recruited through the University of Washington, Division of Rheumatology Biorepository to participate in research studies at the University of Washington Medical Center, Seattle, WA. The study was approved by the institutional ethics board (STUDY00006196), the research was performed in accordance with all relevant guidelines and regulations, and informed written consent was obtained from all participants according to the Declaration of Helsinki.

The 34 patients with SLE included in our study represent the typical spectrum of SLE disease seen in an academic medical centre with standard-of-care medication (Supplementary Table 1). Eleven patients had active disease, defined as having SLEDAI ≥ 6, while 24 had inactive disease with a SLEDAI < 6 at the time of venipuncture. Based on the 95th percentile of the expression of 111 interferon-inducible genes in healthy donors, 26 were considered interferon signature positive. All healthy donors and patients with rheumatoid arthritis were interferon signature negative.

Isolation of neutrophils and peripheral blood mononuclear cells

Polymorphonuclear (PMN) cells were isolated from freshly drawn venous blood (within 1 h of venipuncture) by gradient centrifugation on PolymorphPrep according to the manufacturer’s instructions. The PMN fraction was >95% CD66b⁺ and >85% CD15⁺. The cells were washed and suspended in Hanks’ buffered salt solution at 10⁷ per ml.

RNA isolation and RNA-seq

Extraction of RNA was performed using Quick-RNA Whole Blood kit (Zymo Research, R1201). Cells were lysed in DNA/RNA Shield and RNAs were isolated following manufacturer’s instruction and eluted in RNase-free water. RNA was sent to the UW Northwest Genomics core for strand-specific paired-end RNA-seq at 50 million 101-bp reads per sample.

Alignment to reference genome and alternative splicing detection

After quality control of sequencing quality, FASTQ files were aligned to the reference genome (GRCh38) using STAR (v2.7.9a)⁶¹ with two-pass mode which improves detection of reads mapping to novel splicing junctions. After post-alignment quality control, sorted and indexed BAM files were used to identify alternative splicing events by SplAdder (v3.0.3)⁶² which uses alignment and annotation data to create splicing graphs that are enriched with splicing events based on the provided RNA-seq data, followed by quantification of splicing events. The software filters input alignment reads based on their quality and allows four levels of confidence; we used the highest confidence setting. Additionally, this tool uses event-specific criteria to confirm splicing events based on genomic coordinates and coverage; we only included confirmed events in our analysis. Both annotated and novel splicing events were detected including six types of events: exon skip, alternative 3′ or 5′ splice site, intron retention, mutually exclusive exons, and multiple exon skip.

Discovery and translation of cryptic alternative splicing

We detected cryptic alternative splicing in SLE, defined as splicing events with read support in SLE neutrophils and either absence or rare read support in healthy and rheumatoid arthritis disease controls. We used our established methodology (RiboSplitter²⁴) to translate splicing events into peptide sequences to allow downstream analysis and validation. Splicing events with valid unambiguous translations were used for downstream analyses.

Because an alternative splicing event includes two isoforms, we identified cryptic splicing events (both isoforms are cryptic), as well as cryptic splicing isoforms (one of the event’s isoforms is cryptic). We used strict filtering criteria as follows: first, we required cryptic events and isoforms to have ≥10 supporting reads in at least one SLE sample and no more than 2 supporting reads in each healthy or rheumatoid arthritis sample; or ≥5 supporting reads in at least 1 SLE sample and zero supporting reads in each healthy or rheumatoid arthritis sample. Second, we excluded annotated splicing events and isoforms based on the GRCh38 reference. Third, cryptic peptides were excluded if a match existed among known annotated proteins of the gene. These filtering criteria add specificity in detecting events that are cryptic and unique to SLE neutrophils.

We consistently calculated PSI values for isoform 2, so that a higher PSI always indicates higher proportion of isoform 2, which makes it easier to interpret which group has more of which isoform.

MHC affinity prediction

We identified each class I and class II MHC haplotypes for each patient with SLE from RNA-seq data using arcasHLA (v0.5.0), which was shown to be highly accurate (100% accurate for class I MHC and >99.7% accurate for class II MHC)⁶³. MHC affinity binding was predicted for cryptic SLE peptides using NetMHCpan-4.1 (mers length = 8–11) and NetMHCIIpan-4.0²⁶ (mers length = 14–16) that use neural networks trained on binding affinity as well as mass spectrometry data. Strong affinity binding was defined as percent rank of predicted affinity ≤2% (compared to a set of 100,000 random natural peptides) with predicted binding affinity <500.

HLA immunoprecipitation

Peptides that bind to class I MHC were immunoprecipitated using a method described previously⁶⁴. In brief, HLA-A, HLA-B and HLA-C (W6/32) or isotype antibody was conjugated to protein A magnetic beads (Millipore-Sigma) at 1 mg ml⁻¹. Cell lysate from 10 million cells was mixed with 100 µl of anti-HLA beads (or isotype control beads). After washing, bound peptides were eluted with 0.2% trifluoroacetic acid followed by a C-18 column clean-up and further analysed by LC–MS.

LC–MS and analysis

Mass spectrometry was performed at the University of Washington’s Proteomics Resource (UWPR). Peptide separation was performed on a Thermo EASY-nLC system using a 120-min gradient at 300 nl min⁻¹. Mobile phase B was increased from 6% to 45% over 89 min, followed by a ramp to 100% B and column re-equilibration. Samples were loaded at 2.5 µl min⁻¹. Mass spectrometric analysis was carried out on an Orbitrap Exploris 480 operated in positive ion mode. The full mass spectrometry scans were acquired at 120,000 resolution (at m/z 200) across an m/z range of 375–1575, with a normalized AGC target of 200%. Targeted parallel reaction monitoring MS/MS scans were acquired at 30,000 resolution using high-energy collisional dissociation (HCD) fragmentation with a normalized collision energy of 30%. The isolation window was set to 1.6 m/z. Scheduled acquisition windows of 0–120 min were used for targeted precursor ions listed in the inclusion list. Individual peptide was further analysed in Skyline⁶⁵.

Mass spectra visualization and processing

The Skyline chromatogram data were imported into and visualized with Skyline (64-bit) 24.1.0.199. Data were filtered to include only precursors with an idotp ≥0.85 and mass error <10 ppm. Then the match of peak shape and retention time of the spectra from the daughter ions to those from precursors were visually inspected and confirmed.

ELISA for neoantigen-specific IgG

Selected neoantigen peptides synthesized by Shanghai RoyoBiotechnology were coated onto medium-binding EIA plates at 1 ng μl⁻¹ in 0.1 M bicarbonate coating buffer, pH 9.5. Plates were then blocked in 2% BSA in 1× phosphate-buffered saline (PBS). Sera from patients with SLE and healthy donors were added to the blocked wells followed by anti-human IgG–horseradish peroxidase. TMB substrates were added to the wells, and the reaction was stopped by 2N sulfuric acid solution. Absorbance at 450 nm and 630 nm were measured with a BioTek Synergy 2 plate reader.

Gene expression

We quantified gene expression using Rsubread featureCounts⁶⁶ (v2.12.3). The average expression of 111 interferon-induced genes (ISGs)⁶⁷ was calculated per sample, and SLE samples were considered having high interferon signatures when the average expression of ISGs was higher than the 95% percentile of ISG expression in healthy controls.

Comparison to PBMC splicing

To test whether cryptic splicing events that we discovered in SLE neutrophils were also seen in PBMCs, we used two publicly available RNA-seq datasets of PBMCs from patients with SLE. One study (GSE122459)⁶⁸ included 20 patients with SLE and 6 healthy controls, and the other (GSE211700)⁶⁹ included 20 patients with SLE and 10 healthy controls. We ran the same analysis pipeline as described above to align the data and detect alternative splicing events. Then we searched for neutrophil SLE cryptic splicing in the PBMC datasets.

T cell activation assay

PBMCs (2 × 10⁶ per well) from patients with SLE were cultured in 24-well plates in RPMI 1640 medium (Gibco, 11-875-093) supplemented with 5% autologous patient serum that was syringe-filtered (0.22 μm), heat-inactivated at 56 °C for 20 min, and re-filtered (0.22 μm). Cells were stimulated with 5 μg ml⁻¹ peptide pools containing three to six neopeptides and supplemented with IL-2 (10 U ml⁻¹; R&D, BT-002-050) on days 4 and 9. Positive controls consisted of stimulation with phytohaemagglutinin-M (10 μg ml⁻¹; Sigma, NC2034236) on day 1, followed by restimulation with CD3/CD28 Dynabeads (Gibco, 11131D) on day 14. Tetanus toxoid (Fisher Scientific/List Biological Labs, NC0603970) served as an additional positive control in a subset of experiments. On day 14, PBMCs were restimulated with individual peptides from each pool, and ELISPOT plates were coated with capture antibody for overnight incubation. Following 18–24 h coculture, triplicate aliquots of 3 × 10⁵ cells were resuspended in CTL medium supplemented with 1% L-glutamine (Gibco, 2503008), plated into coated ELISPOT wells, and incubated for 4 h at 37 °C. Spots were developed using the CTL-Human IFN-γ/IL-17 Double-Color ELISPOT Kit according to the manufacturer’s instructions and quantified with a CTL S6 Universal M2 reader. The remaining cells were collected for flow cytometry analysis.

Flow cytometry

Following 14 days of stimulation, cells were transferred to a 96-well round bottom plate and.washed with 1× PBS. To prevent non-specific binding, cells were resuspended in 1× PBS with 1% bovine serum albumin (BSA) and incubated at room temperature for 20 min. Cells were then stained with a cocktail of fluorochrome-labelled monoclonal antibodies, including: CD3 conjugated to Alex Fluor 488 (Biolegend, 300320) at 1:50, CD4 BV510 (Biolegend, 317444) at 1:50, CD8 PE-CY7 (Biolegend, 344712) at 1:50 and CD69 APC (Biolegend, 985206) at 1:200 in 1× PBS with 1% BSA for 30 min at 4 °C in the dark. Flow cytometry was run on a CytoFlex Flow Cytometer (Beckman Coulter) and analysis was done with FlowJo (Becton, Dickinson and Company).

TCRB sequencing

Following the 14 day T cell culture in the presence of no peptide or neoantigen peptides IFI44L-1 or IFI44L-2, 10⁵ cryopreserved cells were sent to Adaptive Biotechnology for the ClonoSeq procedure. The resulting ~20,000 TCRB sequences were sorted by ‘template frequency’ and the top 4 TCRBs in each sample are illustrated in Fig. 5j with CDR sequences of the most enriched one.

TCRA and TCRB sequencing with TIRTL

TCR pairs were sequenced from the same T cell cultures as the TCRB sequencing but using TIRTL³⁵.

TCR retroviral transduction

TCR nucleotide sequences were synthesized (Twist) and cloned in the MSGV1 retrovirus backbone by Gibson Assembly (Telesis Bio). The sequence was codon-optimized for expression in human cells. The human TCR constant regions were exchanged for mouse TCR constant regions with substitutions to enhance surface expression and promote preferential pairing^70,71. TCR retroviral supernatants were generated by co-transfection of 293GP cells with the MSGV1 vector containing the relevant TCR sequence and RD114 plasmid. Two days after transfection the retroviral supernatants were completed and either used fresh or frozen at −80 °C. Transductions were performed on RetroNectin-coated plates (Takara) as previously described using a modified T-FINDER³⁵ cell line modified to express both human CD8 and CD4 coreceptors. Murine TCRβ constant region expression was assessed by flow cytometry to determine transduction efficiency.

TCR cell line functional assays

Transduced T-FINDER cells (5 × 10⁴) were cocultured with 1 × 10⁵ patient B- lymphoblastoid cell lines pulsed overnight with 10 µg ml⁻¹ of the indicated IFI44 15-mer peptides or a pool of 10-mer peptides spanning the IFI44 coding sequence in 96-well U-bottom plate. After overnight culture, the percentage of GFP⁺ cells within the murine TCRβ-expressing population was quantified by flow cytometry.

Immunofluorescence microscopy

Neutrophils from healthy donors, and patients with rheumatoid arthritis or SLE were seeded onto glass-bottom 96-well plates and incubated at 37 °C for 30 min to adhere to the plate. Cells were then fixed with 4% paraformaldehyde for 15 min, followed by 3× PBS washes and blocked with 10% normal goat serum overnight at 4 °C. Neutrophils were incubated with perforin 4G10 antibody (1:20 dilution) overnight, washed three times with PBS, and then stained with goat anti-mouse IgG AF647-conjugated secondary antibody (ThermoFisher Scientific, A-21236; 1:1000 dilution) and DAPI (ThermoFisher Scientific, 62248) for 45 min. After three additional PBS washes, cells were resuspended in PBS. Images were captured using a Keyence BZ-X810 fluorescence microscope and analysed with Keyence Analyzer software. Control conditions included cells stained with primary antibody only and secondary antibody only in healthy, rheumatoid arthritis and SLE neutrophils.

ENCODE shRNA-seq analysis of cryptic splicing regulation

ENCODE short hairpin RNA (shRNA)-seq datasets generated in K562 cells were used to assess potential effects of RNA-binding proteins (RBPs) on cryptic splicing⁷². Only paired-end, strand-specific polyadenylated mRNA libraries were included. Reads were aligned and alternative splicing events were quantified as described above. SLE-defined cryptic events were matched to ENCODE-derived splicing events using genomic coordinates. Events with at least one supporting read were retained for analysis. Each shRNA knockdown (KD) sample was paired with its matched non-targeting control samples. Control matching was performed at the sample level, and only explicitly matched control samples were used for each KD replicate. For each cryptic splicing event and RBP, knockdown effects were evaluated using a replicate-aware, read-based framework. Let KD_i denote the cryptic-associated read metric for the ith KD biological replicate and let C_i denote the set of matched control samples for that replicate. A replicate-level effect size was defined as

Each biological replicate contributed one independent estimate of knockdown-induced cryptic splicing change. For each event–RBP pair, replicate-level effects were summarized by the mean effect size and assessed for directional consistency across replicates. Candidate regulatory effects were defined as those supported by at least two KD biological replicates, with consistent directionality across replicates and a mean effect size >2 reads. Event-level results were subsequently aggregated at the RBP level to summarize the number of cryptic splicing events observed following RBP knockdown.

eCLIP-based enrichment analysis of RBPs on cryptic splicing events

ENCODE eCLIP peak data generated in K562 cells were used to evaluate the physical association between RBPs and cryptic splicing events⁷³. For each RBP, biological replicates were merged to obtain a single reproducible peak set. Peaks were filtered using stringent criteria, requiring both high statistical significance and binding enrichment (−log₁₀ P value ≥3 and log₂ fold enrichment ≥3), and only peaks with valid genomic coordinates were retained.

For each alternative splicing event, event-specific binding regions were defined to capture the genomic window where RBPs are most likely to influence splice-site choice. These regions were centred on the relevant splice junctions or exon bodies depending on event type (for example, donor or acceptor sites for alternative splice sites, exon boundaries for exon skipping, intron boundaries for intron retention), using a uniform ±300 bp window.

For each RBP, the observed number of target cryptic event regions overlapping at least one filtered eCLIP peak was calculated using bedtools intersect⁷⁴, counting each target region only once. To assess whether the observed overlap exceeded random expectation, a shuffle-based null model was constructed. Cryptic event regions were randomly shuffled across the genome 1,000 times using bedtools shuffle. For each shuffle iteration, the number of overlaps between shuffled target regions and the RBP’s filtered peaks was computed, generating a null distribution of overlap counts. From this null distribution, the mean and standard deviation of expected overlaps were estimated. An empirical p-value was calculated as the proportion of shuffled overlap counts greater than or equal to the observed overlap count, with a pseudocount added to avoid zero values. A z-score was computed to quantify the deviation of the observed overlap from the null expectation. To account for differences in overall binding abundance across RBPs, the total number of filtered peaks per RBP was recorded, and additional normalization metrics were derived, including peak density and a normalized enrichment score obtained by scaling the z-score by the square root of the peak count.

Empirical P values were adjusted for multiple hypothesis testing using the Benjamini–Hochberg procedure. RBPs were prioritized based on a combination of empirical significance, enrichment strength, and normalized enrichment. For enriched RBPs, event-level binding was further examined by intersecting filtered peaks with individual cryptic event regions. When multiple peaks overlapped the same event, peak-level statistics were summarized using the maximum log₂ fold enrichment and maximum −log₁₀ P value. These summaries were used for downstream visualization, event–RBP matrix construction, and network-based analyses.

Statistical testing

We quantified differentially expressed genes and calculated gene expression normalized counts using DESeq2 (v1.36.0)⁷⁵ with Wald test with correction for multiple comparisons using the Benjamini–Hochberg method. Independent categorical data were tested with the χ² test, and nonparametric data were tested using the Wilcoxon rank sum test. Continuous data were correlated with Pearson’s test. All P values are two-sided and a cut-off of 0.05 was used for statistical significance.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.