DNA Methylation Changes in Regional Lung Macrophages Are Associated with Metabolic Differences.

A number of pulmonary diseases occur with upper lobe predominance, including cystic fibrosis and smoking-related chronic obstructive pulmonary disease. In the healthy lung, several physiologic and metabolic factors exhibit disparity when comparing the upper lobe of the lung to lower lobe, including differences in oxygenation, ventilation, lymphatic flow, pH, and blood flow. In this study, we asked whether these regional differences in the lung are associated with DNA methylation changes in lung macrophages that could potentially lead to altered cell responsiveness upon subsequent environmental challenge. All analyses were performed using primary lung macrophages collected via bronchoalveolar lavage from healthy human subjects with normal pulmonary function. Epigenome-wide DNA methylation was examined via Infinium MethylationEPIC (850K) array and validated by targeted next-generation bisulfite sequencing. We observed 95 CpG loci with significant differential methylation in lung macrophages, comparing upper lobe to lower lobe (all false discovery rate < 0.05). Several of these genes, including CLIP4, HSH2D, NR4A1, SNX10, and TYK2, have been implicated as participants in inflammatory/immune-related biological processes. Functionally, we identified phenotypic differences in oxygen use, comparing upper versus lower lung macrophages. Our results support a hypothesis that epigenetic changes, specifically DNA methylation, at a multitude of gene loci in lung macrophages are associated with metabolic differences regionally in lung.

INTRODUCTION

Numerous lung diseases display upper lobe predominance, including cystic fibrosis, sarcoidosis, tuberculosis, and smoking-related chronic obstructive pulmonary disease (1). There are well-described physiologic differences in oxygenation and ventilation between the upper and lower lobes of the lung in the upright position (2). The origin of these differences has been largely attributed to gravity’s effects on the ventilation/perfusion ratio (VA/Q) in the lung, which is defined as the amount of inspired air divided by the pulmonary blood flow in an alveolus (3). Metabolic factors, such as lymphatic flow, pH, and blood flow, differ because of these inequalities of the VA/Q in the lung. Additionally, an oxygen gradient occurs in the lung because of this VA/Q from the apex to the base of the lung, with the upper lobe of the lung being notably more oxygen rich than the lower lobe at end expiration in an upright position (1). Hypoxia is known to induce broad effects on various aspects of cell biology, including alteration of the transcription of hundreds of genes that allow cells to adapt to the changing oxygen availability in the local environment (4). However, how metabolic factors contribute to upper lobe predominant lung disease is unclear.

Macrophages are the sentinel innate immune cells of the lung and possess remarkable immune plasticity, with the ability to sense and adapt to the local milieu (5). Lung macrophages play key roles in bacterial recognition and elimination as well as in polarization of innate and adaptive immunity. Depending on the local environment, macrophages sometimes play a role in anti-inflammatory responses, tissue repair, and homeostasis, whereas at other times, they promote inflammatory and phagocytic processes through the complex production of cytokines and cellular interaction (6, 7).

Epigenetics is the study of heritable changes in gene function caused by mechanisms other than changes in the underlying DNA sequence (8). Epigenetic mechanisms have emerged as modulators of host defenses that can lead to a more-prominent immune response and shape the course of inflammation in the host, both driving the production of specific inflammatory mediators and controlling the magnitude of the host response (9). It is now clear that both the innate and adaptive immune responses use epigenetic mechanisms of gene regulation to maintain long-term phenotypes (10). One important mode of epigenetic regulation is DNA methylation (11, 12).

In this study, we investigate whether DNA methylation changes in lung macrophages are associated with the metabolic factors that contribute to the regional phenotypic differences in lung. To examine this question, we initiated an epigenome-wide DNA methylation profiling study of lung macrophages, comparing upper lobe cells to lower lobe–derived cells.

MATERIALS AND METHODS

This study was approved by the Committee for the Protection of Human Subjects at the Geisel School of Medicine at Dartmouth (approval no. 22781). All subjects were healthy adults and nonsmokers without any known underlying lung disease.

Bronchoalveolar lavage and macrophage isolation

Following local anesthesia with nebulized lidocaine and i.v. conscious sedation, subjects underwent flexible bronchoscopy. Briefly, a flexible fiberoptic bronchoscope was inserted transorally and advanced through the vocal cords. Bronchoalveolar lavage (BAL) fluid was obtained from tertiary airways in the right upper lobes (RUL) and right lower lobes (RLL). BAL was performed sequentially in the RUL and RLL with 20 ml of sterile saline followed by 10 ml of air, and this was repeated for a total of five times per airway. Lung macrophages were isolated, as previously described (13, 14). BAL fluid was filtered through two-layer gauze, centrifuged, and washed twice in 0.9% NaCl. Cells were counted with a T10 automated cell counter (Bio-Rad, Hercules, CA). Cytospins were performed using a Shandon Cytospin 3 centrifuge (Thermo Fisher Scientific, Waltham, MA). Briefly, 75,000 cells resuspended in 200 μl of 0.9% NaCl were loaded into a cytology funnel (Fisher Scientific, Pittsburgh, PA) and centrifuged for 10 min. Cells were allowed to air dry and processed for viewing via Hema 3 Stat pack (Fisher). Imaging was done on an Olympus (Waltham, MA) BX41 microscope with DP2-BSW software (version 2007).

DNA methylation array

Epigenome-wide DNA methylation profiling was performed via the Infinium Methylation EPIC Bead Chips (Illumina, San Diego, CA) for the determination of methylation levels of more than 850,000 CpG sites, as previously described (15, 16). DNA was extracted from lung macrophages via Qiagen (Germantown, MD) DNeasy Blood and Tissue Kit. DNA was quantitated on a Qubit 3.0 Fluorometer (Life Technologies, Carlsbad, CA). Bisulfite conversion of DNA was carried out with the Zymo EZ DNA methylation kit (Zymo Research, Irvine, CA), and EPIC array hybridization and scanning were performed at the University of Southern California Molecular Genomics Core.

DNA methylation array data processing

Raw intensity data files from the MethylationEPIC BeadChips were processed by the minfi R/Bioconductor analysis pipeline (version 1.21) (17) with annotation file version ilm10b3.hg19. All samples showed consistent and high-quality profiles. CpGs annotated as single nucleotide polymorphisms, technical probes, and sex chromosome associated, as well as those failing to meet a detection p value < 0.05 in ≥20% samples, were excluded. This preprocessing procedure left 813,096 CpGs with high-quality methylation data in the final data set. DNA methylation microarray data has been deposited into Gene Expression Omnibus under accession no. GSE132547.

Targeted, next-generation bisulfite sequencing of candidate genes

Targeted next-generation bisulfite sequencing (tNGBS) was performed on nine of the original 13 BAL specimens measured in the DNA EPICmethylation array by EpigenDx (Hopkinton, MA). DNA bisulfite modification was done using EZ-96 DNA methylation kit (Zymo), followed by multiplex PCR with Qiagen HotStar Taq and products purified with QIAquick PCR purification kit. Libraries were prepared using the KAPA Library Preparation Kit for Ion Torrent platforms and Ion Xpress Barcode Adapters (Thermo Fisher Scientific). Library products were purified using Agencourt AMPure XP beads (Beckman Coulter, Indianapolis, IN) and quantified using the Qiagen QIAxcel Advanced System. Barcoded samples were then pooled in an equimolar fashion before template preparation and enrichment were performed on the Ion Chef system (Thermo Fisher Scientific) using Ion 520 and Ion 530 Chef reagents. Enriched, template-positive library products were sequenced on the Ion S5 sequencer using Ion 530 sequencing chips (Thermo Fisher Scientific). FASTQ files from the Ion Torrent S5 server were aligned to the local reference database using open-source Bismark Bisulfite Read Mapper with the Bowtie2 alignment algorithm. Methylation levels were calculated in Bismark by dividing the number of methylated reads by the total number of reads.

Extracellular flux assay

Oxygen consumption rate (OCR) in lung macrophages was measured by Seahorse assay (Agilent, Santa Clara, CA), according to manufacturer’s instructions. Briefly, cells were plated at 50,000 per well and allowed to equilibrate overnight. The following morning, cells were washed twice and placed in XF base medium. The Seahorse instrument sensor cartridge was lowered to 200 μm above the cells, creating the transient microchamber. The sensor contains two fluorophores, one measuring oxygen (O2), the other measuring pH. Changes in oxygen concentration and pH were analyzed via XF Mito Stress Test Report Generator software (version 2.0) and reported as OCR in real time, enabling time-resolved, kinetic data of cell metabolism. Seahorse assays were carried out in DartLab, the Immune Monitoring and Flow Cytometry Shared Resource at the Norris Cotton Cancer Center at Dartmouth-Hitchcock Medical Center.

Statistical analysis

To determine the data structure of methylation profiles with respect to sample characteristics, we performed an unsupervised clustering analysis using the recursively partitioned mixture model (RPMM) method (18), assuming two terminal clusters on 10,341 most variable CpGs (all intersample variance > 0.01). RPMM was implemented in the RPMM R package (version 1.25).

To investigate the putative cell-type proportions in all samples, the RefFreeEWAS algorithm (R package version 2.1) (19) was applied to 10,000 most variable CpGs across all samples. The proportion of putative cell types was calculated iteratively for the number of cell types K from 2 to 10. The optimal number of putative cell types K = 2 was selected, as it minimized the variance of the bootstrapped deviance matrix. Subsequently, the two putative cell types by sample were displayed as stratified boxplots for visualization purposes. Putative cell-type 1 proportion was included as a fixed-effect term in the differential methylation analysis linear model, as described below.

The distribution of intersample variance in methylation β value was examined prior to differential methylation analysis: 27,107 CpGs with intersample variance in β values exceeding 0.005 were selected for further investigation. To account for the upper and lower lobe measurements collected from the same subject, the correlation coefficient (0.935) for the 27,107 CpG β values among the unique subjects (n = 13) was first calculated by the duplicateCorrelation function in the limma R/Bioconductor package (version 3.38.3). Differential methylation analysis was implemented by passing logit-transformed β values (i.e., M values), the consensus correlation coefficient for matched subject pairs, into the lmFit and eBayes functions in limma, with adjustment of subject age, sex, and RefFree cell-type 1 as fixed effects and the subject as a random effect in the model, such that Y = β0 + XRUL βRUL + βage Xage + βsex Xsex + βcell type 1 Xcell type 1 + Random Effect (Subject Pair), where Y is the methylation β value, β0 is the intercept, X is a given covariate, β is the respective model coefficient, XRUL is an indicator variable of sample from the right upper lobe, and βRUL is the linear model coefficient associated with it. Data processing, statistical analysis and data visualization R code can be found at https://github.com/Christensen-Lab-Dartmouth/CF_Epigenetics.

RESULTS

Healthy subjects (n = 13) underwent flexible bronchoscopy to obtain BAL fluid from tertiary airways in both the RUL and RLL, and lung macrophages were isolated as described above. Seven female subjects and six male subjects had a mean age of 27.4 y (SD = 5.3). Genome-scale DNA methylation was measured in lung macrophages from both RUL and RLL, and the most variably methylated 10,341 CpGs from the EPIC methylation array were selected for unsupervised model-based clustering with RPMM (Fig. 1). In nearly all patients, matched RUL and RLL samples cluster immediately adjacent to each other, indicating that intersample variance in DNA methylation is larger than intra-sample variance. For all subjects, both the RUL and RLL sample share RPMM cluster membership. The genomic context for the most variable probes was evenly distributed across promoters/nonpromoters, enhancers/nonenhancers, and CpG islands/non-CpG islands.

FIGURE 1.

Unsupervised clustering analysis, upper lobe versus lower lobe macrophages.

RPMM of the 10,000 CpGs (rows) with greatest sample variance across subjects. Individual subjects 1–13 are shown in columns with sample status bar at top: black (upper) and gray (lower). Blue color represents increased sample methylation.

To determine the putative cell-type composition, we applied the RefFreeEWAS algorithm (R package v.2.1, Linux R version 3.5) to 10,000 most variable CpGs, with 500 bootstraps, 10 iterations, and nine presumptive cell types (i.e., K = 2–10), as previously described (16). We chose K = 2 as this value minimized the column variance of the deviance matrix. The distribution of the two putative cell types was very similar between matched RUL and RLL samples (Supplemental Fig. 1). We also explored phenotypic differences between lung macrophages of the upper and lower lobes. Visual evaluation of lung macrophage morphology via cytospin revealed cells similar in presentation from both lobes: typically, a monocytoid appearance with reniform nuclei and ample cytosol (Fig. 2A, 2B), comprising what appeared to be >85% of the total cell population.

FIGURE 2.

BAL cell cytospins.

BAL cell samples were obtained from tertiary airways from both the RUL (A) and RLL (B) of healthy subjects (n = 13). BAL cells from healthy subjects are visually a mostly homogeneous population of cells, primarily lung macrophages with oval to reniform nuclei and abundant cytosol. Images shown are representative of multiple subjects. BAL cells were stained with H&E. Original magnification ×400.

Next, we investigated the relationship between DNA methylation in upper lobe versus lower lobe lung macrophages. Differential analysis comparing RUL and RLL revealeda total of 95 CpGs at a false discovery rate (FDR) < 0.05 (Fig. 3) and 372 CpGs at an FDR < 0.10 (Supplemental Table I). Of the 95 CpGs with an FDR < 0.05, the top 15 CpGs associated with genes of known function are shown in Table I. CpG location varied across multiple University of California, Santa Cruz reference gene groups, including 5′ untranslated region (UTR), TSS1500, TSS200, and body. The observed change in methylation values ranged from −8.2 to 10.7% for upper lobe cells. Upper lobe cells exhibited increased methylation at CpGs tracking to several genes previously implicated in host–pathogen interaction or in inflammatory/immune-related biological processes, including CAP-Gly domain-containing linker protein family member 4 (CLIP4), HSH2D, NR4A1, SNX10, and TYK2.

FIGURE 3.

Epigenome-wide differential methylation in RUL macrophages.

Comparative analysis of DNA methylation in upper versus lower lobe macrophages in healthy subjects identified 95 differentially methylated CpGs (FDR p < 0.05). CpGs hypermethylated in RUL (green), and hypomethylated CpGs (blue) are plotted as log2-fold increase or decrease in methylation M value (x-axis) versus log10 FDR-adjusted p value (y-axis). Statistically significant CpGs associated with specific genes are labeled, and unlabeled points represent CpGs associated with no known gene at that location.

TABLE I.

Top EPIC CpGs with differential methylation in upper lobe macrophages associated with genes of known function

Name	Chromosome	Coordinate	Strand	UCSC Ref Gene Group	Gene	Δβ value	p Value	FDR p Value
cg22265796	chr19	16254043	+	5′-UTR	HSH2D	0.1070	4.05 × 10−6	0.013
cgl9307883	chr19	10484624	+	Body	TYK2	0.1038	8.89 × 10−6	0.017
cg07676849	chr11	71845665	−	TSS1500	FOLR3	0.0747	1.02 × 10−5	0.018
cgl7323256	chr12	52445897	−	5′-UTR/Body	NR4A1	0.0507	1.33 × 10−5	0.019
cg00889217	chr12	9233083	−	Body	A2M	−0.082	1.98 × 10−5	0.026
cgl3427753	chr6	45406009	−	Body	RUNX2	0.049	4.76 × 10−5	0.032
cg25648980	chr12	51719079	+	TSS1500	BIN2	0.0484	4.77 × 10−5	0.032
cg02290197	chr7	26370519	+	5′-UTR	SNX10	0.0508	3.90 × 10−5	0.032
cgl1259727	chr2	175465297	−	5′-UTR	W1PF1	0.0516	4.33 × 10−5	0.032
cgl8255141	Chr12	6487113	−	TSS1500	SCNN1A	0.0742	4.79 × 10−5	0.032
cgl6711835	chr5	59558467	−	5′-UTR	PDE4D	0.0747	3.80 × 10−5	0.032
cg14975420	Chr10	91096289	+	5′-UTR/Body	IFIT3	0.0479	5.95 × 10−5	0.033
cg06086177	chr16	73090646	+	5′-UTR	ZFHX3	0.0465	7.92 × 10−5	0.038
cgl3126638	Chr12	6487080	+	TSS1500	SCNN1A	0.0545	1.01 × 10−4	0.040
cg26118047	chr2	29320360	+	TSS200	CLIP4	0.0590	1.40 × 10−5	0.046

Δβ, increase in the proportion of methylated alleles in upper lobe cells; TSS, transcription start site; UCSC, University of California, Santa Cruz.

Subsequently, we used tNGBS for validation of the results from Infinium MethylationEPIC array. Assays were designed and run for NR4A1 (cg17323256), A2M (cg00889217), RUNX2 (cg13427753), and CLIP4 (cg26118047). A number of the original Infinium EPIC array CpGs were situated within repeat elements (LINE or SINE), and tNGBS assay design was unsuccessful (i.e., HSH2D, TYK2, SNX10, and SCNN1A). Methylation values from tNGBS were consistent with the EPIC array. Percent methylation difference (Δβ) values for select gene-associated CpGs are presented in Table II. Additionally, an example of tNGBS assay design for CLIP4 is illustrated in Supplemental Fig. 2.

TABLE II.

Comparison of percent methylation changes (Δβ) between EPIC array and tNGBS platforms at gene-specific CpGs

Gene	Location	EPICΔβ Value	tNGBSΔβ Value
NR4A1	cgl7323256	0.051	0.011
CLIP4	cg26118047	0.059	0.073
RUNX2	cgl3427753	0.049	0.044
A2M	cg00889217	−0.082	−0.086
PDE4D	cgl6711835	0.075	0.052
ZFHX3	cg06086177	0.046	0.017

At all targeted gene loci, the two methylation platforms, Illumina EPIC array versus tNGBS, showed significant correlations: NR4A1, CLIP4, RUNX2, A2M, PDE4D, and ZFHX3 (Fig. 4). Correlation and p values ranged from 0.769 and p = 1.93e−04 (ZFHX3) to 0.96 and p = 2.79e−10 (PDE4D).

FIGURE 4.

Correlation of DNA methylation analysis between Illumina EPIC methylation array and tNGBS platforms.

Correlation coefficients and linear regression p values were determined for CpGs associated with the genes NR4A1, CLIP4, RUNX2, A2M, PDE4D, and ZFHX3 as methodology validation, directly comparing the Illumina EPIC DNA methylation array and tNGBS assays.

Lung macrophages were also evaluated for differences in metabolic function by measuring basal respiration levels. OCR was assayed as a measure of basal mitochondrial function in cells comparing upper versus lower. The XF Cell Mito Stress Test performed on the Seahorse extracellular flux analyzer revealed a 21.4% difference (p < 0.001) in basal respiration between RUL versus RLL cells, with lower lobe cells exhibiting a higher OCR (Fig. 5).

FIGURE 5.

Oxygen consumption is higher in lower lobe macrophages.

Primary lung macrophages isolated from BAL fluid were plated and allowed to equilibrate in vitro for 16 h. OCR was measured via a two-fluorophore sensor cartridge via the Seahorse Mito Stress assay. Data are pooled from six to eight wells per experiment, with lung macrophages from three separate donors. Values were normalized to percent RUL OCR. ***p < 0.001 by unpaired Student t test.

DISCUSSION

It has been recognized for more than a half century that differences in physiologic and metabolic factors, such as lymphatic flow, ventilation, pH, and oxygenation levels, exist regionally within the lung (2, 3, 20–22). However, the impact of these physiologic differences on disease pathology or the immune response in upper lobe predominant lung diseases is unclear. We hypothesized that regional physiologic differences may contribute to epigenetic variability in lung macrophages isolated from the upper and lower lobes of the lung, as these cells are long lived and would be prone to epigenetic modifications in response to environmental differences. The goal of this study was to determine if epigenetic changes in lung macrophages, specifically DNA methylation changes, exist regionally in the lung (i.e., upper versus lower lobe) that may ultimately predispose upper lobe lung macrophages to respond differently than lower lobe macrophages.

Limited studies to date have examined lung-related DNA methylation profiling or differences in regional innate immune cells of the lung. Many such DNA methylation studies are lung cancer related (23–27). Additionally, several groups have looked at the effects of cigarette smoking on lung-related DNA methylation (28, 29), and various studies have investigated changes in the DNA methylome in the single cell–type setting (30, 31). One of our recent studies has identified differential methylation in cystic fibrosis subjects at 109 gene-associated CpGs in lung macrophages, using the Illumina Infinium MethylationEPIC (850 K) array (16). In addition, we have demonstrated variant regional lung microRNA and cytokine expression associated with hypoxia in lung macrophages and BAL fluid (upper versus lower lobe) (32).

In this study, we report DNA methylation differences from regional lung macrophages, providing a direct comparison of the epigenetic molecular signature of upper lobe versus lower lobe innate immune cells from healthy individuals. Our cell collection methodology (flexible bronchoscopy) is a rare approach used in research studies and affords us a unique opportunity to analyze innate immune cells isolated directly from the airway and alveolar space.

Our initial strategy was to collect lung macrophages of both the upper lobe and lower lobe in healthy subjects to analyze epigenome-wide DNA methylation patterns. DNA methylation profiling via unsupervised clustering analysis showed that matched upper lobe and lower lobe samples cluster immediately adjacent to each other in nearly all subjects. Additionally, we identified 95 differentially methylated CpGs (all FDR < 0.05), of which 80 were hypermethylated in upper lobe cells and 15 were hypomethylated in upper lobe cells as well as 372 differential CpGs (all FDR < 0.10). Differentially methylated CpGs included those associated with genes such as NR4A1, SNX10, SCNN1A, and A2M. The orphan nuclear receptor NR4A1, which has been shown to be rapidly induced in response to stimuli, such as LPS (33), has been linked to the regulation of glucose metabolism and inflammation as well as implicated in hypoxia-induced apoptosis (34, 35). SNX10, a member of the family of sorting nexins, is a protein involved in intracellular protein sorting, trafficking, and signal transduction (36). Moreover, SNX10 has been reported to regulate endosomal morphology, which might be crucial for macrophage function, including phagocytosis and digestion of pathogens, inflammatory response, and Ag presentation (37). The SCNN1A gene encodes for the α subunit of the epithelial sodium channel (ENaC), which plays a critical role in the ion and fluid regulation of the lung. Altered activity of ENaC, which regulates the airway surface liquid layer, might lead to airway dehydration, mucus stasis, and bacterial overgrowth, as can be seen in cystic fibrosis and chronic bronchitis (38); however, an exact role for ENaC in lung macrophages, specifically, has yet to be fully established.

In addition to epigenetic changes in lung macrophages, we have also identified altered oxygen use in upper versus lower lobe cells. Lower lobe macrophages consistently exhibit a higher OCR than upper lobe macrophages. This increased OCR may represent a compensation mechanism to counter the decreased oxygen present in the lower lobe. Mitochondrial and cellular metabolic processes are central to immune cell responses, hypoxia sensing, and apoptosis in addition to their well-known roles of substrate oxidation and ATP production (39–41). Additionally, glucose metabolism is pivotal during macrophage pathogen phagocytosis (42). Although we do not currently have a direct causal link, one gene we have identified as differentially methylated, NR4A1, has been shown to be a transcriptional regulator of multiple glycolytic enzymes (35) and hence could ultimately be connected to changes in oxygen use levels in lung macrophages.

This study has several major strengths. To the best of our knowledge, this is the first study to demonstrate gene-associated DNA methylation differences in lung macrophages, comparing upper versus lower lobe. We have used the Infinium MethylationEPIC array, spanning 850,000 sites across the human genome, and validated these observations via tNGBS. In addition, we provide evidence of phenotypic differences of upper versus lower lobe cells, namely, a divergence in basal metabolism. However, a limitation of this study is that although we have identified epigenetic and a metabolic phenotypic difference between upper lobe and lower lobe cells in healthy subjects, a causative relationship between the DNA methylation changes and cellular phenotype is yet to be ascertained. Future work with larger sample size and power should investigate such a relationship.

In conclusion, through epigenome-wide DNA methylation profiling, we have identified 95 differential methylated CpGs in lung macrophages from the RUL relative to the RLL. Additionally, we have demonstrated energy metabolome-related phenotypical difference (i.e., oxygen consumption) in regional macrophages, upper versus lower. These observations support a hypothesis that epigenetic changes, specifically DNA methylation, at a multitude of gene loci in lung macrophages are associated with the well-established metabolic regional differences in lung and may, at least in part, be contributing to a differential innate immune response in alternate regions of the lung.