Subcutaneous Adipose Tissue Transcriptome Highlights Specific Expression Profiles in Severe Pediatric Obesity: A Pilot Study

The prevalence of pediatric obesity is rising rapidly worldwide, and “omic” approaches are helpful in investigating the molecular pathophysiology of obesity. This work aims to identify transcriptional differences in the subcutaneous adipose tissue (scAT) of children with overweight (OW), obesity (OB), or severe obesity (SV) compared with those of normal weight (NW). Periumbilical scAT biopsies were collected from 20 male children aged 1–12 years. The children were stratified into the following four groups according to their BMI z-scores: SV, OB, OW, and NW. scAT RNA-Seq analyses were performed, and a differential expression analysis was conducted using the DESeq2 R package. A pathways analysis was performed to gain biological insights into gene expression. Our data highlight the significant deregulation in both coding and non-coding transcripts in the SV group when compared with the NW, OW, and OB groups. A KEGG pathway analysis showed that coding transcripts were mainly involved in lipid metabolism. A GSEA analysis revealed the upregulation of lipid degradation and metabolism in SV vs. OB and SV vs. OW. Bioenergetic processes and the catabolism of branched-chain amino acids were upregulated in SV compared with OB, OW, and NW. In conclusion, we report for the first time that a significant transcriptional deregulation occurs in the periumbilical scAT of children with severe obesity compared with those of normal weight or those with overweight or mild obesity.


Introduction
The adipose tissue (AT) is a heterogeneous endocrine organ consisting of several depots distributed throughout the human body [1]. The AT localized in the connective tissue under the skin is termed subcutaneous AT (scAT), while the fat that surrounds the internal organs is the visceral adipose tissue (VAT) [2]. However, the heterogeneity of AT is not limited to its distribution as it also presents differences at the cellular level which ultimately lead to differences in function [1]. In term of function and lifespan, it is possible to differentiate three types of AT: brown adipose tissue (BAT), beige or brite -Triglyceride-glucose index (TyG index), evaluated using the formula (ln[fasting triglycerides (mg/dL) × fasting plasma glucose (mg/dL)/2]) [23].

RNA Isolation and Library Preparation
Total RNA was extracted from the SAT using the Trizol reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's instructions. RNA quantity and integrity were estimated using the Qubit instrument (Invitrogen, Carlsbad, CA, USA) and the Agilent 4200 TapeStation System (Agilent, Santa Clara, CA, USA). RNA libraries were prepared using the CORALL total RNA-Seq library Prep Kit (Lexogen, Vienna, Austria), and the RiboCop rRNA Depletion Kit (Lexogen, Vienna, Austria) was used to remove rRNA. Library quality was assessed with a High Sensitivity D1000 ScreenTape Assay using the 4200 TapeStation System (Agilent, Santa Clara, CA, USA) and quantified using a Qubit ds-DNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA). Libraries were sequenced using 75 bp paired-end reads on the Illumina NextSeq 500 platform (Illumina, San Diego, CA, USA).

Bioinformatic Analysis and Quality Assessment of Raw Data
The quality of the raw data output was examined on a FastQC (last accessed on 31 May 2022). The bioinformatic data analysis pipeline processed with FASTQ was generated using an Illumina NextSeq sequencer via unique molecular identifier (UMI) extraction, trimming, alignment, and quality control steps. Because CORALL libraries contain N12 UMIs at the start of Read 1, UMIs were removed in the first step using the UMI tools software. Adapter sequences, poly(A) sequences at the 3 end of Read 1, and poly(T) sequences at the 5 end of Read 2 were then trimmed using the Cutadapt software. Subsequent steps to assess gene and transcript intensities were carried out using the STAR software.
For each sequenced sample, the total number of reads sequenced, the number of reads mapped on the GRCh38, and the percentage of reads on target are reported in Supplementary Table S1. The four samples (OB_12, OB_15, NW_22, and NW_23) with a percentage of mapped reads below 70% were excluded from further bioinformatic analyses. Gene and transcript abundance were computed using the FeatureCounts software, with the "stranded forward" option.
A differential expression analysis was performed using R package DESeq.2; coding and non-coding genes were considered differentially expressed and retained for further analysis when |log2 group2/group1)| ≥ 1 and FDR ≤ 0.1 [25]. The raw data obtained from the RNA-seq analysis are deposited in the Gene Expression Omnibus repository with the accession number GSE228892.
The R software was used to generate heatmaps (the heatmap.2 function from the R ggplots package), PCA plots (the prcomp function from the R ggplots package), and volcano plots. An enrichment analysis of the differentially expressed genes (DEGs) was performed using g:Profiler. A gene ontology and functional enrichment analysis of the DEGs was carried out using the webtools Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.ad.jp/kegg, accessed on 15 December 2022), Reactome, and WikiPathways. A protein-protein interaction (PPI) analysis was performed using STRING. To evaluate global changes in gene expression, a gene set enrichment analysis (GSEA) was carried out using iDEP (integrated differential expression and pathway analysis, http://bioinformatics.sdstate.edu/idep96/, accessed on 15 December 2022).

Results
Periumbilical subcutaneous adipose tissues (scATs) were collected from Caucasian male children to identify gene expression patterns associated with obesity. All samples were stratified according to the BMI z-scores of the subjects.

Clinical and Biochemical Features of Enrolled Subjects
The clinical and biochemical parameters of the subjects whose AT was included in the RNA-seq analysis are reported in Table 1. PCA plots summarizing the global variability in gene expression levels allowed us to distinguish specific clusters between OW and NW ( Figure S1A), while no difference in cluster distribution was found between OB and OW ( Figure S1B) or between SV and NW ( Figure 1A). In the SV vs. OB ( Figure 1B) and SV vs. OW ( Figure 1C) comparisons, specific clusters were observed. The hierarchical clustering heatmaps show the 60 transcripts most deregulated between SV and NW ( Figure 1D), SV and OB ( Figure 1E), and SV and OW ( Figure 1F). The comparison between OB and NW showed only 1 downregulated DEG: zinc finger and AT-hook domain containing (ZFAT) ( Table S2). In total, three DEGs, including two protein-coding (one downregulated (ZFAT) and one upregulated (odd-skipped related transcription factor 1, OSR1)) and one downregulated lncRNA, were obtained from the comparison between OW and NW (Table S2). There were 178 DEGs between SV and NW, including 166 coding genes (83 downregulated and 83 upregulated) and 12 non-coding genes (4 downregulated and 8 upregulated), as is shown in Figure 2A. Different non-coding biotypes were present, and, specifically, two were long non-coding RNAs (lncRNAs), two were processed pseudogenes, one was a transcribed_unprocessed pseudogene, two were small nuclear RNAs (snRNAs), and four were small nucleolar RNAs (snoRNAs) ( Figure 2B). In Figure 2C, the 12 most downregulated and upregulated DEGs in the SV vs. NW comparison are reported. The comparison between OB and OW revealed four downregulated genes, (two coding and two non-coding genes). Of the two non-coding genes, one was a processed pseudogene and one was an unprocessed pseudogene. As is shown in Figure 2C, out of 537 total DEGs between SV and OB, 503 were coding genes (333 downregulated and 170 upregulated) and 34 were non-coding genes (10 downregulated and 24 upregulated). Regarding the non-coding biotypes, 12 were lncRNAs, 3 were processed pseudogenes, 1 was a transcribed_unprocessed pseudogene, 4 were snRNAs, 9 were snRNAs, 1 was a TEC, and 4 were misc_RNAs ( Figure 2E). In Figure 2F, the 12 most downregulated and upregulated DEGs in the SV vs. OB comparison are shown. The comparison between SV and OW showed that 762 out of the 796 DEGs were coding genes (432 downregulated and 330 upregulated) and 44 were non-coding genes (11 downregulated and 33 upregulated) ( Figure 2G). Regarding the non-coding biotypes reported in Figure 2H, 12 were lncRNAs, 8 were processed pseudogenes, 1 was a transcribed_unprocessed pseudogene, 1 was an rRNA_pseudogene, 7 were snRNAs, 10 were snoRNAs, 1 was a TEC, and 4 were misc_RNAs. In Figure 2I, the 12 most downregulated and upregulated DEGs in the SV vs. OW comparison are reported.
When focusing on lncRNAs, the analysis showed the presence of upregulated and downregulated lncRNAs amongst the different comparisons. Specifically, the SV vs. NW comparison revealed one upregulated lncRNA (OIP5-AS1) and one novel downregulated lncRNA (ENSG00000285756). The SV vs. OB comparison showed six downregulated lncRNAs, of which two had never been reported before (ENSG00000285756 and ENSG00000260267), and five upregulated lncRNAs (ENSG00000261468 and ENSG00000235609 were novel). Finally, in the SV vs. OW comparison, four downregulated lncRNAs were found, of which two were novel (ENSG00000282057 and ENSG00000285756), and seven upregulated lncRNAs were found, of which three were novel (ENSG00000272335, ENSG00000261468, and ENSG00000235609).
Notably, OIP5-AS1 and the novel ENSG00000285756 lncRNA were found to be shared in all three comparisons. In contrast, SNHG5, MAP3K4-AS1, and the novel ENSG00000285756, ENSG00000261468, and ENSG00000235609 lncRNAs were found in both the SV vs. OB and SV vs. OW comparisons. Novel lncRNA DEGs are reported in Table 2.

Functional Enrichment Analysis
To gain insight into the biological pathways, a KEGG pathway analysis ( Figure 3) together with Reactome and WikiPathways analyses ( Figure S2) were performed, comparing the NW, OB, and OW groups with the SV group. The pathways involved in the metabolism of different macromolecules, such as sphingolipids, lipids, and proteins, were found to be deregulated. In addition, the hippo, oxytocin, and apelin signaling pathways were found to be deregulated, as well as the cardiomyopathy pathways ( Figure 3). In the PPI analyses, relevant significant networks were only detected among the deregulated genes in the SV vs. OW comparison. In particular, 762 nodes and 1885 edges were predicted, and the GO was significantly enriched, resulting in the primary metabolic process (biological process; GO:0044238; Figure S3). A KEGG pathway analysis was also conducted for the lncRNAs using the ncPath tool, showing their likely roles in several biological processes. Of particular interest is their involvement in metabolic pathways or pathways linked to events potentially related to excess body weight, e.g., the regulation of lipolysis in adipocytes, the insulin signaling pathway, insulin resistance, type II diabetes mellitus, the GnRH signaling pathway, the estrogen signaling pathway, the thyroid hormone signaling pathway, the TNF signaling pathway, or the IL-17 signaling pathway.

Gene Set Enrichment Analysis
Further examination of the RNA sequencing data was carried out using GSEA for GO and pathway analysis in order to evaluate global gene expression changes. The downregulated and upregulated pathways shared between the comparisons are reported in Table 3. The lipid metabolism pathways were altered in the SV group when compared with the OB and OW groups. In particular, fatty acid degradation and metabolism were upregulated in SV vs. OB and SV vs. OW, while the regulation of lipolysis in adipocytes was downregulated in OB vs. NW and upregulated in SV vs. OB ( Figure 4A,B). Moreover, bioenergetic processes, such as pyruvate metabolism, the TCA cycle, and oxidative phosphorylation, were upregulated in the SV group when compared with the other three groups. In addition, the upregulation of branched-chain amino acid (BCAA) catabolism was observed. Additionally, the NAFLD and myocardiopathy pathways were found to be deregulated in children with severe obesity.  red boxed indicate upregulated genes; gray boxes indicate genes whose expression is not affected in these conditions.

Regulation of Lipolysis in Adipocytes
Regarding the regulation of lipolysis in the adipocyte pathway, the KEGG analysis showed that almost all the proteins involved in this pathway were downregulated (shown in green) in the OB vs. NW comparison ( Figure 4A), except for neuropeptide Y (NPY) and protein kinase G (PKG), which were upregulated (shown in red) ( Figure 4A). In contrast, in the SV vs. OB comparison, the latter proteins were downregulated (shown in green), together with adenylate cyclase (AC), insulin receptor substrate (IRS), and Akt, while all the others were upregulated (shown in red) ( Figure 4B).

Obesity-Associated Diseases
Next, we investigated the gene expression changes in children with severe obesity using the OMIM database in order to assess their possible association with mendelian diseases. Between the SV vs. OB and SV vs.

Discussion
Childhood obesity represents a serious challenge for global public health [8], and it may also contribute to obesity in adulthood, leading to negative health outcomes. It is known that during obesity, dysfunctional AT undergoes immune, metabolic, and functional changes [26]. To date, microarray analyses of blood cells, isolated adipocytes, adipose tissue, or stromal vascular fractions make up nearly all published studies, and, in some cases, children are grouped together without considering the differences between normal weight and obesity [27]. However, a few RNA-seq studies have been performed on the AT of children or adolescents [28,29]. In their work, Sheldon and colleagues compared intra-abdominal AT collected from severely obese adolescents in relation to different stages of NAFLD [28]. In addition, a pilot study by Zingale et al. focused on the expression of neuro-inflammatory markers, proposing a relationship between obesity and neuroinflammation [29].
Here, we explored for the first time the transcriptomic profile of scATs collected from the periumbilical area in four different groups of pediatric subjects, showing differences in expression according to the severity of ponderal excess and using BMI as a reliable indicator of body fatness. However, since BMI scores vary with age and sex [20], it was necessary to use BMI z-scores, i.e., BMI normalizations with respect to the two above indicated parameters. Moreover, in our cohort, only prepubertal male children were considered in order to limit the impact of sex hormones, sex chromosome complements, or developmental hormonal variation. It has been demonstrated that during puberty, males and females accumulate different types of lipids [30] in different body regions [31]. In addition, ethnicity can also have an impact on body fat distribution [32], and for this reason our patients were all of European origin.
It is well known that blood levels of triglycerides, cholesterol, glucose, and insulin greatly differ between children with and without obesity [33], rendering these parameters good biomarkers for metabolic derangement. In line with this, blood tests in our patients revealed a pronounced increase in fasting glycemia, insulin, and triglycerides from NW to SV, showing that the differences were already noticeable between NW and OW or OB. However, such changes observed at the plasmatic level were not reflected in the AT gene expression. In fact, no differences were found in the periumbilical scATs of children with either obesity or overweight when compared with those of normal weight, as is shown by the PCA and heatmap plots. Thus, we hypothesize that this could be due to the fact that this type of adipose tissue is less metabolically active than visceral adipose tissue.
Thus, it could be interesting to investigate visceral transcriptomic profiles, and to increase the sample size. The most marked differences were observed in the comparisons with the SV group, as demonstrated by the hundreds of coding and non-coding DEGs that were found. The pathway analysis of the coding DEGs revealed their associations with several signaling cascades, including the metabolic, hippo signaling, oxytocin, apelin, and HIF-1α pathways, proving that our data are consistent with previous observations reported in the literature [34][35][36][37]. A recent RNA-seq report on periumbilical scAT from children with obesity and normal weight revealed a significant deregulation in several coding DEGs implied in fatty acid and carbohydrate metabolisms, as well as in inflammatory pathways [29]. In our study, we did not observe any differences between the OB and NW groups. Moreover, even though inflammation is a characteristic feature of obesity [7], no inflammatory pathway was found in any of our comparisons. This may be due to the inflammation likely present in NW subjects undergoing surgery. Moreover, we report a significant variation in non-coding DEGs, highlighting the important impact of epigenetic regulation on obesity. Besides coding genes, a number of non-coding DEGs were also found to be deregulated in the SV group compared with NW, OW, and OB groups. Notably, OIP5-AS1 was found to be upregulated in the SV group when compared with the other groups. Although it has been implicated in a wide variety of cellular processes, as well as in chronic diseases (e.g., diabetes, myocardial ischemia) and the progression of several cancers, we suggest that OIP5-AS1 plays a further role in pediatric obesity [38]. However, further studies are needed to explore the molecular mechanisms involved in its association with obesity, as well as to consolidate the evidence for its potential use in distinguishing children with severe obesity from NW, OW, and OB children.
It is worth emphasizing that even if we notice a significant deregulation between SV and NW, this difference is more pronounced between SV and OB or OW. Obesity is a multifactorial disease, and adipose tissue expansion is a dynamic process that may already have begun during intrauterine life [39]. Adipose tissue alterations depend on several factors, such as eating behaviors, dietary components, the gut microbiome, drugs, and physical exercise [40][41][42]. All of these factors can have epigenetic consequences which impact the adipose tissue gene expression. Moreover, lipids have also been associated with epigenetic regulation, since they can regulate chromatin structure and stability [43,44]. To elucidate the reasons behind the differences in diversity among the groups, it could be very interesting and helpful to study lipidomic profiles.
Furthermore, to achieve a global overview of transcriptional alterations, we compared our data with a priori gene sets. In the SV group, three biological processes were found to be deregulated in comparisons with the other three groups considered in this study: lipid metabolism, bioenergetic processes, and pathways related to CVD. Regarding lipid metabolism, it is interesting to note that while the regulation of lipolysis has been observed to be upregulated in biopsies of the SV group, this pathway was downregulated in the comparison between the OB and NW groups. Insulin was associated with fatty acid synthesis and lipolysis inhibition [45]. However, insulin resistance promotes lipolysis [46]. In line with this, the elevated HOMA-IR index obtained from the SV group seems to correlate with the increase in the lipolysis pathway. However, in the OB group, the lipolysis regulation trend is inverted compared with that of the NW group. It has been reported that TNF-alpha-mediated acute inflammation could accelerate lipolysis [47]. It is important to consider that the NW patients enrolled in this study were subjected to emergency chirurgical interventions, and their results were thus possibly obtained under acute inflammatory conditions. Hence, we suppose that, even though insulin resistance is present in NW children, inflammation could have a major effect on lipolysis regulation. However, more evidence is needed to confirm this supposition. Regarding the bioenergetic processes, it is well known that pyruvate metabolism, the tricarboxylic acid (TCA) cycle, and oxidative phosphorylation are linked together. In fact, the metabolism of pyruvate supports the citric acid cycle thanks to its conversion to acetyl-CoA. The reducing equivalent NADH produced in the TCA is subsequently re-oxidized back into NAD+ in the electron transport chain (ETC) or oxidative phosphorylation, coupling this process with the export of protons across the inner mitochondrial membrane [48]. However, it has been established that these mechanisms are impaired in patients with obesity [49]. For instance, Sohn and colleagues demonstrated that TCA metabolites were higher in children affected by obesity before 18 months of weight loss [50]. The catabolism of branched-chain amino acids (BCCAs) seems to play a central role in all the lipid and energetic pathways. In fact, the degradation of valine, leucine, and isoleucine yields acetyl-CoA molecules, thus feeding the TCA and electron transport chain [51]. However, even though a relationship between BCAAs and obesity has been proved, the functional role of BCAA metabolism in the white adipose tissue of patients suffering from obesity still remains unclear and controversial [52]. In fact, while Green et al. demonstrated that BCAA catabolism increased lipogenesis in adipocyte differentiation [51], here we found an increase in lipolysis and fatty acid degradation. This difference might be due to the fact that our patients were affected by severe obesity, and where therefore in an exacerbated condition with respect to the in vitro model described by Green et al. Lastly, it is known that cardiometabolic risk factors are closely associated with obesity, not only in adulthood, but also in childhood [53]. As expected [54], the triglyceride/HDL cholesterol ratios of our patients lead us to suppose that children with severe obesity are more prone to develop CVD. On the other hand, the RNA-seq data revealed that the cardiomyopathy pathways were deregulated in the SV group, contrary to expectations [53]. However, since scAT has been reported to have fewer metabolic implications than VAT [55], we could hypothesize, in line with the findings of Liu and colleagues [56], that scAT is less associated with cardiovascular risk factors.
The restricted number of enrolled subjects could be considered the main limitation of this study, as it led to the inclusion of only a few samples for each group, in particular for the SV subjects. Even if this condition is very frequent, it is challenging to collect periumbilical subcutaneous adipose tissue from children who do not undergo abdominal surgery, and to harvest biopsies from normal-weight patients is even more challenging. For this reason, here we collected biopsies from patients who needed a surgical intervention not associated with obesity, and the same was done for the normal weight patients. The age of NW subjects could be a second limitation. Hence, we considered children from 1 to 12 years old, avoiding alterations due to breastfeeding (infants >1 year old) or due to puberty (>12 years old). Moreover, we stratified the subjects according to BMI z-score, a parameter that already accounts for the age factor. However, in the future it could be interesting to investigate the expression of the most deregulated genes in NW patients whose ages are comparable with those of the patients in the other groups.

Conclusions
In conclusion, we report for the first time that a significant transcriptional deregulation occurs in the periumbilical adipose tissue of children with severe obesity when compared with children with overweight or mild obesity, or those of normal weight, with a specific expression profile affecting lipid metabolism. Although this work is to be considered a pilot study and will need to be applied using a wider cohort, it might be useful to integrate our data with other omic approaches.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/cells12081105/s1, Figure S1: PCA plots of OW vs. NW and OB vs. OW; Figure S2: Reactome and Wikipathway enrichment pathway analysis; Figure S3: PPI analysis network in SV vs. OW comparison; Table S1: Summary statistics of reads mapping and transcripts assembly for each sample; Table S2 Informed Consent Statement: Informed written consent was obtained from the parents and/or legal representatives of the patients after they had received information about the study.

Data Availability Statement:
The raw data obtained from the RNA-Seq analysis have been deposited in the Gene Expression Omnibus repository (GSE228892).