Chromatin accessibility landscape of stromal subpopulations reveals distinct metabolic and inflammatory features of porcine subcutaneous and visceral adipose tissue

Ethics statement

All methods involving animal experiments were performed strictly in accordance with the recommendations for the Care and Use of Laboratory Animals of the National Institutes of Health of China. The protocol employed in this study was reviewed and approved by the Animal Ethical and Welfare Committee of Gansu Agricultural University (approval number: AEWC-GAU-2019-096).

Animals and sample collection

Two 150-day-old females Bama miniature pigs were chosen from the animal breeding facility of the Chongqing Academy of Animal Science in Chongqing, China. Pigs were placed on high-fat diet until body weight reached 165 ± 5 kg, and then the pigs were humanely sacrificed using electro-stunning followed by severance of blood vessels in the neck. All pigs used in this study were allowed access to feed and water ad libitum and lived under the same normal conditions. SVT (back fat) and VAT (retroperitoneal adipose) were excised from each carcass and used for subsequent experiments.

Isolation of adipose SVF

Three grams of WAT (SAT or VAT) from each pig was minced in a 50 mL tube and digested in 20 mL Hank’s balanced salt solution (Life Technologies) containing 1 mg/mL collagenase type I (Life Technologies) and 1.5% bovine serum albumin (Sigma-Aldrich) at 37 °C in a shaker at 200 rpm for 45 min. The digested homogenates were filtered through 100 µm cell strainer to remove undigested tissues and centrifuged at 800 × g for 15 min at 4 °C. The bottom SVF pellets were washed with 10 mL of PBS and passed through a 40 µm cell strainer to remove larger adipocytes and clumps, and recentrifuged at 800 × g for 15 min at 4 °C.

Low-input DNase-seq

We pooled the SVF (S-SVF or V-SVF) from two biological replicates (each one 5.0 × 10⁴ cells). Approximately 1.0 × 10⁵ SVF cells (S-SVF or V-SVF) were used for liDNase-Seq, performed as previously described (Lu et al., 2016). Briefly, SVF cells were digested in 36 µL lysis buffer (10 mM Tris-HCl, pH 7.5, 10 mM NaCl, 3 mM MgCl₂, 0.1% Triton X-100) and incubated on ice for 5 min. DNase I (3.6 U; Roche) was added to final concentration of 0.1 U/µL and incubated at 37 °C for 5 min. Stop Buffer (80 µL; 10 mM TrisHCl, pH 7.5, 10 mM NaCl, 0.15% SDS, 10 mM EDTA) containing 2 µL Proteinase K (20 mg/mL, Life Technologies) and 20 ng circular carrier DNA (pUC19 DNA, Life technologies) was added and incubated at 50 °C for 1 h to stop the reaction. DNA was purified by phenol–chloroform (Sigma-Aldrich) extraction and precipitated with ethanol containing linear acrylamide (Life Technologies). The sequencing libraries were prepared using NEBNext Ultra DNA Library Prep Kit for Illumina (NEB). The libraries were amplified twice using PCR amplification. DNA fragments (160–300 bp) were purified using E-gel, and sequenced on an Illumina Novaseq 6000 platform with single-end 50 bp reads.

DNase-seq data analysis

DNase-seq reads were mapped to the pig genome (Sscrofa 11.1) using Bowtie2, allowing a single mismatch. Unmapped reads were iteratively aligned by trimming 5 bp each pass until reads were less than 26 bp. Redundant reads were removed if it mapped to the same location with the same orientation. The DHSs of each library were identified using MACS2, with a cutoff of P < 1 × 10⁻⁵. The differential DHSs were scored using DEGseq an adjusted P value of P < 0.01, and then assigned to genes using HOMER. Heat maps of the DNase-seq data were generated by deepTools2.

ChIP-seq and data analysis

The frozen SATs and VATs (with two biological replicates) were grounded into powder in liquid nitrogen. The powders were fixed in 1% formaldehyde (Sigma-Aldrich) at room temperature for 20 min, and the reaction was subsequently quenched by adding glycine to a final concentration of 0.2 M and incubating for 5 min. Nuclear lysates were sonicated using Covaris S220 (Adaptive Focused Acoustics^®). Fixed chromatin was then incubated with an antibody against H3K27ac (Abcam) or M-280 sheep anti-rabbit IgG Dynabeads (Thermo Fisher) at 4 °C for 4 h. DNA was extracted using antibody–bead complexes, and ChIP-seq libraries were prepared according to standard Illumina protocols and further sequenced on an Illumina Novaseq 6000 platform with paired-end 150 bp reads. Reads were aligned to the Sscrofa 11.1 genome using Bowtie2 and filtered to remove duplicates with Picard and SAMTools. Peak calling was conducted using MACS2. Peaks were called if they passed a false discovery rate of 0.01 and were enriched over input. The reproducible peaks in two biological replicates were identified using irreproducible discovery rate (IDR) with a cutoff of P-value < 0.01. The differential analysis of peaks was performed using DESeq2 with a cutoff of P < 0.01. The differential peaks were assigned to genes using HOMER. Heat maps of the ChIP-seq reads were generated by deepTools2.

RNA-seq and data analysis

The total RNA of SATs and VATs (with two biological replicates) were isolated using the miRNeasy Mini Kit (Qiagen). Sequencing libraries were constructed using the TruSeq Stranded Total RNA Library Prep Kit (Illumina), and further sequenced on an Illumina Novaseq 6000 platform with paired-end 150 bp reads. Clean data were obtained by removing reads containing adapters, reads containing over 10% of poly-N (unrecognized bases), and low-quality reads (>50% of bases with Phred scores <5). Reads were mapped to the pig genome (Sscrofa 11.1) using Bowtie2. Gene expression levels were evaluated by the RPM (reads per million mapped reads) values of mRNA calculated using featureCounts. Differential gene expression analysis was performed using edgeR 3.31. The differentially expressed genes with a P-value of < 0.01 were identified. The R package RNASeqPower (http://bioconductor.org/packages/release/bioc/html/RNASeqPower.html) was used to estimate the power of differential expression analysis.

Gene ontology enrichment analysis and motif analysis

Gene Ontology (GO) enrichment analyses were performed with Metascape (http://metascape.org/gp/index.html#/main/step1). Enrichment of known motifs in differential DHSs was analyzed using Homer (http://homer.ucsd.edu/homer/motif/).

Establishing chromatin accessibility landscape of S-SVF and V-SVF

To understand the characteristics of chromatin accessibility landscape in S-SVF and V-SVF, we isolated SVF from porcine SAT and VAT and used liDNase-seq method (pooling SVF form two biological replicates) to map DHSs in S-SVF and V-SVF. It is known that the RNA polymerase II complex binds double-stranded DNA and releases single-stranded DNA around the transcriptional start site (TSS) to establish a DNase I hypersensitive region. This region represents open chromatin, and therefore, would have abundant DHS signal intensity (Allen & Taatjes, 2015; Heinz et al., 2015). Thus, we examined the DHS signal intensity around the gene TSS within ± 2 kb. We found strong DHS signal intensity around the TSSs of numerous genes in both S-SVF and V-SVF (Fig. 1A). Using peak calling, 3930 and 2801 DHSs were identified in S-SVF and V-SVF, respectively (Table S1). To characterize these DHSs, we examined their genomic location on annotated genes (Fig. 1B and Table S2). The DHSs are highly enriched in proximal promoters (−2 kb to +100 bp from TSS) in S-SVF (50.9%). This ratio decreased in V-SVF (21.6%), and DHSs at gene bodies (exon and intron) and intergenic regions are increased in V-SVF, showing the different distribution pattern of DHSs between S-SVF and V-SVF. Furthermore, we found that CD34 and ALCAM (CD166), which are the molecular markers of SVF (Bora & Majumdar, 2017; Walmsley et al., 2015), have obvious DHS in distal gene promoters and TSS (Fig. 1C), indicating that our isolation of SVF and corresponding DNase-seq data are reliable.

Heterogeneity in chromatin accessibility landscape of S-SVF and V-SVF

To determine how the chromatin accessibility landscape changes between S-SVF and V-SVF, we performed differential analysis of DHS signal intensity between S-SVF and V-SVF. We identified 1261 significantly different DHSs (P < 0.01) (Fig. 2A and Table S3-1), suggesting distinct chromatin accessibility landscapes between S-SVF and V-SVF. Differential DHSs are mainly located at gene bodies (42.2%) and proximal promoter regions (34.2%) (Fig. 2B). To determine the specific genes that were associated with the DHSs, we assigned DHSs to their closest genes if they were located in a gene body or promoter (−20 kb to +100 bp from TSS) region. One thousand sixty-one differential DHSs were assigned to 790 genes (Table S3-2). By using GO enrichment analysis, we found that genes with differential DHSs were enriched in functional categories related to neuronal development, metabolism of lipids and carbohydrates, immune processes, and inflammation (Fig. 2C). These biological processes are closely related to known differences in metabolic effects, and inflammatory potentials between SAT and VAT (Tchkonia et al., 2013; Tran & Kahn, 2010). For example, preadipocytes in SAT tend to have greater capacity for lipid accumulation than VAT (Tchkonia et al., 2006), and our results determined that some ELOVL fatty acid elongase family genes, such as ELOVL4 and ELOVL6, which are crucial to regulating lipid biosynthesis and insulin sensitivity (Matsuzaka et al., 2007; Shimano, 2012), had differential DHSs at their gene loci. Specifically, ELOVL4 exhibited a higher DHS signal intensity at gene promoter in S-SVF than that in V-SVF (Fig. 2D). This differential DHSs in the ELOVL family genes implies differences in gene transcriptional activity, which in turn, may cause different expression levels of ELOVL family genes, and may directly influence the capacity of lipid accumulation between S-SVF and V-SVF. The promoter of phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit beta (PIK3Cβ), which encodes an isoform of the catalytic subunit of phosphoinositide 3-kinase (PI3K), also shows strong DHS in S-SVF but not in V-SVF (Fig. 2D), indicating that the PIK3Cβ gene may have a stronger transcription activity in S-SVF than V-SVF. These data suggest that there is a reinforced signal transduction of the insulin-induced PI3K pathway in S-SVF, which is consistent with previous studies that found that SAT is associated with improved insulin sensitivity, while VAT is associated with insulin resistance (Lafontan & Berlan, 2003; Zierath et al., 1998). Moreover, the intron 1 of an adaptor protein, phosphotyrosine interacting with PH domain and leucine zipper 2 (APPL2), which is involved in impairing insulin sensitivity and inflammation (Cheng et al., 2014; Yeo et al., 2016), exhibited a strong DHS in V-SVF but not in S-SVF (Fig. 2D). This difference in DHSs in APPL2 may result in differential expression of APPL2, thereby influencing the regulation of APPL2, and subsequently, influence changes in insulin sensitivity and inflammation in S-SVF compared with V-SVF.

As DHSs are always bound by transcription factors (TFs), the differential DHSs could facilitate identification of some key TFs that may potentially regulate the differential DHS-associated genes. Therefore, we performed motif enrichment analysis in these differential DHSs. The results found motifs for the Krüppel-like factor (KLF) family of TFs, which have been reported to play a crucial role in regulating adipocyte metabolism and differentiation (Brey et al., 2009; Hsieh et al., 2019), were significantly enriched (Fig. 2E), suggesting different developmental traits of preadipocytes between SAT and VAT.

Linking accessible DHSs to H3K27Ac histone modification

Open chromatin structure, as indicated in our study by DHSs, are often associated with histone modifications (Heintzman et al., 2007; Wang et al., 2008). The histone modification H3K27Ac is a mark of transcription activation in mammals (Calo & Wysocka , 2013). To investigate whether our DNase-seq DHSs colocalized with the H3K27Ac histone modification, we performed H3K27Ac ChIP-seq using porcine SAT and VAT with two biological replicates, peak calling each ChIP-seq data. Reproducible peaks for the two biological replicates were identified by using the IDR, with 29,950 and 36,031 reproducible peaks identified in SAT (Table S4-1) and VAT (Table S4-2), respectively (IDR <0.01). Next, we examined the relative ChIP-seq signal intensity within 2 kb of the center of DHSs in SAT and VAT. The results found that H3K27ac signals were enriched around the DNaseq-seq DHSs (Fig. 3A). Moreover, we intersected the DNase-seq DHSs with H3k27a peaks and found that 65.4% and 27.3% of all DHSs overlap with H3k27a peaks (overlap >10 bp) in SAT and VAT, respectively (Fig. 3B and Table S5). Notably, the DHSs in gene proximal promoter frequently overlapped with H3k27a peaks in both SAT and VAT (94.9% and 85.5%), indicating that the DHSs in proximal promoters are more likely regulating transcription activation. However, the DHSs sites in intergenic regions rarely overlapped with H3k27a peaks (Fig. 3B). Together, these data suggest that there is a strong correlation between DNaseq-seq DHSs and H3K27ac sites, especially the DHSs in gene proximal promoters, and these data suggest a significant association between the transcriptional activation and H3K27ac histone modification.

To determine whether SAT and VAT have differential H3K27Ac histone modification patterns, we performed differential analysis of H3K27Ac signal intensity between SAT and VAT. Three thousand one hundred sixty-two differential peaks were identified between SAT and VAT (P < 0.01) (Fig. 3C and Table S4-3). To gain insight into the genes regulated by differential peaks of H3K27ac, we assigned these peaks to annotated genes using the same criterion as those used in DNase-Seq. One thousand three hundred fifty-nine peak-associated genes were identified (Table S4-4). We intersected this gene dataset with the differential DHS-associated gene dataset to identify a set of 104 genes that contain both differential DHSs and differential H3K27ac peaks (Fig. 3D and Table S4-5). GO enrichment analysis on these 104 genes were mainly enriched in categories related to energy metabolism of adipocytes, immune processes, inflammation, and neuronal development (Fig. 3E). Notably, an accessible region in the TSS of pyruvate dehydrogenase phosphatase regulatory subunit (PDPR), which encodes the regulatory subunit of pyruvate dehydrogenase phosphatase (PDP), had stronger signal intensity in S-SVF. Consistently, a higher H3K27ac peak was observed at same site in SAT (Fig. 3F), suggesting that the DHS was modified by H3K27ac, and this region may act as an enhancer to promote PDPR gene transcription in S-SVF. PDP is a crucial activator for pyruvate dehydrogenase complex (PDC),which plays a central role in linking glycolysis pathway to the oxidation process of tricarboxylic acid cycle and fatty acid synthesis (Patel & Korotchkina, 2006). Thus, the different transcriptional regulation of the PDPR gene may cause differential activation of PDC, further influencing the differences in the energy metabolism processes between SAT and VAT.

Linking differentially expressed genes to differential DHSs

Gene expression regulation is often influenced by multiple cis-regulatory elements. To investigate the relationship between differential DHSs identified by DNase-seq and differentially expressed genes, we performed RNA-seq on porcine SAT and VAT with two biological replicates. The detailed information and sequencing depth of RNA-seq are listed in Table S6. Differential analysis of global gene expression profiles between SAT and VAT was performed, and 1,784 differentially expressed genes (DEGs) were identified (P < 0.01), including 736 up-regulated genes in SAT and 1,048 up-regulated genes in VAT (Fig. 4A and Table S7). The results of a power analysis calculation shown that the estimated power is 8.12 with two biological repeats in RNA-seq. We compared the DHS signal intensity distribution in DEGs and non-DEGs within 5 kb upstream and downstream in both S-SVF and V-SVF. We found that the DHS signal intensity of DEGs were higher than non-DEGs throughout promoter, gene body, and downstream regions of gene transcription end site (TES) in both S-SVF and V-SVF, especially around the TSS. Specifically, the DHS signals of DEGs is 24.9% and 20.7% higher than that of non-DEGs in S-SVF and V-SVF, respectively (Fig. 4B), indicating that the chromatin around TSS of DEGs has a higher accessibility and contains more cis-regulatory elements than non-DEGs. Thus, we concluded that differential expression in the DEGs was more likely to be explained by the different chromatin accessibility in these DEGs between SAT and VAT.

To further explore key functional genes that could cause differential biological characteristics between SAT and VAT, we intersected the DEG dataset with differential DHS-associated gene dataset and differential H3K27ac peaks-associated gene dataset to identify a set of 20 genes, which showed corresponding differential DHSs, differential H3K27ac peaks, and differential gene expression (Fig. 4C, Fig. S1 and Table S8). GO enrichment analysis of these 20 genes was enriched in immune processes and lipid metabolism (Fig. 4D). Notably, two energy metabolism associated genes were identified, nuclear receptor subfamily 1 group D member 1 (NR1D1) and Crystallin Mu (CRYM) (Fig. 4C, highlighted in red). NR1D1 is a key regulatory gene of adipogenesis and lipid metabolism (Marciano et al., 2014). The DNase-seq data found two DHSs, one at the proximal promoter and one in intron 1, in S-SVF that were not present in V-SVF. Moreover, the H3K27ac ChIP-seq data showed there were strong H3K27ac peaks across the proximal promoter, TSS, and part of intron 1 of NR1D1 in SAT but not in VAT. Simultaneously, the RNA-seq data found higher expression levels of NR1D1 in SAT than VAT (Fig. 4E). CRYM has been reported to regulate insulin sensitivity in human WAT (Serrano et al., 2014). The proximal promoter of CRYM showed a strong DHS in S-SVF but not in V-SVF, at the same location of genome. Additionally, SAT had a strong H3K27ac signal but VAT did not, while RNA-seq data showed that CRYM was highly expressed in SAT, but was almost not expressed in VAT (Fig. 4E). The higher expression of NR1D1 and CRYM in SAT could be explained by higher chromatin accessibility and H3K27ac signal around gene TSS in SAT. Taken together, chromatin accessibility could be a key regulator of gene expression in SAT and VAT, and these differences could lead to the distinct lipid metabolism and immune processes present in these two adipose tissue types.