High-throughput analysis of lung immune cells in a combinedmurinemodel of agriculture dust-triggered airway inflammation with rheumatoid arthritis

Rheumatoid arthritis (RA)-associated lung disease is a leading cause of mortality in RA, yet the mechanisms linking lung disease and RA remain unknown. Using an established murine model of RA-associated lung disease combining collagen-induced arthritis (CIA) with organic dust extract (ODE)-induced airway inflammation, differences among lung immune cell populations were analyzed by single cell RNA-sequencing. Additionally, four lung myeloid-derived immune cell populations including macrophages, monocytes/macrophages, monocytes, and neutrophils were isolated by fluorescence cell sorting and gene expression was determined by NanoString analysis. Unsupervised clustering revealed 14 discrete clusters among Sham, CIA, ODE, and CIA+ODE treatment groups: 3 neutrophils (inflammatory, resident/transitional, autoreactive/suppressor), 5 macrophages (airspace, differentiating/ recruited, recruited, resident/interstitial, and proliferative airspace), 2 T-cells (differentiating and effector), and a single cluster each of inflammatory monocytes, dendritic cells, B-cells and natural killer cells. Inflammatory monocytes, autoreactive/suppressor neutrophils, and recruited/differentiating macrophages were predominant with arthritis induction (CIA and CIA+ODE). By specific lung cell isolation, several interferon-related and autoimmune genes were disproportionately expressed among CIA and CIA+ODE (e.g.Oasl1,Oas2, Ifit3, Gbp2, Ifi44, and Zbp1), corresponding to RA and RA-associated lung disease. Monocytic myeloid-derived suppressor cells were reduced, while complement genes (e.g. C1s1 and Cfb) were uniquely increased in CIA+ODEmice across cell populations. Recruited and inflammatory macrophages/monocytes and neutrophils expressing interferon-, autoimmune-, and complement-related genes might contribute towards pro-fibrotic inflammatory PLOS ONE PLOSONE | https://doi.org/10.1371/journal.pone.0240707 February 12, 2021 1 / 27 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Citation: Gaurav R, Mikuls TR, Thiele GM, Nelson AJ, Niu M, Guda C, et al. (2021) High-throughput analysis of lung immune cells in a combined murine model of agriculture dust-triggered airway inflammation with rheumatoid arthritis. PLoS ONE 16(2): e0240707. https://doi.org/10.1371/journal. pone.0240707 Editor:Michal A. Olszewski, University of Michigan Health System, UNITED STATES Received: September 29, 2020 Accepted:December 18, 2020 Published: February 12, 2021 Peer Review History: PLOS recognizes the benefits of transparency in the peer review process; therefore, we enable the publication of all of the content of peer review and author responses alongside final, published articles. The editorial history of this article is available here: https://doi.org/10.1371/journal.pone.0240707 Copyright: 2021 Gaurav et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: The raw and analyzed data is deposited in the Gene Expression Omnibus (GEO) database with access number GSE155436. lung responses following airborne biohazard exposures in setting of autoimmune arthritis and could be predictive and/or targeted to reduce disease burden.

Introduction previously [11]. All procedures on mice were done under isofluorane to minimize distress. After every instillation or injection, animals were monitored consistently until they regained consciousness and mobility. They were monitored every day once by the investigators and once by the vivarium staff for any signs of discomfort or distress.

Organic dust extract
Organic dust extract (ODE) was prepared as previously reported [12] to model airway inflammatory disease. Briefly, an aqueous extract of organic dust from swine confinement feeding facilities (microbial-enriched agriculture setting) was prepared by incubating 1 g dust in 10 ml sterile Hank's Balanced Salt Solution (Mediatech, Manassas, VA) for 1 hour at room temperature followed by centrifugation for 10 minutes at 2,850 x g and repeated twice. The end supernate was filter-sterilized with a 0.22 μm syringe filter to remove any microorganisms and coarse particles. Constituents of the extract have been well characterized and include both endotoxin and peptidoglycans [11,12]. ODE stock was prepared and stored at -20˚C in batches; aliquots were diluted for each experiment to a final concentration (vol/vol) of 12.5% in sterile phosphate buffered saline (PBS; pH = 7.4). Endotoxin concentrations ranged from 150-175 EU/mL as determined using the limulus amebocyte lysate assay (Lonza, Walkersville, MD). This concentration of ODE has been previously shown to produce optimal experimental effects and is well-tolerated in mice [11,12].

Animal co-exposure model
The protocol for the co-exposure model has been previously described [11]. Briefly, mice were age-matched and randomized to 4 treatment groups: Sham (saline injection, saline inhalation), collagen-induced arthritis (CIA; CIA injection, saline inhalation), ODE (saline injection, ODE inhalation), and CIA + ODE (CIA injection, ODE inhalation). CIA was induced with two subcutaneous tail injections (100 μg) of chick type II collagen (2 mg/ml) emulsified in Freund's complete adjuvant (FCA) on day 1 and in Freund's incomplete adjuvant (IFA) on day 21. Sham injections and saline inhalation were conducted with sterile PBS. Following an established protocol, 50 μl of intranasal saline or 12.5% ODE daily for 5 weeks (weekends excluded) was used to induce airway inflammatory disease [11,12]. Mouse treatment groups were ran in parallel with euthanization occurring 5 weeks after initiation of treatments.

Arthritis evaluation
Arthritis inflammatory scores were assessed weekly using the semiquantitative, mouse arthritis scoring system provided by Chondrex (www.chondrex.com) as previously described [11]. Scores range from 0 (no inflammation) to 4 (erythema and severe swelling encompassing ankle, foot, and digits). Arthritis evaluation was assessed on 8 mice per treatment group from 3 independent studies.

Single-cell RNA sequencing
For these studies, 2 mice per treatment group were euthanized with isoflurane in a desiccator. Lungs were exposed from the thoracic cavity and perfused with 10 ml heparin-PBS [11]. Harvested lungs were dissociated with gentleMACS dissociator (Miltenyi Biotech, Auburn, CA) in a digestion solution (collagenase I, 0.2 μg/μl + DNase I, 75 U/ml + heparin, 1.5 U/ml, in Dulbecco's Modified Eagle's Media; DMEM) and incubated for 30 minutes at 37˚C in a shaking incubator. Digestion solution activity was neutralized with PBS containing 4 mM EDTA. Red blood cells were lysed with 1 ml ammonium-chloride-potassium (ACK) lysis buffer (Quality Biological, Gaithersburg, MD) for 1 minute and neutralized with ice-cold DMEM (Gibco). Cells were processed for RNAseq in FACS buffer (2% fetal bovine serum (FBS) + 0.1% NaN 3 in PBS). All reagents purchased from Sigma unless otherwise specified.
Single cell suspensions generated from whole lung were quantified and viability tested using a LUNA-FL TM Dual Fluorescence Cell Counter (Logos Biosystems, Annandale, VA). Single cells were then isolated from cell suspensions (100-2,000 cells/μl) using a 10x Chromium controller per manufacturer's suggested protocol (10x Genomics, Pleasanton, CA). Following cell capture, the gel beads in emulsion (GEM)/sample solution was recovered and placed into strip tubes. Reverse transcription was performed on a thermocycler (C1000 Touch™ Thermal Cycler, Bio-Rad, Hercules, CA) per recommended protocol followed by cDNA amplification. Amplified products were solid phase reversible immobilization (SPRI) bead-purified and evaluated by Fragment Analyzer (Agilent, Santa Clara, CA). Twenty-five percent of the cDNA volume was subjected to fragmentation and double-sided SPRIselect (Beckman Coulter, Indianapolis, IN) was used for PCR purification and clean-up. After adaptor ligation, SPRI clean-up was performed and PCR amplification using sample specific indexes for each sample was completed. PCR products were purified, quantified and library size distribution determined by Fragment Analyzer. Libraries were sequenced per the manufacturer's suggested parameters on a NextSeq500 sequencer to an average depth of 50,000 reads per cell.

Single-cell RNA sequencing data processing
Basecall files (BCL) were generated through 10xGenomics Chromium Single cell 3' Solution followed by RNA Sequencing using Nextseq 500 and Nextseq 550. Cellranger mkfastq was used for demultiplexing and to convert BCL files into FASTQ files. FASTQ files were run through Cellranger count to perform alignment (using STAR aligner), filtering, and unique molecular identifier (UMI) counting. Chromium cellular barcodes were used to generate gene-barcode matrices, perform clustering, and do gene expression analyses. Cellranger aggr was used to normalize and pool the results from different samples, followed by the application of Principal Components Analysis (PCA) to change the dimensionality of the datasets. t-SNE (t-Stochastic Neighbor Embedding) was used to visualize the data in a 2-D space. Graph-based unsupervised clustering was then used to cluster the cells. We used Loupe browser [13], R packages including cellranger R-kit [14], complex heatmap [15], and Geom_violin [16] for more indepth analysis to compare genes expression in each cluster compared to all the other clusters and plot the data. The data sets have been deposited to the Gene Expression Omnibus (GEO) database with access number GSE155436.

RNA isolation
The 4 cell-sorted populations were counted, assessed for viability by trypan blue exclusion (>95%), washed and lysed with RLT buffer containing -mercaptoethanol for RNA isolation as per manufacturer's instructions with Qiagen RNAeasy Micro Kit (Qiagen, Germantown, MD).

NanoString nCounter system
Quality and quantity of total RNA was evaluated using a Fragment Analyzer (Agilent, Santa Clara, CA) and Nanodrop (ThermoFisher), respectively. Total RNA (25-50 ng) was hybridized and processed per the manufacturer's suggested protocol with capture and reporter probes to prepare target-probe complexes using reagents from the Mouse Autoimmune profiling panel containing 771 genes (NanoString, Seattle, WA). Complexes were purified, immobilized and aligned on a cartridge for counting on the nCounter system and processed as per the manufacturer's instructions.
For NanoString analyses, three independent studies of 2-3 pooled mice per group/experiment (N = 8 total mice/group) were analyzed by ANOVA with Tukey's multiple comparisons test. Arthritis inflammatory scores over time was also analyzed by ANOVA. Gene expression data were normalized to 20 housekeeping genes, treatment groups (CIA, ODE and CIA +ODE) were compared to Sham, and data plotted as fold-change. ANOVA with Tukey's multiple comparison test was used on myeloid-derived suppressor cell (MDSC) posthoc analysis. Bar graphs were used to depict means with standard errors of the ratio change in MDSCs normalized to Sham (percentile of MDSC treatment group divided by percentile of MDSC Sham group). Statistical analyses were performed using the GraphPad Prism software, version 8.4.3 (GraphPad, San Diego, CA), and statistical significance accepted at p-values <0.05.

Agriculture (swine) exposure-related ODE-induced airway inflammation coupled with arthritis modeling
Consistent with the previously published study describing this model system [11], the highest arthritis inflammatory scores were demonstrated for CIA+ODE at 5 weeks ( Fig 1A). Throughout the 5 weeks, arthritis inflammatory scores were recorded weekly, and as early as 1 week, there were increases in arthritis inflammatory scores in CIA+ODE and CIA alone as compared with Sham. At 4 weeks, there were significant increases in ODE alone. Previous lung pathology investigations [11] reported increases in alveolar and bronchiolar compartment inflammation and increases in cellular aggregates with ODE and CIA+ODE as compared to Sham, but the number and size of cellular aggregates were reduced in CIA+ODE as compared to ODE alone. Instead, pro-lung fibrosis mediators including lung hyaluronan and fibronectin were increased in CIA+ODE as compared to ODE alone [11]. In these current studies, lung sections from the previous study were stained with H&E from each treatment group showing increased inflammation in CIA+ODE ( Fig 1B). Sections were also stained with markers to identify macrophages (CD68), neutrophils (MPO), T cells (CD3), and B cells (CD45R) (Fig 1C and 1D) distributed in the parenchyma, peribronchiolar and perivascular region. Infiltration of recruited CCR2 + inflammatory monocytes was also demonstrated in the CIA and ODE groups, and these cells were further potentiated in the CIA+ODE group (Fig 1E and 1F). These studies lay the foundation for current studies delineating cellular and gene determinations through high-throughput analysis.

ScRNA-seq identifies 14 unique immune cell subsets
The 10x genomics platform was utilized to cumulatively capture all lung cells. In total, 16,822 cells were analyzed with a mean of 42,901 post-normalization reads per cell and 956 median genes per cell. Unsupervised clustering was performed on 11,577 CD45 + cells and plotted on t-distributed Stochastic Neighbor Embedding (t-SNE). Projection of cells was colored based on unique molecular identifier (UMI) count to identify level of transcripts among the cells. The average UMI count range was roughly between 2,000 to 12,000. Cells that were distributed in the middle showed the highest level of transcripts while cells at the top showed the lowest level (Fig 2A). Unsupervised clustering on the t-SNE projected cells revealed 14 unique immune cell subsets coded by different colors and arbitrary numbers ( Fig 2B). Clusters 3, 4 and 8 were identified as neutrophil subsets based on distribution of Csf3r (granulocyte colony stimulating factor receptor [21]) in t-SNE analysis ( Fig 2C). Macrophages were distributed in the middle and were identified with Cd11c (ITGAX) expression in clusters 1, 2, 5, 11, 14, and partially in cluster 10. Monocytes (inflammatory monocytes) were identified with F13a1 expression in cluster 12. Cluster 9 showed high levels of Ccl5 expression suggesting the presence of NK cells. Similarly, Cd19 expressing cells in cluster 7 identified B lymphocytes, and Trbc2 expression in clusters 6 and 13 identified T lymphocytes, along with expression in the NK cell population (cluster 9). Dendritic cells (DCs) were located in cluster 10 and were characterized by Siglech expression, particularly evident in cluster 10a ( Fig 2C).

CIA and ODE drive unique distributions of immune cells within identified clusters
The 4 treatment groups (Sham, CIA, ODE, and CIA+ODE) exhibited unique distributions of lung immune cells among the identified clusters (Fig 3A and 3B). Among the neutrophil clusters, Sham was exclusively represented by cluster 4, but not cluster 3 or 8. In contrast, the CIA group almost entirely showed neutrophil distribution in cluster 3. The ODE group demonstrated selective distribution of neutrophils in cluster 8 with overlap into cluster 4. In the combination exposure CIA+ODE group, there was broader distribution of neutrophils with predominance in cluster 3, but also evidence for distribution in cluster 4 and partially in cluster 8 (Fig 3A and 3B).
Based on the shifts observed in the cell populations among treatment groups, a manual subclustering was performed to delineate exact number of cells distributed among the sub-clusters ( Fig 3C). Among the macrophage clusters, the ODE group had prevalence in clusters 5, 1b and 1c compared to the CIA group. Likewise, the ODE group lacked clusters 1a and 2a. A subset of cluster 12 (12b) and cluster 10 (10b) were [22] unique to the CIA group, while clusters 10a and 12a were unique to the ODE group (Fig 3). The combination group with CIA+ODE showed a mixed population representing CIA and ODE, while leaning more towards the CIA group ( Fig  3).
Lymphocyte populations were confined to clusters 6, 7, 9, and 13, and were represented in all treatment groups, although modest shifts in cell population distribution were observed. Particularly, NK cells (cluster 9) and B cells (cluster 7) were differentially expressed in ODE and CIA treatment groups with apparent shifts from cluster 9a in ODE to cluster 9b in CIA and shifts from cluster 7a in ODE to cluster 7b in CIA, respectively (Fig 3). Similar to the macrophage clusters, the CIA+ODE group portrayed CIA and ODE group while inclining more towards the CIA group (Fig 3).
Signature genes were selected to highlight respective neutrophil populations on the t-SNE plot ( Fig 4D). Ccl3 was selected to highlight "inflammatory neutrophil" as Ccl3 enhances recruitment and activation of neutrophils in a paracrine fashion [23,24]. Because Il1r2 gene encodes for type 2 interleukin-1 receptor and is constitutively expressed in mouse neutrophils [35], it identified all subsets of neutrophils in the t-SNE clusters ( Fig 4D). Autoreactive neutrophils/gMDSCs were exclusively positive for Mmp8 and Ifit3 in t-SNE. Mmp8 is a neutrophil collagenase [41, 42] and Ifit3 codes for interferon induced protein with tetratricopeptide repeats 3, as both are highly upregulated in gMDSCs and can suppress immune response [31, 32].
Sham, CIA and CIA+ODE groups showed similar distribution of airspace (particularly cluster 1a) and recruited macrophages (especially cluster 2a) while the ODE group had a substantial reduction in these macrophage populations with segregation towards the center of the t-SNE plot in clusters 1b and 1c (Fig 3).

Lymphocytes segregate in four clusters among treatment groups
Four discrete lymphocyte clusters were found in the analysis (Fig 6A). Cluster 6 was identified as "T lymphocytes", with increased expression (log2 fold-change range: 4.87-7.04) of Lef1, Igfbp4, Tcf7, Cd3d, and Cd3e (Fig 6B-6D). This population favored type 2 CD4 + cells based upon the expression of differentiating or expanding population of T lymphocytic genes [87][88][89][90][91]. Cluster 6a was more represented by Sham whereas cluster 6b was represented by CIA and CIA+ODE treatment groups. The ODE group had sparse distribution between cluster 6a and 6b (Fig 3). In contrast to cluster 6, cluster 13 exhibited increased expression of genes indicative of activated T lymphocytes including Icos, Thy1, Cd3g, Ikzf2 and Maf (log2 fold-change range: 4.87-5.87) and thus were termed as "effector T lymphocytes" (Fig 6B-6D) with upregulation of co-stimulatory and adaptive immune pathways. Subtle differences were observed in the distribution of activated T lymphocytes among the treatment groups (Fig 3).
Cluster 7 was remarkable for increased gene expression characteristic of B lymphocytes such as Ebf1, Cd79a, Ms4a1, Cd79b and Ighd (log2 fold-change range: 7.29-7.77) (Fig 6B-6D). Along with genes implicated in B-cell differentiation, memory, signaling and autoimmunity, this cluster showed striking similarities with upregulated pathways in the T lymphocyte population (S3 Table) [92][93][94][95]. B lymphocytes had an overall distribution in CIA+ODE group, but largely represented as cluster 7b in the CIA group. The ODE group had very few B lymphocytes with sparse distribution (Fig 3).

Differential gene expression of ex vivo sorted lung neutrophils across treatment groups represent disease progression
To understand the relevance of myeloid-derived lung cells in RA and RA-associated lung disease, these studies sought to determine whether gene expression of sorted lung myeloidderived cells corresponded to disease-specific findings among treatment groups. Lung neutrophils were isolated by fluorescence activated cell sorting (FACS) based on traditional cell surface markers as percent of CD45 + cells that were Ly6C + Ly6G high (Fig 7A and S1 Fig). Sham represented resident neutrophils in the lungs at baseline [102]. By NanoString analysis, upregulated genes of isolated neutrophils resembled the gene expression demonstrated in scRNAseq data by respective treatment groups. Neutrophils isolated from lungs of CIA and CIA +ODE groups (as compared to Sham) demonstrated increased transcript levels of genes involved in autoimmunity as well as genes associated with gMDSCs/autoreactive neutrophils. These included (CIA and CIA+ODE): Ifit3 (21.7 and 7.49-fold), Ifit1 (17.2 and 5.9-fold), Oas2 (17.0 and 6.0-fold), Zbp1 (14.6 and 8.5-fold), Cxcl9 (5.1 and 8.5-fold), and Oas1a (12.5 and 7.2-fold) (Fig 7B and 7C). Interestingly, the CIA+ODE group also showed gene expression that paralleled that of the ODE group, including increased transcript levels of Src (CIA+ODE: 6.1-fold and ODE: 7-fold), Pf4 (7.7 and 5.3-fold), and complement cascade genes such as C2 (9.6 and 6-fold-change), and Cfb (9.2 and 2.5-fold-change) (Fig 7B and 7C). In contrast, the ODE group demonstrated exclusive upregulation of Ltf (4.7-fold-change), Ccl4 (3.2-fold), Ccl3 (2.8-fold) and Il1rn (2.5-fold) as compared to Sham, consistent with the inflammatory neutrophil cluster (Fig 7B and 7C). There was a single complement cascade gene (C1s1) that was exclusively upregulated in CIA+ODE (4.3-fold) (Fig 7C).

Discussion
In this study, scRNA-seq analysis was applied to whole lung immune cells from a mouse model of RA-associated inflammatory lung disease with key findings confirmed in sorted lung cell populations and NanoString analysis. Building upon the preclinical model of RA-associated inflammatory lung disease [11], we report a number of key findings in this study including: (a) identification of 3 unique neutrophil populations including inflammatory, transient and immunosuppressive/autoreactive granulocytes among experimental groups, (b) heterogeneity among 5 macrophage populations including metabolically active, proliferative, differentiating, recruited, and residential with classical (M1) and alternatively (M2)-activated genes, (c) identification of 2 stages of T-lymphocytes (differentiating and effector), a B-cell population and a NK cell cluster, (d) variability in the distribution of cellular clusters among the treatment groups representing RA and RA-associated lung disease (CIA and CIA+ODE groups, respectively), (e) identification of gMDSC and mMDSC populations based on cell surface markers, and (f) identification of unique genes (interferon-related/autoimmune and complement cascade) that are found in a mouse model of RA-associated inflammatory lung disease.
Occupational exposures from farming, construction, mechanics, medical and military waste have been associated with increased risk of development of RA and/or RA-associated lung disease [109][110][111][112][113]. However, precise mechanism(s) of the development of RA-associated lung disease with occupational and/or environmental inflammatory exposures is not known. Working towards identifying these mechanisms, previous studies demonstrated that repeated exposure to microbial-enriched ODE and particularly ODE+CIA increases citrullination and malondialdehyde-acetaldehyde (MAA)-adduction of lung proteins with a corresponding increase in circulating autoantibody concentrations, periarticular bone damage, and increased deposition of extracellular matrix proteins with a reduction in classical airway inflammatory markers [11,114]. These findings suggested a transition of inflammatory lung disease towards a pro-fibrotic phenotype. Leveraging this co-exposure mouse model with CIA and ODE, several immune cell populations exhibiting unique gene expression signatures that were differentially distributed across treatment groups were demonstrated that suggest potential roles in the pathogenesis of RA, inflammatory lung disease, and inflammatory lung disease specific to RA.
Neutrophils have been classically related to inflammation and host response to pathogens [115]. Knowledge of their role in inflammation and homeostasis continues to evolve as various subsets of neutrophils have been identified and proposed based upon steady state, inflammatory, or anti-inflammatory programming [116]. Many intermediate phenotypes have also been defined, further complicating classifications and the proposed roles in both disease and homeostasis [117]. The 3 different neutrophil populations currently identified aptly signify an inflammatory (in ODE), anti-inflammatory (or autoreactive) (in CIA and CIA+ODE) and homeostatic (transient) (in Sham). In the cell sorting/NanoString studies, neutrophils were identified and limited to Ly6C + Ly6G + cells based upon equipment limitations of simultaneous 4-cell capture ability. Future studies to include the Ly6G + Ly6Ccells that were predominately observed with ODE and CIA+ODE could also be informative. Further studies are needed to delineate precisely how these sub-populations program the lung immune response towards inflammatory and pro-fibrotic disease states.
MDSCs have been implicated in RA [118] as well as ILD [119]. However, their relationship to the RA-associated lung disease is not well-established. MDSCs are transient populations representing myeloid cells at various stages of differentiation that suppress immunity and are subdivided into granulocytic (g) or monocytic (m) origin [120][121][122]. MDSCs are identified by their high expression of Nox2, calprotectin (S100a8/S100a9), Mmp8, Mmp9, Cd33, and multiple interferon-inducible genes such as Ifit3, Ifit1, Oas2, Zbp1, Ifi44, Ifi44l and Oas1a [123,124]. By scRNA-seq analysis, resolution of gMDSC was high with population segregation in cluster 3, which was driven by systemic arthritis induction (i.e. CIA and CIA+ODE groups). This finding was further strengthened by the RNA analysis of the sorted Ly6G + Ly6C + neutrophil population of corresponding treatment groups. Posthoc analysis confirmed that the sorted group contained gMDSCs with Ly6G + Ly6C + CD11b high SSC high gating [107,108], but there was no difference across treatment groups. Unsupervised clustering of immune cells did not segregate mMDSCs. However, gene expression of sorted monocyte-macrophage populations based on CD11b and CD11c expression suggested the presence of mMDSC-like properties based upon the immunosuppressive genes that were elevated in CIA and CIA+ODE groups. Using a classical gating strategy (Ly6G -CD11b + Ly6C high SSC low ) for mMDSCs [108], mMDSCs were identified with FACS and found to be decreased in combination (CIA+ODE) exposure group. These findings are consistent with a recent report that showed that the expansion of MDSCs following tofacitinib treatment is inversely related to the progression of ILD in the SKG mouse model of RA-ILD [125]. These collective findings would suggest a potential protective role for lung MDSCs (particularly mMDSCs) in the development RA-related lung disease, and future studies are warranted to understand their role in disease manifestations to potentially develop novel targets for therapeutic interventions.
Macrophages are one of the most versatile immune cells with immense population heterogeneity and diverse functions [126,127]. Macrophages are increasingly appreciated for their role in fibrosis, wound repair and resolution [47,68,128]. In this current study, the metabolically active airspace macrophages, resident, recruited and differentiating macrophages (clusters 1, 11, 5 and 2, respectively) contribute to inflammation and resolution, while the proliferative airspace macrophages (cluster 14) signify self-renewing properties to maintain a steady population in the lungs. These studies potentially open avenues for hypothesis generation based on various non-traditional genes (interferon-related/autoimmune and complement cascade) expressed in the macrophage subsets that have not been previously investigated in health and disease. The unique distribution of various macrophage clusters among the treatment groups, particularly clusters 1b, 1c, 5 and 11 in ODE group, and clusters 1a, 2a and 12b (inflammatory monocytes) in CIA and CIA+ODE groups signify their importance in disease transition from RA to RA-associated lung disease. Correspondingly, infiltrative CCR2 + inflammatory monocytes were demonstrated with CIA and ODE with highest expression in combined exposures. Future studies could determine whether targeting this cell population results in a reduction of disease manifestations.
Subtle differences among lymphocyte populations (clusters 6, 7, 9 and 13) support earlier work demonstrating that B lymphocytes are skewed towards an autoreactive response following airborne biohazard exposure [114]. B lymphocytes have been recognized as one of the major drivers of autoimmunity [129] and are the target of highly effective RA therapies such as rituximab [130]. Colocalization of MAA with autoreactive B lymphocytes in lung tissues of RA-ILD patients [131] further signifies their potential role in the pathogenesis of RA-ILD. While NK cells are considered a bridge between innate and adaptive immune responses [132], targeting cluster 9b could be of interest. Similarly, autoreactive T lymphocytes (perhaps cluster 6b) and cellular phenotypes supporting fibroproliferation with increase in activated fibroblasts with extracellular matrix deposition could be of particular interest in RA-ILD.
In addition to confirming the upregulation of several interferon-induced genes implicated in autoimmunity, complement cascade genes such as C1ra, C1qa, C1qb, C2, and Cfb representing classical and alternative pathways [133] were identified. C1s1 was highlighted among the complement cascade genes due to its high expression, and was invariably upregulated in the CIA+ODE group in all 4 cell-sorted neutrophil and monocyte-macrophage populations. C1 acts as a sensor for self and non-self-recognition and thus plays a major role in self-tolerance [134]. Complement cascade genes are recognized in RA [135] and ILD [136][137][138], but remain overlooked as therapeutic targets [139]. Holers et. al. [135] suggested that RA-ILD has a complement connection, but to our knowledge, this is among the first reports to experimentally identify C1s1 or complement cascade genes in a RA-inflammatory lung disease model.
There are limitations of this study. In general, animal modeling of arthritis-lung disease interactions are overall limited. The over-expressing TNF-alpha transgenic mice that develop an array of connective tissue diseases has been associated with interstitial lung disease and pulmonary hypertension, particularly in female mice [22]. The arthritic SKG mice develop a cellular and fibrotic interstitial pneumonia [140], but SKG mice do not develop compelling evidence of autoimmunity or arthritis following inhalation injury (i.e tobacco smoke or bleomycin) [140]. Our model system utilizing complex agriculture dust exposure-induced airway inflammation in the setting of CIA would be representative of populations whereby inhalant airborne toxicants are recognized risk factors for RA and RA-lung disease. Cigarette smoke exposure is the most-well-environmental established risk for RA, and occupational inhalant exposures (e.g. work exposures in farming, construction, mining, warehouse environments) have been increasingly associated with risk of disease development, particularly among men [112,141,142]. Moreover, it has previously been reported that ODE exposure induced citrullinated and malondialdehyde-acetaldehyde (MAA) modified proteins in the lung tissue [114], antigens that have been implicated in RA pathogenesis. The advantage of this modeling system is that disease development is dependent on airborne biohazard exposures and can be readily modified to investigate other types of exposures. These current studies were focused on male mice because previous work demonstrated that male mice were profoundly more susceptible to CIA+ODE induced adverse effects as compared to female mice [11]. Moreover, it is also recognized that non-arthritic male mice are susceptible to inhalant endotoxin-induced bone loss whereas female mice were protected and that this protection was dependent upon estrogen [143]. Future studies may be warranted to investigate how single airborne biohazard exposures (e.g. endotoxin, pollutants) in the setting of arthritis induction affect lung-arthritis disease as well as specific hormone factors including testosterone, progesterone, and estrogen. Although lung DCs were noted to be limited in this mouse modeling system, future studies could also investigate the role of lung DCs in lung inflammatory disorders associated with RA.
In conclusion, application of scRNA-seq to an animal model combining systemic arthritis induction and environmental inhalant-induced lung inflammation (i.e. RA-associated lung disease) identified unique populations of lung immune cell clusters differentially ascribed to individual treatment conditions. Neutrophil subpopulations and heterogeneous macrophagemonocyte populations were identified in addition to unique genes (interferon-related and complement cascade) that could be contributing to the pathogenesis of RA-associated lung disease. Additionally, this information might inform potential candidates that could be exploited in future investigations examining targeted interventions and the identification of informative disease biomarkers. Granulocytic (g) MDSCs were identified as live, singlets, CD45 + , non-lymphocytes that were Ly6C + Ly6G + CD11b + SSC high . Whereas monocytic (m) MDSCs were identified as live, singlets, CD45 + , non-lymphocytes that were Ly6G -CD11b + Ly6C + SSC low . (TIF) S1