Single‐cell RNA‐seq reveals altered NK cell subsets and reduced levels of cytotoxic molecules in patients with ankylosing spondylitis

Abstract Ankylosing spondylitis (AS) is an autoimmune disease with unknown aetiology. To unravel the mechanisms mediating AS pathogenesis, we profiled peripheral blood mononuclear cells (PBMCs) from AS patients and healthy subjects using 10X single‐cell RNA sequencing. The frequencies of immune cell subsets were evaluated by flow cytometry. NK cells were purified from PBMCs using isolation kit and were examined for gene expression by RT‐qPCR. Plasma levels of cytolytic molecules were examined by enzyme‐linked immunosorbent assay. Compared to healthy controls, AS patients showed a significant decrease in total NK cells as well as CD56dim NK subset, whereas CD56bright NK cells were increased. Additionally, impaired expression of cytotoxic genes in NK cells of AS patients was observed by bioinformatics algorithm and verified by RT‐qPCR and flow cytometry. Consistent with changes in transcriptomics, we found decreased plasma levels of granzymes, but not granulysin, in AS patients. Furthermore, Pearson correlation analysis revealed a negative correlation between plasma GZMB levels and disease activity (r = −0.5275, p = 0.0358). No correlation was observed between plasma cytolytic molecules and biochemical indexes (ESR and CRP). Our findings uncover altered NK cell subsets and cytotoxic profiles in peripheral circulation of AS patients at single‐cell resolution.


| INTRODUC TI ON
Ankylosing spondylitis (AS) is a branch of spondyloarthritis (SpA), characterized by long-term rheumatic inflammation of unknown origin. The current diagnostic system tends to classify SpA into two categories: peripheral SpA (mainly affecting the extremities, associated with psoriasis, inflammatory bowel disease or preceding infection) and axial SpA (mainly affecting the spine, such as AS). 1 Chronic inflammatory back pain is a characteristic complaint of most patients with AS, presenting as a dull or vague ache that often worsens in the morning or evening. 2 Over time, osteophyte formation and ligament ossification induced by prolonged inflammatory irritation contribute to irreversible structural damage. 3 In some severe cases, complete fusion of the spine causes kyphosis, limited thoracic mobility and complications such as cardiopulmonary and digestive dysfunction.
For now, there is no cure for AS, but reasonable physical exercise and specific medication, such as JAK inhibitor, can effectively alleviate pain and prevent the condition from worsening. 4,5 Although there is no satisfactory explanation for the disease mechanism, many hypotheses have been proposed from clinical practice or animal experiments. HLA-B27 is present in up to 70-80% of patients and its association with AS is still considered to be important. 6 Unlike other HLA-B alleles, due to the specificity of the molecular structure, heavy chain of HLA-B27 readily forms dimers and oligomers that not only bind to immune-related receptors (e.g. KIR3DL2) on the cell surface but also trigger unfolded protein responses, resulting in inflammatory responses. 7 Imbalance of immune cell populations is also thought to be a hallmark of AS. 8 Previous studies have shown that the proportion of CD4 lineage T cell subsets, such as Th17 and Th22, that secrete inflammatory cytokines is elevated in the peripheral circulation of AS patients. 9,10 In another line of research, some scholars linked AS to the changes in cytotoxic cell profiles. A large-scale genotyping of immune-related loci involving 10,619 AS cases and 15,145 controls identified four CD8 + lymphocyte-associated SNPs, including EOMES, IL7R, RUNX3 and ZMIZ1. 11 By analyzing epigenetic, RNA sequencing and protein expression data, Li et al. 12 showed that AS-associated loci were enriched in immune cell types (e.g., monocytes and T cells) and that proteins encoded by genes downregulated in AS patients were enriched in CD8 + T cells and natural killer (NK) cells. Gracey et al. 13 found a changed cytotoxic cell profile in AS patients whose joint fluid had a significantly increased number of CD8 + T cells with activated phenotype. In addition to AS, studies of other diseases in SpA family (such as psoriatic arthritis and SpA with Crohn's disease) consistently support a pathogenic role for CD8 + T cells, suggesting that cytotoxicity-related factors may play a neglected role in the development of SpA. 14,15 Despite the updated knowledge of AS in recent years, research on the role of NK cells in the pathogenesis and pathophysiology of AS has been less satisfactory.
Understanding gene expression differences between phenotypes is crucial for transcriptomics studies. As a novel method for reliable assessment of transcript abundance, single-cell RNA sequencing (scRNA-seq) focuses on gene expression in individual cells, thus providing higher resolution and better understanding of the status of different cell populations. Here, we isolated peripheral blood mononuclear cells (PBMCs) from AS patients and healthy controls and examined a subset of participants (n = 3 per group) using 10X scRNA-seq technology to explore the heterogeneity in immune cell populations and cytotoxic profiles between the two groups. Pilot data generated from scRNA-seq were followed up with flow cytometry, RT-qPCR and enzyme-linked immunosorbent assay (ELISA) in our larger cohort.

| Study subjects
In this study, two cohorts totaling 29 AS patients and 29 healthy controls were enrolled. Detailed characteristics of the participants were summarized in Table S1. All patients with AS fulfilled the modified New York criteria, and active disease was defined by Bath Ankylosing Spondylitis Disease Activity Index (BASDAI) ≥4. 16 This study was approved by the Ethics Committees of Hangzhou Xiaoshan Hospital of Traditional Chinese Medicine. All participants were informed of the details of this experiment and signed consent forms.

| Isolation of PBMCs from fresh peripheral blood
We isolated PBMCs using Lymphoprep density gradient medium (#07801, Stem Cell Technologies) according to the user guide.

| Data processing pipelines
After sequencing, the original BCL files were converted to FASTQ files using Cell Ranger (version 3.1.0). Upon completion of gene expression quantification, we transferred the output data to Seurat R package (version 3.2.3) for subsequent analysis. 17 To further exclude unwanted cells, we set strict quality control criteria: the number of genes identified in a single cell was between 500 and 4000; the number of UMIs in a single cell <20,000 and the percentage of mitochondrial gene expression in a single cell <10%.
Principal component analysis was performed to reduce the number of gene dimensions, and UMAP (uniform manifold approximation and projection) algorithm was run to visualize cells in a two-dimensional space. Thereafter, marker genes of each cluster were identified using "FindAllMarkers" function embedded in Seurat R package. To become a cluster-specific differentially expressed gene (DEG), there were three conditions to be met. First, |log2 fold-change| ≥0.25.
Second, the gene was expressed in more than 25% of cells in the target cluster. Third, adjusted p-value ≤0.05. Cell-type annotation was automatically inferred by "SingleR" package 18 and then manually checked according to known cell-linage-specific genes.

| Function enrichment analysis of marker genes
Gene Ontology (GO) and KEGG are internationally standardized databases that identify enriched biological functions and pathways by comparison with genome-wide background. We used clusterProfiler R package 19

| Flow cytometry
Peripheral blood mononuclear cells were isolated from fresh blood and resuspended at 10 × Data analysis was performed using FlowJo software (version 10).

| Isolation of NK cells and RNA extraction
Natural killer cells were harvested by using human NK cell isolation kit (#17955, Stem Cell Technologies). To assess the efficiency of NK cell sorting, enriched cells were analyzed by flow cytometry. Total RNA was extracted from NK cells using RNA-quick purification kit (#RN001, Yi Shan Biotechnology), followed by RNA quality assay (Nanodrop, Thermo Scientific) and first-strand cDNA synthesis (#K1622, Thermo Scientific). Amplification of cDNA was performed on ABI-7500 system (Applied Biosystems). Primer sequences for RT-qPCR were shown in Table S2. The results of gene expression were calculated using the 2 −ΔΔC T method.

| Detection of cytolytic molecules in plasma
Plasma levels of cytotoxic granules, including GZMA, GZMB and granulysin, were quantified using ELISA kits (#EK1162, #EK1114 and # EK1280, BOSTER Biological Technology). Each plasma sample was diluted with equal volume of dilution buffer prior to testing, and duplicate analyses were performed to ensure the reliability of the assay.

| Statistical analysis
Statistical analysis was performed using GraphPad Prism (version 7.0). Comparisons between AS patients and healthy controls were made using unpaired t-test, and results were shown as mean ± SEM.
Receiver operating characteristic (ROC) curve was adopted to estimate the predictive ability of a specific gene. Area under the ROC curve (AUC) >0.8 implies that the gene is able to distinguish between patients and healthy controls and may be a valuable diagnostic biomarker for AS. Correlations between plasma levels of cytotoxic particles and BASDAI were calculated using Pearson correlation analysis. All tests with p-value <0.05 were considered significant.

| scRNA-seq identified major immune cell populations in peripheral blood
We harvested PBMCs from three AS patients and three healthy subjects by gradient density centrifugation and profiled them based on 10X Genomics platform ( Figure 1A). A total of 53,823 cells were recovered according to the quality control criteria mentioned in "Methods" section, of which 24,810 cells were from the AS group and 29,013 cells from the control (Table S3) Table 1). Apart from the widely known identity markers mentioned above, we obtained substantial novel clusterspecific genes that will be informative for future single-cell studies (Table S4).

| NK cells were depleted in AS patients
To understand which cell population distinguish AS patients from healthy controls, we first attempted to compare the percentage of each cell type between the two groups. Due to the limited number of samples examined by scRNA-seq, we failed to identify the cell types that were uniquely enriched in AS (Table S5). Previous evidence, however, suggests that patients with AS do have multiple immune perturbations or dysregulation, such as marked expansion of certain T cell subset 22 and monocytes 23 in the peripheral circulation.
In the current study, we mainly focused on NK cells and probed their alterations in AS condition, considering that the role of NK cells in the pathogenesis of AS is poorly understood but well studied in some other autoimmune diseases, such as systemic lupus erythematosus (SLE). 24 For this purpose, we carried out flow cytometry and observed a significant reduction of CD3 − CD56 + NK cells in PBMCs, with the mean percentage of NK cells in patients (4.96%) being much lower than in controls (10.53%) (p < 0.001) (Figure 2).

| CD56 dim NK cell subset was reduced in AS patients
Natural killer cells, also known as large granular lymphocytes, are cytotoxic lymphocytes that have the ability to recognize and kill harmful cells without the involvement of MHC molecules and antibodies and are therefore critical to the innate immune system. In the initial clustering atlas, clusters 6 and 13 were authenticated as NK cells (3535 in total), of which 1182 cells were from AS patients and 2353 cells were from controls ( Figure 1B and Table S5). To analyse the NK cell population at a finer scale, we extracted all cells from clusters 6 and 13 and performed a secondary clustering analysis using the Seurat R package. As a result, two subsets of NK cells (named NK0 and NK1) were generated, exhibiting heterogeneous transcriptome properties with distinctive marker genes ( Figure 3A, B and Table S6).
In terms of the NK0 subset, the preferential expression of cytotoxic genes (such as GZMB and NKG7) as well as FCGR3A (an important receptor for initiating antibody-dependent cellular cytotoxicity) implied a strong killing capacity, much like that previously reported for CD56 dim NK cells 25 (Figure 3C, D). NK1 subset was not dominant in cell numbers; however, this subset was enriched in genes such as GPR183, IL7R, SELL and TCF7, which are important for lymphocyte activation, migration and functional regulation, indicating its identity as CD56 bright NK cells 26 (Figure 3C, E).
Reviewing previous studies, we discovered that patients with certain immune diseases, such as SLE and multiple sclerosis, have

| Impaired expression of cytotoxic genes in NK cells of AS patients compared to healthy controls
To understand whether the reduction in NK cells was accompanied by altered transcripts, we analysed gene expression levels of NK cells between the two groups. Compared to NK cells from healthy controls, 114 genes were upregulated and 59 genes were downregulated in NK cells of AS patients (Table S8). Notably, genes associated with MHC molecules, which are known to be responsible for antigen presentation in immune responses, showed higher levels in AS patients ( Figure 4A). In contrast, genes encoding cytotoxicity-related molecules or receptors were at low levels, such as granzymes (GZMA, GZMB and GZMM) and killer cell lectin-like receptors (KLRB1, KLRC1 and KLRC3); some transcription factors related to immunity and inflammatory regulation (e.g. CEBPB, MAF and JUNB) were also downregulated (Table S8).
Among these DEGs, we mainly focused on changes in cytotoxic profiles and summarized top five cytotoxic genes in Figure 4B (ranked by adjusted p-value). Moreover, we isolated NK cells from by scRNA-seq, with the most pronounced suppression of GZMB expression ( Figure 4C and Figure S2A). ROC curves were plotted to assess the ability of these genes to differentiate between AS patients and healthy controls ( Figure 4D). Consistent with the RT-qPCR results, GZMB displayed optimal predictive power with AUC = 0.9333 (p = 0.0048). To verify the reduction of GZMB expression at the protein level, we performed flow cytometry analysis on AS patients and healthy controls (n = 10 in each group) by intracellular staining with GZMB antibody. Significantly, the percentage of GZMB + NK cells in the total NK cell population was significantly lower in AS patients than in controls (p < 0.01) ( Figure 4E, F).
We further studied the up-and downregulated DEGs by dropping them into GO and KEGG databases, respectively. As shown in Figure S2B, antigen processing and presentation, interferonγ response and T cell activation were significantly upregulated in patients with AS. By contrast, downregulated biological processes such as protein targeting to ER, protein targeting to membrane and translational initiation were identified ( Figure 4G). Additionally, we evaluated affected KEGG pathways in AS patients. Consistent with GO analysis, the upregulated DEGs in NK cells from AS patients were enriched in antigen processing and presentation, T cell receptor signalling pathway and helper T cells differentiation ( Figure S2C). For the downregulated DEGs, we noted that these genes were impactful in a number of pathways with immune or inflammatory context, and in particular, NK cell-mediated cytotoxicity was hampered, which may result from impaired expression of cytotoxicity-related genes encoding granzymes and NK cell receptors ( Figure 4H).

| Plasma levels of granzymes were declined in patients with AS
To learn more about secretory cytolytic products in peripheral blood, we examined the plasma levels of GZMA, GZMB and granulysin in AS patients and healthy controls using ELISA kits. Consistent  (Figure 5D, E). We did not observe either an association between plasma granulysin levels and BASDAI ( Figure 5F) or correlations between the three cytolytic products and biochemical indexes (ESR and CRP) (data not shown).

| DISCUSS ION
In this study, we isolated PBMCs from peripheral blood and pro- The reduced expression of cytotoxicity-related genes and skewing of the NK cell compartment towards CD56 bright phenotype (NK1) prompted us to question whether AS patients also exhibit altered circulating cytotoxicity. In our study, we found decreased plasma levels of granzymes (GZMA and GZMB) in AS patients using ELISA assay (no gross differences in granulysin). In agreement with our findings, Gracey et al. 13  Our study has some limitations. We mainly focused on NK cells and cytotoxic profiles without detailing the transcriptome features of other members (e.g. T cells and monocytes) as well as changes in inflammatory cytokines. We are also aware of the limitations in sample source. Due to the scarcity of human skeletal and ligamentous specimens, we relied on peripheral blood throughout the study.
In conclusion, we demonstrate the differences in NK cell subsets and cytotoxic profiles between AS patients and healthy controls by

ACK N OWLED G M ENTS
We thank Academy of Chinese Medical Sciences, Zhejiang Chinese Medical University for the professional technical support to this work. This study was financially supported by the National Natural Science Foundation of China (No. 81904053).

CO N FLI C T O F I NTE R E S T
The authors confirm that there are no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw data of single-cell RNA sequencing in this study have been deposited with the Sequence Read Archive (SRA) under accession number PRJNA749866.