Gene pathway development in human epicardial adipose tissue during early life

The prevalence of childhood obesity has risen sharply during the past 3 decades, resulting in more than 42 million overweight children around the world (1). Animal and cell culture models have elucidated the pathogenic mechanisms of early life obesity and associated metabolic pathways (2, 3). However, to date, clinical trials based on these findings have not been translated into successful interventions (4). A potential reason for this is that the origin of metabolic diseases is rarely a consequence of an abnormality in a single gene but instead reflects perturbations in the complex web of intragenomic regulatory activity, organ systems, and their relationships to the environment (5, 6).

Elucidating the dynamic and functional changes that are occurring in adipose tissue at different stages during development is pivotal to understanding the detrimental consequences of childhood obesity. In early life, there is a major transformation of body energy usage as heat-generating brown adipose tissue (BAT) cells are lost and are gradually replaced with white adipose tissue (WAT) (7). This process is especially important, as the magnitude in reduction of thermogenic capacity is one of strongest phenotypic associations with metabolic-related disorders (8). The effect of BAT on metabolic health primarily resides in its unique ability to rapidly produce heat and oxidize fat and glucose by the tissue-specific uncoupling protein 1 (UCP1) (9). Multilocular brown adipocytes, with cell numbers declining by age, were first identified in pericardial (presumably epicardial) fat in autopsies of infants and children (10). In adults, epicardial adipose tissue (EAT) consists of beige unilocular adipocytes, which are smaller than white adipocytes and express abundant UCP1 mRNA and protein (11). Currently, no UCP1 immunohistochemistry, gene expression profiles from microarray, or gene module analyses have been described for any adipose tissue depot during early stages of human development that are associated with congenital heart disease (CHD), although the molecular signatures have been described in a range of neonatal fat depots (12).

Following the rediscovery of BAT in human adults, there has been an extensive search for UCP1-expressing cells (13), including in EAT, which has been shown to exhibit an intermediate (or “beige”) phenotype in adults (11). However, one of the most challenging aspects in unraveling the genomic changes that occur during early adipogenesis remains the practical difficulty in obtaining relevant samples from newborns and infants (14). To address this, we collected EAT biopsies from pediatric patients through collaborations with cardiac surgeons. To analyze this data set, we used weighted gene coexpression network analyses to identify precise age-related specific changes in gene expression within EAT. Based on previous studies using animal models (15), we hypothesized that human EAT undergoes pronounced genotypic and phenotypic adaptations during the transition from the neonatal period through infancy and childhood.

Patient characteristics. EAT samples were collected from children with CHD requiring cardiac surgery for clinical indications (Supplemental Table 1; supplemental material available online with this article; doi:10.1172/jci.insight.87460DS1). For our analyses, we studied a total of 46 pediatric patients: 11 neonates (aged 6–24 days), 28 infants (aged 40 days to 12 months), and 7 children (aged 2–7 years). Physiological traits, such as weight and height, which allow calculation of Z scores for comparable parameters, such as weight for age (WAZ) and height for age (HAZ), were recorded (Table 1 and Supplemental Table 1). Growth Z scores were negative in the majority of participating children, reflecting the clinical status of the children with CHD approaching surgical intervention in our study cohort. Although the intrauterine growth of the neonates was largely unaffected by their cardiac conditions and early postnatal growth was appropriate in 14 of 28 of the infants (defined as WAZ [–1.5 < Z score]), there was greater variability in the WAZ and HAZ in children over 2 years, reflecting, in part, the inherent difficulties of maintaining adequate ongoing nutrition.

Table 1

Participant clinical and anthropometric characteristics

Histological and transcriptomic properties of EAT during the transition from the neonatal period to childhood. UCP1, the product of the primary BAT gene, was expressed in all EAT samples irrespective of age, whereas the histological appearance of UCP1 changed with age. In neonates, EAT was uniformly composed of multilocular, UCP1-positive adipocytes (Figure 1), whereas, in infants and children, it was predominantly composed of unilocular, UCP1-negative cells. Importantly, in 9 of 35 of those over 28 days of age, discrete islands of multilocular, UCP1-positive adipocytes were interspersed between white adipocytes.

Figure 1

Representative immunohistochemical detection of uncoupling protein 1 in a histological section of epicardial adipose tissue (EAT) from a neonate. (A) Light microscopy (scale bar: 400 μm) of UCP1 immunostained EAT typically observed in a neonatal sample, demonstrating an area of positively stained brown adipocytes (boxed region). (B) Light microscopy (scale bar: 25 μm) of UCP1 immunostained EAT, demonstrating the histological differences between BAT and WAT during the neonatal period.

Microarray gene expression analyses in EAT from infants (i.e., those aged 1–12 months) and children (aged 2–7 years) were similar (Supplemental Tables 2–5) but differed significantly from those of neonates (aged 0–24 days: FDR ≤ 0.05), with the difference being most marked between the neonatal and infant groups (Figure 2). Of the 1,074 differentially expressed genes in EAT during the transition to infancy (Supplemental Table 3), 103 transcripts were mitochondrial components involved in the activation of thermogenesis (P = 5.8E^–08), of which 83 were downregulated. Notably, the main marker of cell thermogenesis, UCP1, exhibited the largest reduction in gene expression with age (up to 3 log₂ fold change; Figure 2C). However, there were only 127 differentially expressed genes between the EAT biopsies from neonates and children (Figure 2I), of which 107 genes were upregulated or downregulated in the infant group (Figure 2A and Supplemental Table 4). Neither gene set included UCP1 (Figure 2G) or was significantly enriched with mitochondrial genes, indicating a possible repopulation of prothermogenic cells in this older age group.

Figure 2

Differential gene expression at distinct developmental periods in epicardial adipose tissue (EAT). (A) Venn diagram summarizing the number of differentially regulated genes between neonatal and infant EAT (red) versus that between infants and children (green). (B) Top 5 differentially expressed genes during the transition from the neonatal period to infancy (eBayes moderated t statistics) and (E) from the neonatal period to childhood. (C) Genes with the largest decline in expression between the neonatal period and infancy and (F) between neonates and children (gray background: downregulated genes with more than 2.5 log₂ fold change and FDR q > 0.05. (D) Top 5 genes with the greatest increases in gene expression between the neonatal period and infancy and (G) between the neonatal period to childhood samples. (H and I) Volcano plots showing changes in EAT gene expression between (H) neonatal period and infancy and (I) infancy and childhood (cutoffs: log fold change = 1, log P = 1.3, q < 0.05, FDR correction; eBayes moderated t statistics).

Unsupervised clustering and principal component analyses revealed that a majority of the samples from the neonatal age group form a distinct gene cluster, possibly highlighting the existence of different gene-to-gene interactions (Figure 3). Intriguingly, both cluster analyses concluded that, although there were large changes in gene expression, the neonatal and infant groups each clustered closely together but away from the majority of the samples from children. Consequently, this identified a distinct transcriptome pattern undetected by standard statistical procedures. Therefore, the conjunction of statistical and cluster analysis indicates the potential existence of transcriptional changes as well as the replacement of thermogenic cells occurring during the first year of postnatal life. These findings remain consistent with the hypothesis that, in infancy and childhood, adipose tissue passes through discrete developmental stages characterized by increased white adiposity (16) and reduced thermogenic capacity (10).

Figure 3

Comparison of gene expression analysis in epicardial adipose tissue (EAT) between age groups. (A) Heat map and unsupervised hierarchical clustering dendrogram for the top 100 differentially expressed gene transcript comparisons identified on microarray analyses of EAT samples from neonates (blue, n = 11), infants (red, n = 28), and children (green, n = 11) (average linkage, Euclidean distance metric; eBayes moderated t statistics, q < 0.05). Gene expression was transformed to a Z score. Blue indicates an increase and red a decrease in gene expression by age. Principal component analysis (PCA) of EAT gene expression patterns for (B) 3 age groups and individual interactions between (C) neonatal and infant EAT samples and (D) neonatal and child EAT samples. Every participant is represented by a sphere, and the specific age group of each participant is indicated by color coding. The PCA was performed using probes (n = 2,454), with P < 0.05 (without FDR correction).

Construction and description of the neonatal and infant EAT functional gene networks. To understand the shift of gene function in EAT during the early transition to infancy, we performed weighted gene coexpression network analysis (WGCNA). This unsupervised and unbiased algorithm allocated coexpressed transcripts to a specific gene module (17). By applying this algorithm, we increased the likelihood of identifying the activation of discrete biological pathways (17). It yielded 7 distinct clusters in neonates and 6 distinct clusters in infants (Figure 4, A and B, and Supplemental Tables 6 and 7), each having variable degrees of connectivity (k), most markedly in the neonatal networks, affecting their topography (Figure 4, A and B). Further analysis of the neonatal network demonstrated that it was organized in a distinct hierarchy in which the top 30 interconnected genes were allocated to the same module (represented by the neonatal turquoise module in Figure 4A; gene modules and associated colors are shown in Figure 5). This was characterized by a strong intramodule connectivity (kIN_average = 108.3). By contrast, connectivity patterns of the 6 modules from the infant network were much weaker, as illustrated in the infant blue module (kIN_average = 64.5). Strikingly, 76 of the 362 neonatal turquoise module genes were functional mitochondrial genes, including UCP1 (P = 2.06E^–22; Figure 5A), emphasizing the importance of mitochondrial regulation during early development.

Figure 4

Summary of the strength of intramodular connectivity for coexpressed genes within neonatal and infant epicardial adipose tissue (EAT) and the relationships between the gene modules and clinical characteristics for each age group. Pearson correlations were generated for the 1,041 genes coexpressed in neonates (n = 11) and infants (n = 28) and compared. All pairs of genes were transformed into network connection strengths (denoted by color intensity from yellow to red). Each cell represents an individual gene for either (A) neonates or (B) infants. Each cell shows the correlation coefficient (and P value after Bonferroni correction) between the module eigenvector (rows) to clinical characteristic (columns) for (C) neonates and (D) infants. Significant values (P < 0.05) are highlighted in red. Colors in cells identify gene modules. WAZ, weight for age; HAZ, height for age; WHZ, weight for height; BMI-Z, body mass index Z scores.

Figure 5

Biological processes and cellular components enriched in the epicardial adipose tissue (EAT) in each age network. Hypergeometric enrich- ment lists, using the DAVID database, of biological processes or cellular components found in significant (q < 0.05) modules of (A) neonatal and (B) infant gene networks using the 1,041 probe sets allocated to each gene network. The enrichment analysis was performed using each gene network module probe separately. Normalized enrichment score and prediction of fraction of genes affected in each term is shown (additional details, P values, and full gene lists are provided in Supplemental Table 8).

Integration of EAT network coexpression analysis and the relationship with neonatal and infant physiological outcomes. Next, we assessed the physiological relevance of each module by examining the overall correlation of the module genes with anthropometric indicators from each patient involved in our study. In the entire neonatal network, only the red module was correlated with any of the physiological traits; specifically, it was negatively correlated with BMI Z score (Figure 4C). This neonatal module contained 78 genes and included those with the lowest connectivity and cluster membership values in the network (kME_average = 0.69). The neonatal red module was an exception in this network, as other modules showed a strong interconnectivity between all the gene members. These results lead us to conclude that the collective expression and interactions between genes may play a strict role in the specific developmental physiological responses (e.g., thermogenic activation) in the neonatal age group. In contrast, half of infant modules (3 of 6) were correlated with 4 different physiological traits (Figure 4D). The two strongest correlations were negatively related to the HAZ score for the infant red (P = 0.003, R² = –0.54) and infant yellow (P = 0.001; R² = –0.59) modules. In addition, the infant red module in this network was positively related to child age (P = 0.01; R² = 0.46) and negatively related to WAZ score (P = 0.02; R² = –0.43). The infant green module correlated positively with WAZ (P = 0.03; R² = 0.41) and HAZ (P = 0.02; R² = 0.44) scores (Figure 4D). These results may demonstrate that the infant gene network captures the state-dependent interaction between genes and their expression in response to environmental factors, such as nutrition and physiological stress, such as is reflected in a negative Z score. Additionally, the lower connectivity observed in the infant network may be indicative of a specialized function within each module, thus increasing the importance of the interactions between peripheral modular and hub genes as potential translators of environmental adaptations.

Changes in the biological significance of the networks, capturing a transition from a programmed to an environmental responsive state. The specialization of cell function is an important factor in adipose tissue maturation during the transition from the neonatal period to infancy. In order to address this process, we examined each module for overrepresentation of relevant functional gene classifications (Figure 5 and Supplemental Table 8). Closer examination of the relationship between functional enrichment and network analysis indicated striking differences (Figure 6, A and B). For instance, in the neonatal group, 32 of 37 significant (q < 0.05) functional gene classifications were allocated to the neonatal turquoise module, while only 2 classifications were allocated to either the neonatal red or neonatal tan gene modules (Figure 5A). In contrast, for the infant group, with its relatively lower intramodule connectivity coefficients for intergene connections, a more even assignment of gene functional classes was possible between each module (varying from 1 to 15 classifications per module; Figure 5B). We were unable to assign a gene function in 2 of 6 modules.

Figure 6

Cytoscape pathway analysis of the change in primary networks within epicardial adipose tissue (EAT) with age and schematic summary of changes in adipose tissue composition, structure, and primary regulatory pathways. Comparative subnetwork visualization for hub gene modular topology in (A) neonates and (B) infants. Complete networks for each age group are shown in the smaller diagram in each panel. Layouts are based on “organic distribution,” and gene node shading depicts gene expression fold change. Red and green denote positive and negative gene expression, respectively. Relative numbers of direct gene interactions are indicated by node size and edge color (from positive [red lines] to negative [green lines]). These describe the extent (strength, intensity, or capacity) of interactions between two genes (weighted correlation). (C) Summary of effects on transcriptional architecture that accompany age from the neonatal period into childhood.

For the neonatal network, however, no functional category was overrepresented in 4 of 7 modules. The disproportionate distribution of gene function allocation between the modules in neonates might be attributed to the strength of their network organization. The infant gene modules exhibited a more specialized enrichment indicative of an enhanced internal genetic reorganization and transition to a more committed cellular function (i.e., lipogenesis). We therefore evaluated the effect of the gene enrichment analysis on the internal organization of both networks and applied the betweenness centrality coefficient [Cb(n)]. Remarkably, acetylation [Cb(n) = 0.453] had the highest value enrichment list during the neonatal period, while, in the infant network, there was a shift to protein phosphorylation [Cb(n) = 0.503]. The transition in the dominant roles between these different biochemical processes that occur with age is demonstrated when these features are displayed as a functional enriched gene network (Figure 6, A and B).

Next, we explored the function of the hub genes in each module (or the genes with largest number of connections to other genes, known in graph theory as “degree”) and its location in the networks. The genes with the highest degrees in the neonatal network were located on the neonatal turquoise module and included leucine zipper-EF-hand–containing transmembrane protein 1 (LETM1; kIN = 119.2). Further interrogation by the propinquity (i.e., nearness index) of gene connections to LETM1 determined that this list included genes such as solute carrier family 25 member 1 (SLC25A1; kIN = 102.6) and cytokine-induced apoptosis inhibitor 1 (CIAPIN1; kIN = 104.6), which are key mitochondrial regulatory transcripts whose critical functions expand from organelle maintenance to respiration, influencing the activation of thermogenesis (18–20). In addition, UCP1 (kME = 0.88) was allocated to the turquoise module in this age group, indicating the importance of the genetic regulation of thermogenesis immediately after birth (Figure 6A).

A similar analysis of the infant network revealed that the gene hub in this network is adaptor-related protein complex 2 (AP2S1; kIN = 68.5), located in the infant blue module. A brief inspection of the genes nearby to AP2S1 in the blue module indicated that these are mainly RNA transcripts of proteins performing similar functions (e.g., GANAB; kIN = 63.4 and COPE; kIN = 63.2). The phosphorylated protein products of these genes have crucial roles in all phases of the protein translation and degradation process, including cargo selection and formation of vesicles (21, 22).

The activation of thermogenesis is closely linked to the child’s ability to develop normally. We evaluated changes in the infant green module, including associated thermogenic and other mitochondrial genes as well as UCP1. Although we cannot exclude the potential that hypoxia affected gene expression patterns and profiles in our cohort (Table 1), we did not have sufficient the statistical power to examine this possibility. Critically, the majority of infants in the study were underperforming in growth (Table 1). Therefore, the positive relationships of the infant green module to WAZ and weight for height Z score (WHZ) (Figure 4D) suggest that infants with a higher Z score show less pronounced declines in UCP1 as well as other prothermogenic genes. Inhibition of cell thermogenesis and adipocyte differentiation occurs during periods of extreme food deprivation (23), indicating that these responses can be considered to be a normal response to environmental stress. Indeed, our analysis indicated that thermogenic genes are progressively replaced by other active genes associated with lipogenesis and tissue remodeling on the infant red module (24), such as PCOLCE2, ANG, and ADIRF. However, the overactivation of these transcripts may have long-term adverse health consequences, as these genes are also associated with the progression of cardiovascular diseases and with insulin resistance (24). Our results indicate that gene expression patterns are not uniquely linked by biochemical pathways or cellular compartment location, but that other, yet unidentified, factors play an important role in the formation of gene modules.

Gene expression patterns are also influenced by a diverse range of factors, including transcriptional binding sites, transcription factors, and mRNA prodegradative elements (25). Thus, the infant red module contained a number of important genes, such as SETD7 (kME = 0.89), TWIST2 (kME = 0.89), and EIF4EBP2 (kME = 0.8), which are transcriptional factors playing central roles in the transcriptional activation of genes involved in angiogenesis and insulin regulation (26–28). These results are in broad agreement with previous research on the physiological functions of EAT and its role in the development of cardiovascular disease (29). A widely accepted assumption is that mature EAT both acts as a store of free fatty acids for the myocardium and displays some thermogenic capacity (11). However, it can also be subjected to a wide range of pathological circumstances, varying from lipodystrophy, poor cardiac contractility, and obesity, with abnormal secretion of inflammatory and fibrotic factors and white adipocyte hypertrophy, leading to deterioration in cardiovascular function (30, 31).

Here, we demonstrate, for the first time to our knowledge, how contemporary systems biology approaches, e.g., through the construction of gene networks, can be used to elucidate the transcriptional function of EAT. This computational methodology, through the development of open source packages (e.g., WGCNA, DAVID, and Cytoscape), represents a complementary method in the understanding of adipose tissue biology. Our results demonstrate that this approach can be used to subdivide genes into pathways on the basis of the regulation of their expression. By segregating genes into functional groups, we demonstrate that changes in the topology intramodular networks can be very useful for identifying biologically important genes within these pathways. One key example is the neonatal network, which was characterized by high connectivity hub genes regulating cellular activities, particularly those associated with cell differentiation and mitochondrial function. Therefore, one of the most striking observations from this study was the intricate coordination of multiple processes implicated in the control of thermogenesis in EAT of neonates and the way in which these gene-to-gene interactions shift towards lipogenesis through infancy. Although our data in older children will need to be validated in larger studies, there appears to be a distinct thermogenic pathway, which is in a state of latency in older children. These findings raise the possibility that this pathway may be amenable to recruitment in studies of the prevention of childhood obesity and insulin resistance (32). Our analysis highlights the complexity and time-dependent nature of these adaptations, and it is apparent that simply targeting single genes is likely to be of limited benefit in reactivating thermogenesis through UCP1 in humans.

Another finding observed in the infant network is the identification of a cluster of genes (the infant red module), which are potent stimulators of vascular remodeling and metabolic dysfunction (24, 26–28). This module revealed striking correlations, both positive and negative, between patient age and negative Z scores for growth, and may ultimately link EAT dysfunction to environmental factors that may promote cardiovascular diseases. However, as we were unable to define a functional gene category, it remains to be determined whether this is only a statistically significant relationship or whether a common mechanistic association linking these genes exists. Interestingly, our analysis suggests that the EAT transcriptome in children is sensitive to longer term reduction of cardiac performance, as shown in adults with advanced coronary artery disease (31).

Despite limitations inherent in the sample size in our cohort, participants clustered into relatively homogenous age groups, allowing elucidation of distinct gene networks reflecting distinct biological entities. Therefore, our data support the interpretation of the existence of time-dependent trajectories in postnatal developmental pathways in adipose tissue. These results may reflect the unique endocrine environment during the transition from birth to early childhood (7, 10, 16, 19), which induces a complex series of events in the cells constituting EAT, including cell differentiation, proliferation, and the change in primary function between heat production and energy dissipation to energy storage.

While we have demonstrated variations in cellular heterogeneity intrinsic to EAT in the 3 age groups, representative histological examinations (e.g., Figure 1) demonstrate that the tissue samples taken span all cell types, minimizing any influence on expression analyses. As the cells composing the EAT respond to both functional and environmental challenges, they appear to apply different transcriptional regulatory mechanisms. These lead to a reduction or more controlled thermogenic activity in some cells, while others simply become energy stores (19). Such heterogeneity in genetic responses may underlie the separation of genes into modules as well as the differences observed in the topology of each network. However, as we noted in our cohort, this transition needs a precise temporal coordination of gene expression, which in turn is dependent on the nutritional environment.

In summary, the structure and architecture of EAT differs among the neonate, infant, and child, with profound regulatory effects on UCP1 (Figure 6C). EAT retained discrete islands of UCP1-positive cells persisting into childhood. These thermogenic cells exhibited substantial differences in their transcriptional architecture, which could alter gene activation at different age ranges in early life. In the neonate, EAT has limited flexibility and responsiveness. With age and maturation, the depot changes to become more responsive to environmental stimuli. Furthermore, we expect that the intensity of these changes is the product of both cellular programmable events and external environmental factors. The range of processes in EAT that we have identified as age- and growth-dependent, including acetylation, protein phosphorylation, and mitochondrial function, have the potential to regulate adipose tissue function. These unique, developmentally regulated gene interactions and metabolic-related pathways are potential targets for intervention strategies, including diet and/or pharmaceuticals designed to promote BAT function.

Tissue collection. With regulatory approval and informed consent, as detailed below, samples of adipose tissue overlying the right ventricle were collected between April 2012 and October 2013. Due to the variable aggregation and amount of EAT overlying the right ventricle, samples were collected from the region that was readily and safely accessible to sampling. EAT samples were taken immediately prior to the child undergoing cardiopulmonary bypass to eliminate any potential effects of intraoperative therapeutic cooling on EAT. All children were cared for by their families prior to surgery at thermoneutrality; notably no child’s residence in this central England location was equipped with air conditioning (the environmental temperature being 10.1^oC). External environmental temperatures prior to surgery undertaken for children with CHD were not different between age groups.

Height and weight measurements were undertaken using standardized, validated equipment in regular clinical use. Data on participant age, weight, height, medications, and other medical conditions were recorded prospectively. Anthropometric assessments were performed using the Emergency Nutrition Assessment Tool and the World Health Organization Anthro (v 3.2.2).

Immunohistochemistry. Formalin-fixed tissue samples were rehydrated, loaded into Histosette II (Simport) cassettes, and processed by ethanol dehydration followed by xylene rinse. The samples were then embedded in paraffin wax by the Shandon Excelsior tissue wax processor (Thermo Scientific). Five-micrometer tissue sections were cut from the samples using a sledge microtome (Anglia Scientific) and mounted on to Superfrost Plus slides (Menzel-Gläser Inc.). At least 2 slides per participant alongside a negative control were labeled with a random identifier, loaded into the Leica BondMax IHC slide processor (Leica Microsystem), and run on an automated software program (Vision Biosystems Bond version 3.4A) using bond polymer refine detection reagents (Leica Microsystem) and using a 1:500 dilution of preoptimized primary rabbit polyclonal antibody to UCP1 (ab10983, Abcam) for the antibody-positive samples.

RNA isolation, quantification, and quality control. For RNA isolation, 100–1,000 mg of EAT biopsy was dissected for analysis into 2 ml of TRI reagent (Sigma-Aldrich). Total RNA was extracted using the RNeasy Plus kit (Qiagen) according to the manufacturer’s instructions. Total RNA quantity was measured with a NanoDrop ND-1000 Spectrophotometer (Thermo Scientific). Optical density ratios (260 nm/280 nm) were >1.9 for all samples. Total RNA quality was assayed by the Agilent BioAnalyzer RNA 6000 Nano Kit (Agilent Technologies). Quality control required intact, distinct ribosomal peaks (RIN > 7).

Transcriptome profiling with Affymetrix GeneChip. Sample mRNA was labeled and hybridized onto Human Genome U133A plus 2 arrays according to manufacturer’s recommendations using the GeneChip 3’ IVT Express kit (Affymetrix). Detection was performed using a GeneChip Scanner 3000 7G (Affymetrix).

Expression array analysis. Normalization and network analyses were performed, unless otherwise stated, using free and open source packages from the R project (http://cran.r-project.org/). The open source Bioconductor community (http://www.bioconductor.org/) provides specific resources for microarray analysis, including the Human Genome U133A plus 2 arrays (Affymetrix). We used the function “gcrma” embedded in the “affy” package for preprocessing data stages, including background correction, normalization, and probe match verification. For statistical analysis of gene expression, we used the “limma” library, which enabled us to perform empirical Bayesian statistical modeling (33) between the 3 age groups. For all other statistical analysis, unless otherwise stated, we applied the FDR approach (34) and considered q ≤ 0.05 as a significant result. The R-package “gplots” was used for fold change, heat map, and volcano expression plots.

Accession numbers. All original microarray data were deposited in the NCBI’s Gene Expression Omnibus (accession GSE82155).

Gene coexpression network construction and analysis. The core of network construction analysis was performed using the WGCNA (17) R-package for weighted correlation network analysis as previously described (17, 35, 36), using methods developed at the David Geffen School of Medicine.

The analysis was restricted to differentially expressed genes, removing those with low expression (i.e., log₂ fold change ≥ 0.5 or ≤ –0.5; (uncorrected) P ≤ 0.01, and an average gene expression of >4.5). Each probe set was treated as a unique gene. This method yielded the two gene networks subsequently used in our analysis. We identified gene transcripts outliers by applying a canonical Euclidean distance–based network [A(uv) = 1–||S(u)–S(v)||²/maxDiss] (where A(uv) represents gene adjacency matrix values u to v; S(u), the coexpression similarity value of gene u; S(v), the coexpression similarity value of gene v, and maxDiss, maximal distant), transformed to a standardized connectivity index: A(uv) = 1–||S(u)–S(v)||^2/maxDiss [Z.k_u = (k_u – mean(k))/(sqrt(var(k)))], where Z.k_u represents the standardized connectivity index; k, the connectivity of a gene node or nodes; k_u the connectivity of gene node u; mean, the average and sqrt(var(k)), the square root of the variance (var) in k.

To distinguish between gene repression and activation, thereby retaining biologically relevant information, we implemented signed coexpressed network algorithms. To maintain the scale free topology criterion for this analysis, we choose a β (soft threshold) value of 16, where we noted that R² reached a saturation point. Then, as described by weighted correlation network analysis (WGCNA), β served as a power function for the Pearson’s correlation formula [a(ij) = |0.5 + 0.5 × cor(x(i), x(j))|β)] (where “cor(x(i), x(j)) represent the absolute Pearson correlation between gene expression profiles of genes i and j) to transform the similarity matrix (–1;1) into an adjacency matrix (0;1). The value of a(ij) indicates the specific relationship or connection strength between an individual pair of genes. The generated adjacency matrix was then applied as a measure of node dissimilarity on which the topological overlap matrix was built.

We used the dynamic branch cutting method (0.25%) to define modules as branches of the hierarchical clustering tree. To allocate the 1,041 selected genes into modules, we used “moduleManualhybrid” from the WGCNA R-package (17). By calculating the module and clinical characteristics eigenvectors, we could evaluate similarity or proximity of a specific gene to physiological trait(s) and derive module membership eigenvectors (kME) on interesting modules using Pearson correlation. Within-module connectivity (kIN) for each probe set was determined by summing the connectivity of that probe set with other individual gene members in the cluster of interest (17).

Functional annotation charts. The overexpression functional annotation analysis was performed individually for each gene module list and analyzed using DAVID (Database for Annotation, Visualization and Integrated Discovery, https://david.ncifcrf.gov/ DAVID combines results from various data-mining tools, including gene ontology, Kyoto Encyclopedia of Genes and Genomes, and SwissProt, to attribute term and fold enrichment scores for each gene list, filtered by statistical significance (i.e., FDR; q ≤ 0.05)(37).

Network visualization. For network visualization, WGCNA R gene expression profiles were imported into Cytoscape (open source v.3.0) (38). Each WGCNA file was subjected to statistical filters for fold changes in gene expression, kIN values, and gene ontology enrichment analysis. Network figures shown here use the spring-loaded organic layout from Cytoscape to aid visualization.

Statistics. For statistical analysis of gene expression, we used the “limma” library, which enabled us to perform empirical Bayesian statistical modeling (33) between the 3 age groups. For all other statistical analysis, unless otherwise stated, we applied the FDR approach (34) and considered q ≤ 0.05 as a significant result. The R-package “gplots” was used for fold change, heat map, and volcano expression plots.

Study approval. The collection of EAT samples in all participants in this study was conducted with ethical approval from the UK National Research Ethics Committee (East Midlands: Nottingham 2 REC ref: 12/EM/0214). Additionally, institutional approval was granted by the University Hospitals of Leicester NHS Trust. Informed written consent was received from parents/guardians of participating children as well as from the young person themselves if they were deemed to possess the capacity to give informed consent. Tissue samples were transferred to Academic Child Health, University of Nottingham, a facility approved under the Human Tissue Act (HTA) UK 2004. All procedures were conducted within the HTA regulatory framework.