Multidimensional assessment of alveolar T cells in critically ill patients

Pneumonia is the leading cause of death from an infectious disease in the United States and the most common infection identified in patients admitted to intensive care units (1, 2). Severe pneumonia is also the dominant cause of acute respiratory distress syndrome (ARDS), a devastating respiratory syndrome for which there are no targeted therapies and for which the mortality rate approaches 40% (3). In mouse models of pneumonia and ARDS, regulatory T (Treg) cells, a subset of CD4⁺ T cells that expresses the Foxp3 transcription factor, resolve lung inflammation (4, 5) and orchestrate tissue repair (6, 7). Effector-memory T (Teff-mem) cells are an additional CD4⁺ T cell subset that traffics to inflamed tissues following infection and can promote damaging inflammation. In mice, distinct transcriptional profiles and epigenetic phenomena, particularly DNA methylation, govern the identity and differential function of Treg and Teff-mem cells during homeostasis and lung infection, resolution, and repair (5, 8). We lack a detailed understanding of the T cell response to severe human pneumonia and ARDS.

In cross-sectional clinical studies, analysis of bronchoalveolar lavage (BAL) fluid obtained from patients with ARDS has identified Treg cells by flow-cytometric surface marker profile (4, 9). However, an understanding of the molecular mechanisms that drive alveolar T cell identity and function in human pneumonia requires high-fidelity multidimensional assessment protocols that can be safely repeated over the course of pneumonia and its recovery. Safety concerns involved with alveolar sampling as well as the low numbers of T cells found in BAL fluid present difficulties that have limited the ability to study the pathogenesis and recovery of pneumonia and ARDS at a molecular level. Here, we describe what we believe is a novel method that facilitates safe retrieval followed by immunologic, transcriptional, and epigenetic profiling of T cell subsets from the alveolar space of critically ill patients. We use these technologies along with a bioinformatic approach in T cells obtained from critically ill patients to validate functional genomic elements (Treg cell–specific super-enhancers [Treg-SEs]) predicted in murine systems.

Our approach utilizes a lavage (wash) of the alveolar space collected via nonbronchoscopic BAL (NBBAL) (Figure 1A and Supplemental Video 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.123287DS1). During this bedside procedure, respiratory therapists advance a 16-French catheter via the patient’s endotracheal tube into the distal airways, where nonbacteriostatic saline is instilled and then withdrawn (10). In critically ill patients with suspected pneumonia, neutrophils are the most common cell type found in the alveolar space, and lymphocytes are consistently present (Figure 1B). In order to credential NBBAL as a safe procedure, we conducted a review of 1,224 consecutive NBBALs performed on intubated, mechanically ventilated patients with suspected pneumonia at our institution. We found that the rate of minor complications was 2.3% and the rate of major complications was 0.3% (Table 1).

Figure 1

Nonbronchoscopic BAL (NBBAL) technique and T cell sorting procedure. (A) Schematic representation of the NBBAL procedure. (B) BAL fluid cell characteristics from a 5-year review of 2,695 consecutive NBBALs performed at our institution. (C) Fluorescence-activated cell sorting strategy used to obtain purified regulatory T (Treg) and effector-memory T (Teff-mem) cells from BAL fluid. The percentage of cells in the indicated gate is shown for a representative sample. (D) Sorted Teff-mem and Treg cell numbers used for downstream next-generation sequencing. (E) Nucleic acid mass purified from sorted Teff-mem and Treg cells. Lines represent mean values. Data in D and E represent 3 replicates for each cell type obtained from 3 individual participants. SSC-A, side scatter area; FSC-A, forward scatter area; FSC-H, forward scatter height; see Supplemental Table 2 for fluorochrome abbreviations.

Table 1

Safety profile of the nonbronchoscopic BAL (NBBAL) procedure

To confirm that our workflow could be integrated into routine clinical care, thus obviating the risks and costs associated with research-specific sampling procedures, we performed all analyses with NBBAL fluid (up to a volume of 15 ml not used for clinical diagnostics) collected from mechanically ventilated patients in the medical intensive care unit (MICU) undergoing NBBAL to investigate suspected pneumonia. We identified Treg cells by their CD127^loCD25⁺ status (4). Non-Treg cells were found to be primarily CD62L^loCD45RA^–, consistent with a Teff-mem phenotype (Figure 1C) (11). BAL fluid cellular viability measured by flow cytometry ranged from 86.8% to 98.1%. Fluorescence-activated cell sorting (FACS) yielded purified yet small populations of Teff-mem and Treg cells for downstream processing (Figure 1D). Total RNA and genomic DNA were extracted from purified cell populations (Figure 1E) and prepared for RNA-sequencing and DNA methylation profiling, respectively.

3′ mRNA-sequencing libraries were generated, amplified, and then sequenced using single-end reads on an Illumina NextSeq 500 instrument. We selected the 3′ mRNA-sequencing procedure for its low input requirement and resilience to degraded RNA. Despite low nucleic acid input, high-quality RNA-sequencing data were obtained with this protocol (Supplemental Figure 1, A and B). Differential gene expression analysis and hierarchical clustering demonstrated clear transcriptional differences between Teff-mem and Treg cells, with canonical lineage-characterizing genes segregating by cell type (Figure 2, A and B). Examples include IL7R (encodes CD127), IRF8, and IRAK4 (Teff-mem cells), and FOXP3 and IL2RA (encodes CD25) (Treg cells). Additionally, genes encoding important immune regulatory molecules, including IL10, CTLA4, PDCD1 (encodes PD-1), and TNFRSF18 (encodes GITR), were observed in Treg cells. We further validated our isolation protocol by running gene set enrichment analysis (GSEA) with 2 available curated human T cell gene sets that discriminate conventional T (non-Treg) and Treg cell transcriptional programs (Figure 2C). Thus, transcriptional profiling confirmed the cellular identity of Teff-mem and Treg cells sorted on the basis of their surface marker profile, albeit with significant between-sample heterogeneity.

Figure 2

Transcriptional profiling of bronchoalveolar lavage (BAL) fluid T cell populations. (A) MA plot comparing the gene expression profile of Teff-mem and Treg cells. Genes of interest are annotated. (B) Hierarchical clustering of 1,320 genes with an absolute log₂(fold-change) ≥ 2. Log₂-transformed counts per million are scaled as Z scores across rows. Genes of interest are annotated. FC, fold-change. (C) Gene set enrichment plots testing for enrichment of the gene sets GSE22045_TREG_VS_TCONV_DN and GSE22045_TREG_VS_TCONV_UP, with genes ordered by log₂ difference in average expression comparing Teff-mem versus Treg cells. Normalized enrichment scores (NES) and false discovery rate (FDR) q values are shown. Data represent average values of 3 replicates for each cell type obtained from 3 individual participants.

In parallel with transcriptional profiling, we performed genome-scale single-nucleotide-resolution DNA methylation profiling using a modified reduced-representation bisulfite sequencing (mRRBS) procedure (Figure 3A) (12) on genomic DNA isolated from sorted BAL T cell subsets (Figure 1, D and E). Following a trimming, alignment, and methylation-calling computational pipeline, we found cytosine-phospho-guanine (CpG) methylation to be uniform across read positions (Supplemental Figure 1C). Raw quantification of CpG methylation provided data on the methylation status of nearly 75% of gene promoters in addition to coverage of other genomic elements (Figure 3B). CpG methylation in the data set displayed the classic dip-and-plateau behavior upstream of transcriptional start sites (promoters) with relative hypermethylation across gene bodies (Figure 3C) (13). Genome-wide, CpG methylation within Teff-mem and Treg cells was moderately-to-strongly correlated (Figure 3D). In order to explore differential methylation between cell types and characterize cell type–specific functional genomic elements, we employed Bayesian hierarchical modeling with the Dispersion Shrinkage for Sequencing (DSS) procedure (14) and used this quantification for downstream analysis (Figure 3E).

Figure 3

Methylation profiling of BAL fluid T cell populations. (A) Schematic outlining the modified reduced-representation bisulfite sequencing procedure. An example Agilent Tapestation 4200 image of DNA fragments following MspI digestion and size selection of fragmented DNA is shown. (B) Coverage of promoters, CpG islands, and CpG island shores. (C) Quantification trend plot of DNA methylation (β scores) across gene bodies defined by transcriptional start site (TSS) and transcriptional end site (TES) with 5 kb of flanking sequence on either side. β Scores continuously range from 0 (unmethylated) to 1 (methylated). (D) Scatter/density plot comparing β scores of effector-memory T (Teff-mem) and regulatory T (Treg) cells with raw quantification (SeqMonk bisulfite methylation pipeline over individual reads). (E) Scatter/density plot comparing β scores of Teff-mem and Treg cells with quantification performed using the Dispersion Shrinkage for Sequencing (DSS) procedure. Data represent merged average of 3 replicates for each cell type obtained from 3 individual participants.

A set of 66 hypomethylated Treg-SE elements has been described in mouse T cell populations (15), and we successfully mapped 62 of these regions onto the human genome (Supplemental Table 1). We found many canonical Treg cell genes near these mapped regions, and their expression favored Treg cells over Teff-mem cells for the majority of expression values (Figure 4A). DSS quantification of CpG methylation demonstrated relative hypomethylation across the mapped Treg-SEs within Treg cells (Figure 4, B and C), consistent with increased transcriptional activity of these loci in Treg cells compared with Teff-mem cells. Collectively, sorted T cell subsets from the alveoli of critically ill patients validated functional genomic elements predicted by a murine system.

Figure 4

Analysis of Treg cell–specific super-enhancer elements. (A) MA plot comparing gene expression for genes within 5 kb of Treg cell-specific super-enhancer (Treg-SE) elements. The red line and shaded area represent the loess fit line with 95% confidence interval. Genes of interest are annotated. FC, fold-change. (B) Cumulative distribution function of β scores (DSS quantification) comparing Teff-mem and Treg cells across Treg-SE. (C) β Scores (DSS quantification) for Teff-mem and Treg cells for CpGs contained within Treg-SE and 20 kb of flanking sequence. Data represent merged average of 3 replicates for each cell type obtained from 3 individual participants. *P < 0.05 by 2-way ANOVA with Sidak’s post hoc test; F(DFn = 1, DFd = 80) = 4.096. NS, not significant.

Our method combines the application of a safe, repeatable, on-demand, nonbronchoscopic sampling procedure with FACS, low-input next-generation sequencing techniques, and computational analysis. This multidimensional assessment procedure permits insights into T cell biology within the context of severe pneumonia, a disease with important public health relevance. Additionally, the technology can validate predictions from mouse models, as we have shown with the Treg-SEs.

The NBBAL approach offers a number of potential advantages over traditional bronchoscopic BAL. First, NBBAL is inexpensive; at our institution, NBBAL is over 5 times cheaper than bronchoscopic BAL. Second, NBBAL can be performed by respiratory therapists without direct physician oversight, facilitating timely alveolar sampling for patients admitted to the intensive care unit or who develop new acute lung pathology at night. This availability facilitates rapid diagnosis and guides the initiation of appropriate antimicrobial therapy (16). Third, the NBBAL procedure is repeatable, permitting clinical and biological insight into disease evolution. Finally, NBBAL is safe. Our consecutive series of over 1,000 NBBAL procedures demonstrated favorable complication rates, which are lower than the reported complication rates of bronchoscopic BALs performed for research purposes in critically ill patients (17).

We noted significant transcriptional heterogeneity within FACS-isolated Teff-mem and Treg cell populations (Figure 2B), which is consistent with the well-described diversity and plasticity of individual T cell subtypes (18). One explanation is the insufficiency of surface marker profiles in selecting homogeneous human T cell populations, particularly with regard to FOXP3 expression (19). However, even with substantial between-sample heterogeneity, the GSEA approach enabled external validation of our sorting strategy at the transcriptional level. These findings highlight the limitations of flow-cytometric analyses based on surface marker profiles when investigating T cell subpopulations during lung infection and recovery and the need to perform multidimensional, unsupervised analysis of their molecular profile (20).

Our study is limited by lack of clinical data on the study participants whose T cells underwent detailed molecular profiling. As the goal of our technical feasibility study was to develop and refine a technical workflow, our study would be underpowered to detect differences based on clinical characteristics. As such, we felt the most practical and ethical approach was to use deidentified residual clinical samples, thus avoiding any risk to patients in terms of additional procedures and confidentiality. Ongoing work examining larger cohorts of critically ill patients with pneumonia will begin to address the clinical correlates of molecular profiles.

Combination of a safe bedside clinical procedure with FACS, transcriptional and epigenetic profiling, and a computational biology toolset provides an opportunity to understand the host response to severe pneumonia on new, detailed levels. In addition to T cells, our sampling procedure allows for the simultaneous identification and sorting of other inflammatory cell populations (e.g., monocyte and macrophage populations and neutrophils) for transcriptional and epigenetic profiling, and the residual fluid can be used for pathogen isolation and sequencing and shotgun metagenomic analysis of the microbiome. The ability to obtain high-resolution, multidimensional data from very low-input samples could prove to be a boon to the study of rare but physiologically important cell populations in the lung and other tissues.

Study participants. For this report, BAL fluid was obtained from mechanically ventilated adult patients in the MICU of Northwestern Memorial Hospital undergoing NBBAL to investigate suspected pneumonia. Typical signs and symptoms that prompt NBBAL include fever, worsening hypoxemia, purulent endotracheal tube secretions, leukocytosis, and new opacities on chest radiograph. In contrast to the consecutive retrospective safety cohort, the 3 participants whose T cells underwent detailed molecular profiling represented a convenience sample based on availability. As these were deidentified clinical samples, they were not selected based on any prespecified clinical criteria.

NBBAL. At our institution, NBBALs are performed by respiratory therapists using a 16-French sampling catheter (Mini-BAL Sampling Catheter, Halyard). Exclusion criteria for NBBAL include positive end-expiratory pressure > 10 cmH₂O, ventilator fraction of inspired oxygen > 0.6, internal diameter of the endotracheal tube < 7.0 mm, and increased risk for procedural bleeding (defined by an international normalized ratio > 2 or a platelet count < 50,000 platelets/μl). The NBBAL procedure can be performed at any hour of the day or night.

During the procedure, the ventilator fraction of inspired oxygen is temporarily increased to 1.0. The NBBAL protocol does not require sedation beyond that administered by the bedside nurse as part of routine ICU care or topical anesthesia. The sampling catheter is inserted into the endotracheal tube and directed into the right or left lung as determined by the treating clinician to preferentially sample the site of presumed infection. The catheter is advanced until gentle resistance is felt, suggesting that the catheter tip has wedged into a distal airway. A volume of at least 90 ml (typically 120 ml) of nonbacteriostatic saline is instilled through the catheter. Suction is then applied to retrieve fluid from the distal airways and alveoli into a trap device. Lavage fluid obtained following the instillation of the first 30 ml of saline is routinely discarded to avoid contamination from the upper airway. The remainder of the aliquot reflects the cellular composition of the small airways and alveolar spaces as measured by clinical automated cell counting (10). Return volumes typically range from 10–30 ml, allowing 5–20 ml for research purposes. Patients are monitored after the procedure and adverse events are recorded in a standard data collection form in the electronic medical record. For our study, a small aliquot (up to 15 ml) of lavage fluid not required for clinical purposes was collected in a 15-ml conical tube (Corning) and immediately refrigerated at 4°C for downstream processing. This excess fluid would normally be discarded by the clinical laboratory.

BAL fluid processing and cell sorting. Lavage fluid was brought to the laboratory on ice and processed to prepare a single-cell suspension for FACS. Fluid was filtered through a 40-μm nylon cell strainer and erythrocytes were lysed using lysing buffer (Pharm Lyse, BD Biosciences). Following staining with a viability dye (Fixable Viability Dye eFluor 506, eBioscience) for 15 minutes at room temperature, cells were incubated with a human Fc-blocking reagent (Miltenyi) for 30 minutes at 4°C. Cells were then incubated with a mixture of fluorochrome-labeled antibodies (Supplemental Table 2) for 30 minutes at 4°C before washing and re-suspending in sorting buffer (21).

Cells were sorted at the Northwestern University Robert H. Lurie Comprehensive Cancer Center Flow Cytometry Core facility using a 100-μm nozzle at 40 psi with a 4-way purity precision mode (Special Order Research Project FACSAria III, BD; see Supplemental Table 2 for instrument configuration). Treg cells were identified using sequential gating as singlets/live/CD3ε⁺/CD4⁺/CD127^lo/CD25⁺ cells. Teff-mem cells were identified as singlets/live/CD3ε⁺/CD4⁺/non-Treg/CD62L^lo/CD45RA^– cells (11).

Cells were sorted directly into 350 μl of cell lysis and nucleic acid extraction buffer (RLT Plus, Qiagen) supplemented with 1% β-mercaptoethanol to inhibit RNases. Samples were stored at –80°C until nucleic acid extraction. Total RNA and genomic DNA from sorted cells were then extracted (AllPrep DNA/RNA Micro Kit, Qiagen) and quality and quantity were assessed (TapeStation 4200, Agilent and Qubit 3.0 instrument).

RNA-sequencing. Libraries were prepared using a 3′ mRNA-sequencing library prep kit (QuantSeq FWD, Lexogen) and a PCR add-on kit (Illumina) on an automated liquid handling platform (Bravo, Agilent). Libraries were multiplexed and sequenced using single-end reads (NextSeq 500, Illumina) with a NextSeq 500/550 V2 High Output reagent kit (1 × 75 cycles).

mRRBS. Figure 3A diagrams the mRRBS procedure (12). Genomic DNA was quantified with a fluorometer (Qubit 3.0) and underwent enzymatic fragmentation with MspI (New England BioLabs). After selecting for 100- to 250-bp fragments with solid-phase reversible immobilization beads (MagBio Genomics), DNA was bisulfite converted with the EZ DNA Methylation-Lightning Kit (Zymo Research) per the manufacturer’s protocol. Bisulfite conversion efficiency averaged 99.2% (SD 0.098%), which was calculated by the global CpG methylation of unmethylated λ-bacteriophage DNA (New England BioLabs N3013S) spiked in at a 1:200 mass ratio to each sample. Random priming and adapter ligation were performed with the Pico Methyl-Seq Library Prep Kit (Zymo Research) using Illumina TruSeq indices. The resulting nondirectional libraries were run on a high-sensitivity screen tape (TapeStation 4200, Agilent) to assess the size distribution and overall quality of the amplified library. Six libraries were pooled in an equimolar ratio and sequenced using single-end reads (NextSeq 500, Illumina) with a 500/550 V2 High Output reagent kit (1 × 75 cycles).

Statistics and bioinformatic analysis. Computation-intensive analysis was performed using Genomics Nodes on Northwestern University’s High-Performance Computing Cluster (Quest, Northwestern University Information Technology and Research Computing).

For RNA-sequencing analysis, reads were demultiplexed using bcl2fastq (version 2.17.1.14). Reads were trimmed using BBDuk (version 36.85) with polyA tail removal. Quality was assessed using FastQC (version 0.11.5). Reads were then aligned to the Homo sapiens reference genome (hg38 assembly) using the STAR aligner (version 2.52). Raw read counts were generated by counting uniquely mapped reads over exons using SeqMonk (version 1.41.0) and filtered to protein-coding genes and genes with at least 0.5 counts per million in at least 2 samples. Downstream expression analysis was conducted with the edgeR R/Bioconductor package (version 3.16.5). GSEA was performed using the GSEA 3.0 GSEAPreranked tool testing for enrichment of the gene sets GSE22045_TREG_VS_TCONV_DN and GSE22045_TREG_VS_TCONV_UP, with genes ordered by log difference in average expression.

mRRBS libraries were demultiplexed as above and trimmed using Trim Galore! (version 0.4.3), removing 10 bases from the 5′ end to remove the random priming sequences. Quality was assessed using FastQC. Sequences were then aligned to the Homo sapiens reference genome (hg38 assembly) with Bismark (version 0.16.3). Bismark was also used for methylation extraction, ignoring 1 base at the 3′ end. Bismark coverage (counts) files for cytosines in CpG context were analyzed with the DSS R/Bioconductor package (version 2.26.0) and quantified using the SeqMonk platform with the bisulphite feature methylation pipeline (raw quantification) or DSS as noted in the text and figure legends. Promoters were defined as 1 kb of sequence flanking transcriptional start sites; CpG island shores were defined as 2 kb of sequence flanking CpG islands.

Mapping of Treg-SE elements (15) to the human genome was performed with the UCSC Batch Coordinate Conversion (liftOver) tool (https://genome.ucsc.edu/cgi-bin/hgLiftOver). Regions were mapped to the hg38 genome using default stringency settings (minimum ratio of bases that must remap = 0.1) and not allowing multiple output regions. The published mm9 reference genome coordinates were mapped to the mm10 reference genome before the liftOver to the hg38 reference genome.

All analysis was performed using R (version 3.4), SeqMonk (version 1.41.0), and GraphPad Prism (version 7.04). Statistical tools for each analysis are explicitly described with the results or detailed in the figure legends. Box-and-whisker plots are presented using Tukey’s method. Unless noted, all tests are 2-sided. P values and false discovery rate (FDR) q values ≤ 0.05 are considered to be significant.

Data availability. Raw counts tables (RNA-sequencing) and Bismark coverage files (mRRBS) are available in the online supplement.

Study approval. All procedures involving human participants were approved by the Northwestern University Institutional Review Board (STU00082848). As the study only involved the use of deidentified, discarded clinical samples, waiver of written informed consent was granted by the Institutional Review Board. For Supplemental Video 1, consent was obtained from the patient’s legally authorized representative prior to filming.