Global hypermethylation of intestinal epithelial cells is a hallmark feature of neonatal surgical necrotizing enterocolitis

Necrotizing enterocolitis (NEC) remains one of the overall leading causes of death in premature infants, and the pathogenesis is unpredictable and not well characterized. The aim of our study was to determine the molecular phenotype of NEC via transcriptomic and epithelial cell-specific epigenomic analysis, with a specific focus on DNA methylation. Using laser capture microdissection, epithelial cell-specific methylation signatures were characterized by whole-genome bisulfite sequencing of ileal and colonic samples at the time of surgery for NEC and after NEC had healed at reanastomosis (n = 40). RNA sequencing was also performed to determine the transcriptomic profile of these samples, and a comparison was made to the methylome data. We found that surgical NEC has a considerable impact on the epigenome by broadly increasing DNA methylation levels, although these effects are less pronounced in genomic regions associated with the regulation of gene expression. Furthermore, NEC-related DNA methylation signatures were influenced by tissue of origin, with significant differences being noted between colon and ileum. We also identified numerous transcriptional changes in NEC and clear associations between gene expression and DNA methylation. We have defined the intestinal epigenomic and transcriptomic signatures during surgical NEC, which will advance our understanding of disease pathogenesis and may enable the development of novel precision medicine approaches for NEC prediction, diagnosis and phenotyping.


Introduction
Necrotizing enterocolitis remains a leading cause of death in premature infants [1][2][3] for which there remain limited treatment options, no reliable biomarkers and the pathogenesis remains unpredictable. NEC can rapidly take its toll on a growing premature infant to an urgent need of cardiorespiratory support leading to intestinal resection for necrotic bowel a few hours later. The precise disease mechanisms leading to surgical NEC (sNEC) remain ill defined, likely in part because few intestinal samples are available for analysis. In seeking to advance the understanding of NEC pathogenesis, several authors have performed whole genome microarray analysis [4] or bulk RNA sequencing on intestinal resections from infants with NEC [5,6]. While these studies have offered some additional insights into the gene expression profiles of infants with surgical NEC, there is considerable variability in the samples. To combat this limitation, we developed the NEC biorepository [7,8], which provides us with the unique opportunity to advance the field of NEC research with a larger number of samples for sequencing and analysis. These intestinal samples will help to further identify the precise molecular signaling pathways involved in the pathogenesis of sNEC, with the hope to identify infants at the highest risk for the development of this devastating disease.
High-throughput genomic sequencing approaches have enormous potential with regard to unraveling complex heterogeneous phenotypes at the molecular level. One approach is to evaluate the DNA methylation of cytosine-guanine dinucleotides (CpGs) in the promoter, exons, introns, enhancers and CpG island shore regions of genes as they can have a significant impact on the regulation of gene expression. Alterations in the epigenomic molecular phenotype at the level of DNA methylation have been frequently associated with profound transcriptional changes leading to disease pathobiology [9,10]. As we continue toward the era of precision medicine, we sought to apply this approach to infants in the neonatal intensive care unit (NICU) with NEC. Since epithelial cell signaling during NEC has been an important factor in deciphering the pathogenesis of the disease [11,12], we focused this study on the epithelial cell methylome signature. Here we determined the epigenomic molecular phenotype of surgical NEC via genome-wide DNA methylation analysis. Specifically, we identified DNA methylation differences in intestinal epithelial cell genomic DNA within resected ileum or colon that may provide new insights into the biology and mechanisms of NEC pathogenesis.

Results
We prepared a total of n = 40 bisulfite sequencing libraries from individual tissue samples following laser capture microdissection of intestinal epithelial cells in each case. These consisted of n = 12 from non-NEC colon samples, n = 10 from sNEC colon, n = 11 from non-NEC ileum samples and n = 7 from sNEC ileum. Representative images of neonatal gut epithelial cells before and after LCM are shown as an example in Fig. 1a. The resulting DNA samples were subjected to whole-genome bisulfite sequencing (WGBS) in which we sequenced a total of 5,162,557,910 aligned read pairs across all samples. Summary statistics for WGBS sequencing data are presented in Additional file 1: Table S1A.

Hypermethylation of gut epithelium is observed in sNEC
We first sought to identify broad differences in sNEC versus non-NEC samples in both the large intestine (colon) and small intestine (ileum). Multidimensional scaling of these data revealed a clear separation between sNEC and non-NEC samples in colon, but less pronounced global differences between sNEC ileum and non-NEC ileum (Fig. 1b). To further examine the broad patterns of CpG methylation, we generated density plots of CpG methylation levels across the entire genome. These report the global distribution of CpG methylation levels without reference to genomic location. Density plots revealed a global shift toward higher levels of methylation in sNEC colon compared to non-NEC colon (Additional file 2: Figure S1A). This relationship was observed across CpG sites in promoters, exons, introns, intergenic regions, CpG island shores and enhancers, but was less clear in CpG islands (CGIs). We also identified similar, albeit less pronounced, global differences when sNEC ileum was compared to non-NEC ileal samples with the exception that CGI shores and promoters showed little or no difference in the overall methylation level between sNEC and non-NEC ileal samples (Additional file 2: Figure S1B).
On further analysis, we identified CpG sites falling into low methylation (LM) (< 20%), intermediate methylation (IM) (20-80%) and high methylation (HM) (> 80%) categories and found that for both non-NEC ileum and non-NEC colon, approximately 2% of sites were categorized as LM with the remainder being distributed relatively evenly between IM and HM sites. These distributions were similar for both ileal and colonic sNEC samples, although in both cases, we found that proportions of LM and IM sites were reduced and HM sites were increased in NEC (p ≤ 0.0001) (Fig. 1c, Tables 1,2). Further analysis of distinct genomic elements revealed that these distributions are also seen in exons and introns and similarly altered between sNEC and non-NEC samples. In contrast, LM sites account for far greater proportions of CGIs, promoters, CGI shores and, to a lesser extent, enhancers in both colon and ileum. Although the distributions of these were altered in both sample types between sNEC and non-NEC samples for promoters and CGI shores (p ≤ 0.0001), no such differences were observed between sNEC and non-NEC CGIs (Fig. 1c, Tables 1, 2).
We next mapped CpG methylation across each autosome, which allows CpG methylation patterns to be visualized spatially. We plotted methylation levels across all human autosomes for each group of samples (sNEC and non-NEC colon and ileum) (Fig. 2, Additional file 2: Figure S2). These analyses confirmed the relative global hypermethylation of sNEC colon when compared to non-NEC colon. We also confirmed relative hypermethylation of sNEC colon when considering CpG sites in promoters, exons, introns, intergenic regions, enhancers, CGIs and CGI shores (Fig. 2). Similar, though less pronounced, trends were noted when sNEC ileal samples were compared with non-NEC ileal samples (Additional file 2: Figure S2).

Identification of NEC-specific differentially methylated regions: pathway analysis and tissue specificity
Our analysis identified differentially methylated CpG regions (DMRs) that differ between sNEC and non-NEC colon and ileum. Specifically, we identified 5265 autosomal DMRs with a difference in average methylation rate of at least 0.1 (10%) between sNEC and non-NEC colon samples (Fig. 3a). These included 4212 in promoter regions (Fig. 3b, Additional file 1: Table S2A) and 1627 DMRs in gene bodies (introns/exons) (Fig. 3c, Additional file 1: Table S2B). Furthermore, we identified 2785 autosomal DMRs in sNEC versus non-NEC ileum samples with a difference in average methylation rate of at least 0.1 (Fig. 3d). These included 2223 promoter in regions (Fig. 3e, Additional file 1: Table S2C) and 840 in gene bodies (Fig. 3f, Additional file 1: Table S2D). Further analysis revealed that 197 DMRs were shared between colon and ileum (Additional file 1: Table S2E).

Identification of sNEC-specific differentially methylated CpG sites: pathway analysis and tissue specificity
We identified differentially methylated single CpG sites (DMSs) and found that there were 38,809 of these between laser-captured epithelium from sNEC and non-NEC colon with an adjusted p value (q value) of < 0.05 (p value < 0.00043) (Additional file 1: Table S4A). Notably, there were far fewer of these (n = 652) at the same level of significance (q = 0.05) when comparing sNEC and non-NEC ileum (Additional file 1: Table S4B). The 30 DMSs displaying the most highly significant differences (q value) between sNEC and non-NEC colon and sNEC and non-NEC ileum are shown in Tables 3 and 4, respectively.
We validated a subset of DMSs via multiplex PCR and next-generation sequencing of bisulfite converted DNA. As shown in Additional file 2: Figure S5C, we found high correlation between primary WGBS data and the targeted follow-up analyses. Genomic coordinates and methylation levels for the sites assayed are shown in Additional file 1: Table S5.
Given the clear pattern of sNEC-associated hypermethylation in our data, we explored the methylation levels of a number of key DNA methyltransferases and associated factors. Notably, we found that CpG sites located within the upstream region of DNMT3A were broadly hypermethylated in sNEC compared to non-NEC controls. Similar observations were made upstream of TET2 and within an intron of TET3. We also found that CpG sites within intronic regions of DNMT3B and DNMT3L were significantly less methylated in sNEC than in non-NEC control (Additional file 1: Table S7B).

Analysis of NEC-specific transcriptional changes in neonatal intestine
We used RNA sequencing to evaluate the transcriptional signature of human colon ( Fig. 4a-d) and ileum ( Fig. 4eh) during active sNEC compared to non-NEC controls. Using a significance cutoff of p ≤ 0.05, least squares (LS) mean > 1, we identified 1760 mRNAs (37.3%) whose expressions were elevated in sNEC colon versus non-NEC control colon samples and 2596 mRNAs (62.7%) whose expressions were reduced (Fig. 4a, Additional file 1: Table S8A). In contrast, we identified 649 mRNAs (75.7%) whose expressions were elevated in sNEC ileum versus non-NEC control ileum samples and 208 mRNAs (24.3%) whose expressions were reduced (Fig. 4e, Additional file 1: Table S8B). sNEC-specific mRNA expression changes were confirmed for a number of transcripts by quantitative RT-PCR ( Fig. 5a-d).
Functional pathway analysis of differentially expressed transcripts identified in the colon revealed enrichment of genes in pathways for "hepatic fibrosis/hepatic stellate activation" (p = 3.68 × 10 -10 ), "granulocyte adhesion and diapedesis" (p = 8.76 × 10 -9 ) and "axonal guidance" (p = 4.98 × 10 -8 ). These results are summarized in Additional file 2: Figure S11A-C, Figure S12A-C and Additional file 1: Table S9. Functional pathway analysis of differentially expressed transcripts in ileum revealed enrichment of genes in pathways for "acute phase response signaling" (p = 3.84 × 10 -8 ), "granulocyte adhesion and diapedesis" (p = 4.24 × 10 -8 ) and "hepatic fibrosis/hepatic stellate activation" (p = 1.46 × 10 -6 ). These results are summarized in Additional file 2: Figure  S11D-F, Figure S13A-C and Additional file 1: Table S10. It is notable that the enriched pathway, "granulocyte adhesion and diapedesis, " was also found to be significantly enriched in genes containing ileal promoter DMRs (see above). Importantly, in the setting of inflammation, leukocytes can migrate through the blood vessel endothelium to the site of intestinal injury.

DNA methylation changes in sNEC are associated with altered RNA expression patterns
We next searched for overlap between the protein coding RNAs that are differentially expressed (p ≤ 0.05) in sNEC versus non-NEC colon and the corresponding DMRs present in our sNEC versus non-NEC colon data. Similar comparisons were made for the ileal samples. We found 739 such mRNAs encoded by genes containing sNEC DMRs in colon (Additional file 1: Table S11A). Fewer such mRNAs were identified in sNEC ileum (n = 71) (Additional file 1: Table S11B). Functional analysis of colon data using IPA revealed enrichment in a number of pathways including endothelin 1 signaling (p = 2.70 × 10 -4 ), corticotropin releasing hormone signaling (p = 1.37 × 10 -3 ) and hepatic cholestasis (2.62 × 10 -3 ) (Additional file 2: Figure S14, Additional file 1: Table S12). Importantly, HNF4 was again identified as a potential upstream regulator of the genes identified in colon and the HNF4 promoter was found to be hypermethylated in NEC colon, while the expression of its corresponding RNA was reduced (eightfold reduction in sNEC colon relative to non-NEC). Notably, HNF4A was identified as being a likely upstream regulator of genes containing sNECspecific differentially methylated CpGs (see above). Other potential upstream regulators of genes whose expressions and DNA methylation levels are altered in sNEC colon include HNF1A/B and TNF (Additional file 2: Figure S14, Additional file 1: Table S13).
In the majority of colon-specific genes, in which differential expression was correlated with differential methylation, expression was reduced in sNEC colon compared to non-NEC colon (505/739 [68.3%]). With respect to these down-modulated genes, we found that their corresponding promoter regions were more frequently hypermethylated in NEC colon tissues compared to non-NEC controls (428/505 [84.8%]). In contrast, hypermethylation was less frequent in genes with elevated expression in sNEC colon (153/234 [65.4%]). Thus, sNEC colon-specific decreases in gene expression were associated with hypermethylation in corresponding promoters and, concomitantly, promoter hypomethylation was associated with increased gene expression in sNEC.
These findings were confirmed when we compared colon sNEC-specific mRNA expression data to colon sNEC-specific DMS data. We found 7087 differentially methylated CpG sites (q ≤ 0.05) that were within or adjacent to genes encoding mRNAs whose expressions were altered (p ≤ 0.05, LS mean > 1) between sNEC and non-NEC colon (Additional file 1: Table S14). Of these, approximately 63% DMSs (4474/7087) were found in genes whose expressions were reduced in sNEC colon compared to non-NEC controls. Furthermore, 92% (183/198) of DMSs located within promoters of genes whose expressions were decreased in NEC colon were found to display sNEC-specific hypermethylation, whereas only 66% (63/95) of those located in promoters of genes with increased expression in sNEC colon displayed hypermethylation. A similar trend was observed with respect to CGI shores with ~ 92% (428/467) of DMSs in CGIs found to be hypermethylated in genes whose expressions were decreased in sNEC colon compared to ~ 73% (120/164) within genes whose expressions were increased. As with DMRs, these DMS-specific results suggest that hypermethylation of promoter and CGI shores may influence the regulation of sNEC-specific RNA expression.
To further explore these relationships between DNA methylation and gene expression in sNEC, we plotted the distribution of the Pearson's correlation between promoter or gene body methylation and gene expression in colon samples. In Fig. 6, the red lines represent genes for which the methylation difference between sNEC and non-NEC colon samples > 0.2 (20%), and the blue represents genes for which the difference ≤ 0.2. Clearly for genes with large differences in promoter and gene body methylation between sNEC and non-NEC colon, there was a significant negative correlation between gene expression and promoter (Fig. 6a) and gene body (Fig. 6b) methylation. This suggests that sNEC-specific expression changes are likely to be caused (in part) by changes in promoter and gene body methylation levels in some genes.     We next plotted Pearson's correlation between promoter methylation and gene expression in sNEC and non-NEC colon samples (X-axis in Fig. 6c) and Pearson's correlation between promoter methylation and gene expression in sNEC and non-NEC ileum samples (Y-axis in Fig. 6c). This identified a small, but statistically very significant correlation between the two sets of correlations, suggesting that when gene expression and promoter methylation are positively correlated in colon samples, there is a slightly higher chance that the gene expression and promoter methylation of that gene are also positively correlated in ileum samples. Furthermore, when gene expression and promoter methylation were negatively correlated in colon samples, we found it was more likely that the gene expression and promoter methylation of that gene were also negatively    Fig. 6d) and Pearson's correlation between gene body methylation and gene expression in sNEC and non-NEC ileum samples (Y-axis in Fig. 6d). We found a clear correlation between the two sets of analyses, suggesting that when gene expression and gene body methylation are positively correlated in colon samples, there is a higher chance that the gene expression and gene body methylation of that gene will also be positively correlated in ileum samples. Finally, to identify genes whose expressions are most closely associated with DNA methylation, we plotted promoter methylation rate against log2 normalized gene expression. We selected a number of genes meeting the following criteria: (1) their p values of the empirical Bayes t test (using R package limma) of logittransformed promoter methylation between sNEC and control colons ≤ 0.05; (2) the difference in logittransformed promoter methylation between sNEC and control colons ≥ 1; (3), the adjusted p values of the t test for Pearson's correlation between gene expression and promoter methylation ≤ 0.1; and (4) those with adjusted p values ≤ 0.05. Selected genes are shown in Fig. 7. Of the genes we identified, some have already been described as being functionally important in complex intestinal disease phenotypes. For example, the expression of TINAG is highly restricted to the kidney, colon, duodenum and small intestine. Neurexophilin and PC-Esterase Domain Family, Member 4 (NXPE4), which has biased expression in colon, is potentially associated with ulcerative colitis [16], as is P21 activated protein kinase 1 (PAK1) [17,18]. Taken together, these findings reveal that the intestinal epithelial cell methylome and the transcriptome are significantly altered in the pathogenesis of sNEC.

Discussion
We present a comprehensive analysis of the impact of sNEC on the epigenetic signatures of neonatal colonic and ileal epithelial cells. In addition to providing global insight into the structure and organization of the DNA methylome in the context of promoters, CGIs, CGI shores, enhancers, introns, exons and promoters, we present novel findings regarding the biological pathways that may be impacted by sNEC-associated epigenetic changes and the association between sNEC-specific DNA methylation and sNEC-specific RNA expression signatures. Our data are significant because whole genome bisulfite sequencing was performed in an epithelial cell-specific manner after laser capture microdissection of neonatal sNEC and non-NEC colon and ileal epithelium. Thus, the data are cell type-specific and minimally influenced by dilution effects that can be caused by the analysis of genome equivalents from complex mixtures of multiple cell types that may be present in differing proportions during bulk RNA sequencing. Such cell type-specific analyses are particularly important because DNA methylation is influenced significantly by cell lineage [19,20].
Our data show that sNEC-derived epithelial cells are broadly hypermethylated in a genome-wide fashion a b c d relative to their non-NEC counterparts. This trend was evident when data were analyzed via density plots to examine overall levels of CpG methylation, when data were analyzed spatially and when methylation levels were separated by CpG site into low, intermediate and high categories. The hypermethylation of sNEC versus non-NEC colon was clearly evident when considering CpG sites present in each of introns, exons, intergenic regions, promoters and, to some extent, enhancers. Interestingly, however, it was not evident when considering only CpGs located within CGIs and CGI shores, which displayed no obvious difference. Not surprisingly, we identified a multitude of sNECspecific differentially methylated regions (DMRs) and differentially methylated sites (DMSs) in both colon and ileum. Further prospective investigations are needed to determine the timing of onset of the methylation changes seen during NEC and whether these modifications are detectable in the stool prior to NEC onset. These future studies will also examine the association of infant sex and age on DNA methylation signatures. The detection of hypermethylation prior to NEC onset may serve as a potential biomarker for the detection of NEC. Moreover, it is possible that other exposures or medications that premature infants receive could alter the methylation patterns, but additional studies are necessary to answer these questions. Functional analysis of our data using IPA software revealed a number of results that both fit with our current understanding of NEC pathobiology and also demonstrate interesting novel themes which deserve further investigation. These analyses focused on both promoter DMRs and DMSs. The latter approach is designed to capture the broadest possible functional insight. It has the advantage of including single CpG sites that may exist in highly differentially methylated and disease-specific states. If a promoter contains one such site, then it is included in the analysis. In contrast, the use of promoter DMRs between sNEC and non-NEC samples likely captures more subtle effects that may be contributed by multiple CpG sites. In this study, with the goal of focusing on robust biological effects, we only included promoters with > 3 informative CpG sites and an average methylation difference of > 10%.

Pearson's Correlation between Gene Body Methylation and Gene Expression in NEC and non-NEC samples
Given the likely significance of the microbiome on the development of NEC [21], it is notable that colonic sNEC-specific promoter DMRs are enriched for genes within the pathway, "role of pattern recognition receptors in the recognition of bacteria and viruses. " A number of genes in our DMR promoter data encode proteins previously shown to be associated with NEC pathobiology including TLR4 [12], PRKCZ [22], IL18 [23], IL17A [24]. Other genes shown in Additional file 1: Table S3, however, have not yet been associated with NEC, and these are worthy of further investigation. We also discovered that genes containing promoter DMRs were enriched for the "leukocyte extravasation signaling" and "interferon signaling" pathways, which suggests a role for epigenomic regulation of inflammatory processes in the colon affected by NEC. Differences in methylation signatures between sNEC and non-NEC ileum samples were less pronounced but still revealed enrichment for genes within pathways that are aligned with our emerging understanding of the pathobiology of NEC.
Notably, we demonstrated that the most highly significant predicted upstream regulator of the identified genes containing differentially methylated promoters in colon (DMR and DSS analyses) was HNF4A, a nuclear receptor transcription factor that regulates genes involved in intestinal epithelial cell development and function [25]. Accordingly, HNF4A RNA expression levels were also downregulated in our sNEC ileum samples as well as others [22]. HNF4A has also been shown to be important in the pathogenesis of inflammatory bowel disease as variants in human HNF4A are associated with increased risk of disease [26][27][28]. Taken together with our findings, this suggests not only that the expression of the transcriptional regulator HNF4A is significantly altered by the development of sNEC, but also that its gene regulatory activity may be influenced via epigenomic mechanisms.
More focused exploration of functional themes associated with promoter DMRs (Additional file 1: Table S2B) revealed some significant findings. For example, the GPR35 promoter was found to be the most hypermethylated promoter in sNEC colon compared to non-NEC colon. GPR35 is highly expressed in immune and intestinal epithelial cells, and its loss of function has been shown to promote colitis in an experimental animal model [29]. Similarly, the third most hypermethylated promoter in our data when comparing sNEC colon to control is located upstream of ring finger protein 186 (RNF186), which maintains gut homeostasis by controlling ER stress in colonic epithelium [30]. Furthermore, RNF186 gene variants have been shown to be associated with ulcerative colitis in humans. Other examples of promoters upstream of gut-associated genes that may advance our understanding of NEC pathobiology include, CDK5RAP3, which is known to act as a tumor suppressor in gastric cancer [31], PAK1, which is known to contribute to inflammatory conditions of the gastrointestinal tract including colitis [17,32] and NR5A2, which is believed to contribute to the pathogenesis of inflammatory bowel disease [33][34][35].
In a similar fashion, we identified a number of functionally significant genes in which average promoter CpG methylation was significantly hypomethylated in sNEC colon versus non-NEC colon. Examples include, SOX9, an important regulator of cell proliferation in the intestinal epithelium [36,37] and PDE4D, which is known to be associated with ulcerative colitis [38]. Furthermore, the loss of DOK2 has been shown to cause severe colitis in an animal model [39] and is altered at the level of DNA methylation in ulcerative colitis and Crohn's disease in humans [40]. CD244 is a marker of intestinal immune cells, contributing to the protection against enteric pathogens and whose expression is dependent on the presence of the gut microbiota [41]. However, the role that CD224 plays in the context of epithelial cell signaling requires further investigation. Our study has limitations relative to sample procurement and availability of clinical data. Due to the nature of one of our IRB protocols, a waiver of consent precluded us from obtaining clinical data from these infants, including the corrected age of the infant at the time of surgery. Infants with surgical NEC at the University of Pittsburgh typically have their reanastomosis surgery performed when they are near term-corrected age [42]. We acknowledge that these samples may have a component of intestinal adaptation that occurs after a bowel resection. However, it is a rare event that healthy term infants undergo bowel resection or intestinal biopsies and we have tried to capture as many non-NEC pathologies that are available in our biorepository. Thus, it is important to note that control tissue samples are in fact "healed NEC" tissue, obtained during surgical reanastomosis. These samples are the only available source of human tissue for this study design, and the practice of using these samples as non-NEC comparison tissues is common in NEC research. Furthermore, because of the manner in which samples were obtained and the nature of the de-identified tissue acquired, we were not able to obtain complete information regarding neonatal sex for all samples. We have therefore focused only on autosomal data to avoid the complication of X-inactivation to minimize the impact of sex differences.
In summary, we present the first systematic analysis of sNEC-specific DNA methylation changes in human infant gut epithelium. In addition to providing important reference data for future research, we have identified numerous genes and pathways that may be dysregulated at the level of the epigenome in sNEC. The data we present are comprehensive and hypothesis generating and will catalyze further investigation into epigenomic dysregulation and its consequences and potential therapeutic significance in NEC. Furthermore, given the rapid onset of NEC and its silent progression, there is an urgent need to identify biomarkers of NEC that may have utility for identifying at-risk individuals.

Methods
All authors had access to the study data and had reviewed and approved the final manuscript.

Declaration: ethics approval for study population and selection criteria
Intestinal samples were obtained from this study in accordance with the University of Pittsburgh anatomical tissue procurement guidelines and was approved by the University of Pittsburgh Institutional Review Board (IRB) Protocols (PRO09110437 or PRO14070508). Premature infants were recruited under Protocol PRO09110437 at either Children's Hospital of Pittsburgh (CHP) of University of Pittsburgh Medical Center (UPMC) or Magee-Womens Hospital Neonatal Intensive Care Units (NICUs), and consent was obtained by a parent or legal guardian on behalf of their infant. Sample procurement included intestinal tissue if resected for NEC or other non-inflammatory indications such as reanastomosis, spontaneous intestinal perforation or anorectal malformation. In some instances, deidentified intestinal samples were obtained with a waiver of consent and approval of University of Pittsburgh IRB (PRO14070508). In these cases, the only clinical information able to be obtained was location of the resected intestine and the indication for surgery. Intestinal resections were snap frozen and stored at -80 °C until laser capture microscopy.

Laser capture microscopy
Snap-frozen specimens were mounted on appropriate embedding molds (Large, Thermo Scientific # 2219; or Small, Sakura Tissue-Tek # 4566) with clear OCT compound (Optimal Cutting Temperature Embedding Medium) (Fisher HealthCare # 4585) and sectioned on a cryostat instrument (Leica CM 1850 UV), 7 microns). These sections were mounted on membrane slides (Leica PEN-Membrane Slide, 2,0 microns # 11505158), stained with toluidine blue (toluidine blue 0.1% Aqueous, Newcomer Supply #14027) and air-dried (Sanpla Dry Keeper, Sanplatec. Corp). Laser dissection and capture of the chosen tissue region of interest (ROI) were cut by the laser capture microdissection (LCM) Instrument (Leica LMD 7000 microscope operated by the LMD software V7.3). The ROI laser capture tissue samples were subsequently collected into appropriate PCR tubes (Nushbaum, Inc. Large #110614-247 or small 110823-383) for further DNA extraction and subsequent studies.

Whole genome bisulfite sequencing (WGBS)
DNA was extracted from LCM tissue using the Nucleospin Tissue XS Kit (Macherey-Nagel) with the modifications suggested for LCM. Extracted DNA was quantified by using the KAPA hgDNA Quantification and QC Kit (Roche). DNA was sheared with Covaris to a size of ~ 175 bp. WGBS libraries were prepared using the KAPA HyperPrep Kit (Roche). Libraries were bisulfite converted postligation using the EZ DNA Methylation-Direct Kit (Zymo). Post bisulfite conversion, libraries were amplified 12 cycles. Libraries were sequenced on an Illumina HiSeq using 150-bp paired-end reads. DNA sequence reads were quality trimmed, and adaptor sequences were removed using Trim-Galore (https ://www.bioin forma tics.babra ham.ac.uk/proje cts/trim_ galor e/). The reads were aligned to the human reference sequence (GRCh38/hg38) using Bismark in paired-end Bowtie 2 modes. Unaligned paired-end reads were then processed in single-end mode. Read duplicates were removed using Bismark. Methylation was called on paired-end and single-end files and then merged. Differentially methylated sites were confirmed via multiplex PCR and next-generation sequencing of bisulfite converted DNA as previously described [14].

RNA sequencing
Bulk RNA sequencing was performed on ileum and colon samples by the Genome Technology Access Center (GTAC) at the Washington University School of Medicine in St. Louis. Total RNA was isolated using TRIzol (Thermo Fisher). Ribosomal RNA was removed from 1μg of total RNA with RiboErase (KAPA). mRNA was then fragmented and reverse transcribed to yield cDNA using SuperScript III RT enzyme (Life Technologies, per manufacturer's instructions) and random hexamers. A second-strand reaction was performed to yield ds-cDNA and Illumina sequencing adapters ligated to the ends. Ligated fragments were then amplified for 14 cycles using primers incorporating unique index tags. Fragments were sequenced on an Illumina NovaSeq using paired end reads extending 150 bases to a target of 30 M reads. Basecalls and demultiplexing were performed with Illumina's bcl2fastq software and a custom python demultiplexing program with a maximum of one mismatch in the indexing read. mRNA expression was analyzed using Partek Flow software. Adapters were removed and reads were aligned to Genome Reference Consortium Human Build 38 with Bowtie 2 and quantified to Ensembl Transcripts release 96 with an average coverage of 31 and depth of 13. Reads were normalized to total counts per million. Features containing fewer than 10 total normalized reads or a lowest average coverage of 1 across all samples were not included the analysis. Partek flow gene-specific analysis was performed using a multi-model approach based on limma trend, which uses an empirical Bayes method to estimate gene expression.

Quantitative real-time PCR
Tissue was stored in RNA Later (Thermo Fisher), and RNA was extracted with TRIzol (Thermo Fisher). RNA was quantified with a NanoDrop Spectrophotometer (Thermo Fisher). cDNA was made using the QuantiTech Reverse Transcription Kit (QIAGEN) according to the manufacturer's instructions. Quantitative real-time PCR (qRT-PCR) was performed using IQ SYBR Green Supermix and the CFX Connect ™ Real-Time PCR Detection System (Bio-Rad). The expression of genes assessed by qRT-PCR was quantified relative to the housekeeping gene RPLP0.

Data analysis
To identify CpG sites differentially methylated between NEC and control samples, we performed logistical regression with correction for overdispersion, as implemented in the R package methylKit, on the number of methylated cytosines and unmethylated cytosines reported by Bismark. For differentially methylated (CpG) region (DMR) analyses on gene promoters, we defined the promoter region of a gene from 1500-bp upstream to 500 downstream of the transcription start site of that gene. To get the promoter and gene body methylation rate of each gene, we averaged the methylation rates of all CpG sites located inside, respectively, in the promoter region or the gene body. We further applied logit transformation to the methylation rates and used the empirical Bayesian method to test the difference in the promoter/gene body logit-transformed methylation rate between the NEC samples and non-NEC samples.
For analysis of associations between RNA expression and DNA methylation, we used the median of ratios method implemented in the R package DESeq2 to normalize the gene counts in each RNA-seq library and then performed log transformation on the normalized gene counts. We selected genes with the most significant difference (unadjusted p value ≤ 0.05, difference in logit-transformed methylation rate ≥ 1) in the promoter methylation level between the NEC and non-NEC samples. Then for each gene, we calculated the Pearson's correlation between the log2 normalized gene expression across the samples and the corresponding promoter methylation rates. The Student's t test was used to test the significance of the Pearson correlation coefficients. The p values of the tests were then adjusted using the Benjamini and Hochberg's method to control the false discovery rate. A similar analysis was performed on the genes with the most significant difference in the gene body methylation level between the NEC and non-NEC samples.