Abstract

Background

   Type I interferon signature is one of the most important features of
   systemic lupus erythematosus (SLE), which indicates an active immune
   response to antigen invasion. Characteristics of type I
   interferon-stimulated genes (ISGs) in SLE patients have not been well
   described thus far.

Methods

   We analyzed 35,842 cells of PBMC single-cell RNA sequencing data of
   five SLE patients and three healthy controls. Thereafter, 178 type I
   ISGs among DEGs of all cell clusters were screened based on the
   Interferome Database and AUCell package was used for ISGs activity
   calculation. To determine whether common ISG features exist in PBMCs
   and kidneys of patients with SLE, we analyzed kidney transcriptomic
   data from patients with lupus nephritis (LN) from the GEO database.
   MRL/lpr mice model were used to verify our findings.

Findings

   We found that monocytes, B cells, dendritic cells, and granulocytes
   were significantly increased in SLE patients, while subsets of T cells
   were significantly decreased. Neutrophils and low-density granulocytes
   (LDGs) exhibited the highest ISG activity. GO and pathway enrichment
   analyses showed that DEGs focused on leukocyte activation, cell
   secretion, and pathogen infection. Thirty-one common ISGs were found
   expressed in both PBMCs and kidneys; these ISGs were also most active
   in neutrophils and LDGs. Transcription factors including PLSCR1, TCF4,
   IRF9 and STAT1 were found to be associated to ISGs expression.
   Consistently, we found granulocyte infiltration in the kidneys of
   MRL/lpr mice. Granulocyte inhibitor Avacopan reduced granulocyte
   infiltration and reversed renal conditions in MRL/lpr mice.

Interpretation

   This study shows for the first time, the use of the AUCell method to
   describe ISG activity of granulocytes in SLE patients. Moreover,
   Avacopan may serve as a granulocyte inhibitor for treatment of lupus
   patients in the future.

   Keywords: Single-cell RNA sequencing, Interferon, Granulocytes,
   Systemic lupus erythematosus, Lupus nephritis
     __________________________________________________________________

Research in context.

Evidence before this study

   The etiology of systemic lupus erythematosus (SLE) is related to
   comprehensive factors, including internal and external causes.
   Ethnicity, gene susceptibility, and environmental predisposing factors
   are highly relevant to the incidence of SLE. Pathogens or autoantigens
   can induce an aberrant immune response in susceptible SLE patients.
   Moreover, type I interferon response is a key contributor to effective
   antiviral responses, which have been found to be highly active in SLE
   patients. Interferon-stimulated genes (ISGs) are a gene set whose
   expression is regulated by interferon. However, the characteristics of
   type I ISGs in SLE and LN have not been clearly described.

Added value of this study

   Our study described granulocyte ISG activity in SLE patients and found
   that transcription factors like PLSCR1, TCF4, IRF9 and STAT1 were
   associated to ISGs expression. Furthermore, Avacopan may serve as a
   granulocyte inhibitor for treatment of lupus patients in future.

Implications of all the available evidence

   Granulocytes display highly active ISG expression in both peripheral
   blood mononuclear cells and in the kidneys of lupus patients. Drugs
   that target granulocytes, such as Avacopan, may reverse
   lupus-associated kidney damage. Regulatory pathways involving
   granulocytes or ISGs may be considered as potential targets for lupus
   treatment.

   Alt-text: Unlabelled box

1. Introduction

   Systemic lupus erythematosus (SLE) is a chronic autoimmune disease that
   can affect many organs, including the kidney, which often indicates
   poor outcome of disease progression [[55]1,[56]2]. Lupus nephritis (LN)
   is a major risk factor for morbidity and mortality in SLE, and 10%–20%
   of patients will progress to end-stage renal disease (ESRD) within 5
   years of diagnosis [[57]3,[58]4]. The etiology of SLE is related to
   comprehensive factors, including internal and external causes.
   Ethnicities, gene susceptibility, and environmental predisposing
   factors are highly relevant to SLE incidence [59][5]. Pathogens are the
   most common inducing factor in the environment, which could induce an
   aberrant immune response in SLE susceptible patients, and type I
   interferon response is a key contributor to effective antiviral
   responses, which have been found to be highly active in SLE patients
   [[60][6], [61][7], [62][8]].

   Type I interferons can be produced by multiple cell types, such as
   monocytes, neutrophils, plasmacytoid dendritic cells, macrophages, and
   epithelial cells, when pathogens invade the human body [63][9].
   Interferon-stimulated genes (ISGs) are a gene set whose expression is
   regulated by interferon. ISGs can be divided based on the activating
   class of interferon, namely type I, type II, or type III
   [[64]10,[65]11]. Specifically, type I interferons are antiviral
   cytokines that trigger ISGs that combat viral infections and are also
   involved in bacterial infections [66][12]. However, the characteristics
   of type I ISGs in SLE and LN have not been clearly described.

   In this study, single-cell RNA sequencing data of peripheral blood
   mononuclear cells (PBMCs) from patients with SLE and microarray
   transcriptomics profile data from the kidneys of LN patients were used
   to investigate the expression characteristics of type I ISGs and the
   potential regulatory mechanism of LN. Furthermore, MRL/lpr and C57BL/6
   mice were used for validation of the results.

2. Methods

2.1. Ethics

   This study was approved by the Ethics Committee of the Chinese People's
   Liberation Army General Hospital (no.S2019-095-01). Informed consent
   was obtained from all participants. For animal experiments, the
   experimental protocol was carried out in accordance with the approved
   guidelines of the Institutional Animal Care and Use Committee at the
   Chinese PLA General Hospital.

2.2. Blood preparation for scRNA-seq

   Fresh peripheral blood was obtained from healthy volunteers and
   patients who were newly diagnosed with SLE, prior to using any
   immunosuppressive drugs. The Human Leucocyte Separation Kit (P8670,
   Solarbio, Beijing) was used for leucocyte separation. After obtaining
   cell pellets, red blood cells were removed using red blood cell lysis
   buffer (R1010, Solarbio, Beijing). A dual fluorescence cell counter
   (LUNA-FL^TM, Logos Biosystems, South Korea) was used for cell count and
   viability test. Single-cell transcriptomic amplification and library
   preparation were performed using single Cell 3′ v3 chemistry (10X
   Genomics) according to the manufacturer's instructions. Approximately 8
   000 cells were captured for each sample following the standard 10X
   capture, and then sequenced using a NovaSeq 6 000 sequencing system
   (PE150, Illumina).

2.3. Single-cell raw data quality control

   The sequencing data from 10 x Genomics were aligned and quantified
   using the CellRanger software package (version 3.1) against the human
   reference genome (hg19). R (version 4.0.2) and the Seurat R package
   (version 3.2.1) were applied for further analyses [67][13]. The
   MergeSeurat function was used to merge multiple datasets. According to
   the median number of genes and the percentage of mitochondrial genes in
   the kidney samples, cells with <200 and >2,500 genes, and a
   mitochondrial gene percentage of >10% were filtered (Supplemental Fig.
   S1A). Gene expression matrices were normalized to the total cellular
   unique molecular identifiers (UMI) count. The normalized expression was
   scaled by regressing the total cellular UMI counts. After data
   normalization, all highly variable genes in single cells were
   identified after controlling for the relationship between average
   expression and dispersion. Subsequently, we used principal component
   analysis (PCA) with variable genes as the input and identified
   significant principal components (PCs) based on the jackStraw function.
   Twenty PCs were selected as the input for uniform manifold
   approximation and projection (UMAP).

2.4. Removing multiplets

   During single cell capture, two or more cells could be captured within
   the same microfluidic droplet and tagged with the same barcode. Such
   artifactual multiplets can confound downstream analyses. We applied the
   DoubletFinder R package to remove these multiplets [68][14]. After
   multiplet removal, 35,842 cells were retained for further analyses
   (Supplemental Fig. S1B).

2.5. Batch effect correction and cell type identification

   Since these data were obtained from different samples, R package
   “Harmony” was used for batch correction (version 1.0) [69][15]
   (Supplemental Fig. S1C) in order to avoid the batch effect disrupting
   downstream analysis. At a resolution of 0.4, cells were clustered by
   the FindClusters function and classified into 20 different cell types.
   Next, we used the FindAllMarkers function to find differentially
   expressed genes of each cell type. Cluster identification was performed
   based on the top differentially expressed genes (DEGs) of each cluster,
   and each gene was manually checked on the CellMarker database to match
   cell types.

2.6. ISGs score

   Differentially expressed genes (logFC > 0.25) of each cluster were set
   as input to generate type I interferon-related ISGs based on the
   Interferome database, and 178 ISGs gene sets were obtained for ISGs
   scoring using the AUCell R package (Version 1.12.0) [70][16]. These 178
   ISG sets were set as input data for area under the curve (AUC) value
   calculation. According to the AUC value, gene-expression rankings were
   built for each cell. The AUC estimates the proportion of genes in the
   gene set that are highly expressed in each cell. Cells expressing many
   genes from the gene set will have higher AUC values than cells
   expressing fewer genes. Function “AUCell_exploreThresholds” was used to
   calculate the threshold that could be used to consider the current
   gene-set active ([71]Fig. 2a). Then, cell clustering UMAP embedding was
   colored based on the AUC score of each cell to show which cell clusters
   were active in the ISG gene set ([72]Fig. 2b).

Fig. 2.

   [73]Fig. 2
   [74]Open in a new tab

   ISG score of SLE PBMC cell clusters. (a) Score of 178 screened ISG
   sets. The threshold was chosen as 0.12 and the ISG score of 12 974
   cells exceeded the threshold value. (b) UMAP plots based on the ISG
   score of each cell. High ISG score cell clusters are highlighted. (c
   and d) GO and pathway enrichment analysis of DEGs in high ISG score
   cell clusters.

2.7. Transcriptomic data obtained from the GEO database

   Raw data from the expression profiles of [75]GSE113342,
   [76]GSE104954/104948, [77]GSE99325/99340, and [78]GSE32592/32591 were
   obtained from the GEO database. Thereafter, differentially expressed
   genes (DEGs) were calculated using the Biobase, GEOquery, and limma R
   package. The DEGs of the dataset with an absolute log2 fold change (FC)
   > 1 and an adjusted P value of <0.05 were considered for subsequent
   analyses.

2.8. Animals

   Female C57BL/6 mice were purchased from Vital River Laboratory
   (Beijing, China) and MRL/lpr mice from the Jackson Laboratory (Bar
   Harbor, ME, USA). Animals were housed in cages in a room with constant
   temperature on a 12 h day/night cycle with access to drinking water and
   food ad libitum.

2.9. Administration of Avacopan to mice

   The Avacopan solution (HY-17627, MedCheExpress, Monmouth Junction, NJ,
   USA) was dissolved in a vehicle solution consisting of
   PEG-400/solutol-HS-15 (42966, Sigma-Aldrich, St. Louis, MO, USA)
   according to the manufacturer's instructions. C57BL/6 mice were divided
   into normal and vehicle control groups. MRL/lpr mice were divided into
   the disease group and Avacopan group. A minimum of 3 mice were
   allocated to each group, and 40 mice in total were used in this study.
   At the beginning of observation, all mice were 12 weeks old. After
   fasting for one night, Avacopan (30 mg/kg) or vehicle
   (PEG-400/solutol-HS-15) were administered to mice every morning by
   gavage for two weeks. Mice were sacrificed at 14, 15, and 16 weeks of
   age, and the kidneys, blood, and urine of mice were collected for
   further study.

2.10. Measurement of mouse urinary protein and creatinine

   Urine albumin was determined using a commercial enzyme-linked
   immunosorbent assay kit (C035-2-1, Jiancheng, Nanjing, China). Urine
   creatinine levels were measured in the same samples using the
   Creatinine Colorimetric Assay Kit (C011-2-1, Jiancheng, Nanjing, China)
   according to the manufacturer's instructions.

2.11. Kidney histology

   Renal tissues fixed in 4% paraformaldehyde were embedded in paraffin,
   and 2 -µm-thick sections were cut for morphological analysis. These
   sections were stained with hematoxylin and eosin (H&E) and periodic
   acid-Schiff (PAS) stains. Images were captured using an Olympus BX53
   microscope. Cell numbers around microvessels (3 fields per mouse (n=3)
   were observed) were counted using ImageJ (NIH, Bethesda, MD, USA).

2.12. Immunofluorescence staining

   Immunofluorescence staining was conducted on paraffin kidney sections
   using standard procedures. Briefly, kidney sections were incubated with
   primary antibodies at 4°C overnight. The primary antibodies used were
   anti-CD16b (Thermo Fisher Scientific Cat# MA1-7633, RRID:AB_2103889).
   After washing with PBS, sections were incubated with fluorescein (FITC)
   or cyanine dye 3 (Cy3)-labeled secondary antibodies at room temperature
   for 1 h. Slides were mounted with 6-diamidino-2-phenylindole (ab104139,
   Abcam). Images were acquired using a confocal microscope (Leica-SP8).
   Positive cell numbers per field (>3 fields were observed) were counted
   using ImageJ (NIH, Bethesda, MD, USA).

2.13. Western blotting and co-immunoprecipitation

   Western blotting was performed on kidney tissue lysates to detect
   IFNα1, IFNAR2, PLSCR1, and TCF4. Proteins were separated using sodium
   dodecyl sulfate/polyacrylamide gel electrophoresis (SDS-PAGE),
   transferred onto a polyvinylidene fluoride (PVDF) membrane, and
   subjected to immunoblot analysis. Blotting was performed using
   antibodies against IFNα1 (LS-C330899, LifeSpan BioSciences, Seattle,
   WA, USA), IFNAR2 (Abcam Cat# ab56070, RRID:AB_880736), PLSCR1
   (ab180518, Abcam), and TCF4 (Abcam Cat# ab185736, RRID:AB_2801300).
   After rinsing with Tris-buffered saline supplemented with 0.1% Tween 20
   (TBST), the membranes were incubated with a horseradish peroxidase
   (HRP)-conjugated anti-rabbit antibody for 60 min at room temperature.
   For co-immunoprecipitation, equal amounts of protein from groups were
   immunoprecipitated with saturating amounts of STAT1 (Abcam Cat# ab2071,
   RRID:AB_2302805) overnight at 4°C, followed by a 1 h incubation with 50
   μL protein A-sepharose beads at 4 °C. The beads were then washed,
   boiled, and the supernatants were immunoblotted with antibodies against
   IRF9 (LS-C666494, LifeSpan BioSciences). Specific signals were detected
   using enhanced chemiluminescence (ECL, Amersham Biosciences, Little
   Chalfont, UK) on the ChemiDoc™ Touch Imaging System (Bio-Rad, Hercules,
   CA, USA). Quantification was performed by measuring the intensity of
   the gels using ImageJ.

2.14. Transmission electron microscopy

   Kidneys were cut into tissue blocks (1 mm^3) and fixed in 2.5%
   glutaraldehyde in 0.01 mol/L phosphate buffer at 4°C, followed by 2%
   osmium tetroxide. They were then dehydrated in a series of graded
   ethanol solutions. Ethanol was then substituted with propylene oxide,
   and the tissue was embedded in epoxy resin. Ultrathin sections were
   double-stained with uranyl acetate and lead and examined under a
   JEM1200EX transmission electron microscope (JOEL) at 80 kV.

2.15. Statistical analyses

   Data are expressed as mean ± standard deviation (SD). The statistical
   significance of comparisons between two groups was analyzed using the
   Student's t test. Comparisons between multiple groups were performed
   using one-way or two-way ANOVA followed by Bonferroni's post hoc test.
   The Shapiro-Wilk test was used to test normality and Levene's test was
   used to assess the equality of variances for all measurement data. A
   minimum of three biological replicates were included for all mice
   experiments, including for histopathology, immunofluorescence, and
   western blotting, all of which were performed once; representative
   experiments are shown. Values of P < 0.05 were considered to indicate
   statistical significance.

2.16. Sample size estimation

   To cover rare cell types in the PBMCs of lupus patients a larger number
   of cells need to be sequenced. An online web tool,
   [79]https://satijalab.org/howmanycells/, was used to calculate how many
   cells we would need to sequence to cover rare cell types in the PBMCs
   of lupus patients. We set the number of cell type as 20, minimum
   fraction of rarest cell type as 0.02 and minimum cells per type as 5,
   and it showed that at least 734 cells in total would have to be sampled
   to obtain a minimum of 5 cells from each of those cell types (99%
   confidence level). For animal experiments, we used software G*Power for
   sample size calculation[80][17]. We used serum creatinine values of
   16-week-old mice from a previous published article to calculate the
   sample size[81][18]. Mean serum creatinine of the MRL/lpr group was
   input as 180 and of the normal group, vehicle group and Avacopan group
   was input as 70. The results showed that 3 mice of each group would
   obtain a power of 0.99.

2.17. Data deposition

   Data generated during this study have been deposited in the Gene
   Expression Omnibus (GEO) with the accession code [82]GSE162577.

3. Results

3.1. Single-cell transcriptomic data revealed the complexity of SLE PBMC

   We analyzed the transcriptomic data of 38 748 PBMC cells from three
   healthy controls (one of which was obtained from our own data, and the
   other two datasets were obtained from the 10 × Genomics official
   website) and five SLE patients (two of which were obtained from our own
   data, and the other three datasets were obtained from the Gene
   Expression Omnibus (GEO) database ([83]GSE142016)) using 10 × Genomics
   sequencing ([84]Fig. 1a). After stringent raw data processing and
   filtration, 35 842 cells were retained for further analyses
   (Supplementary Fig. 3). After normalization of gene expression and
   principal component analyses (PCA), we used graph-based clustering to
   partition the cells into 20 clusters ([85]Fig. 1b). These clusters
   could be assigned to known cell lineages through marker genes or top 3
   differentially expressed genes ([86]Fig. 1c, d). The proportion of each
   cell lineage varied greatly among different individuals ([87]Fig. 1e,
   Supplemental Table S1). Clusters, such as monocytes, B cells, dendritic
   cells, and granulocytes, were significantly increased in SLE patients
   and subsets of T cells were significantly decreased in SLE patients, as
   confirmed by individual quantification of the cell composition
   ([88]Fig. 1f).

Fig. 1.

   [89]Fig. 1
   [90]Open in a new tab

   Single-cell RNA sequencing revealed the complexity of SLE PBMC. (a)
   Pipeline of single cell RNA sequencing data processing. (b) A UMAP plot
   representing the 20 clusters across 258 868 PBMCs from eight
   individuals (five SLE patients and three healthy controls (HC)). (c)
   Violin plots showing expression of marker genes for 20 distinct cell
   types. (d) Heatmap showing expression of the top 3 DEGs in each cell
   type. (e) Bar plots showing the proportion of cell types in each
   sample. (f) Ration comparison of each cell type in SLE and HC groups.
   The red arrow represents significantly elevated cell types in the SLE
   group, while the blue arrow indicates the opposite. Data were analyzed
   using the Chi-Squared Test. * P < 0.01 (For interpretation of the
   references to color in this figure legend, the reader is referred to