RNA-seq Characterization of Sex-Differences in Adipose Tissue of Obesity Affected Patients: Computational Analysis of Differentially Expressed Coding and Non-Coding RNAs

Obesity is a multifactorial disease presenting sex-related differences including adipocyte functions, sex hormone effects, genetics, and metabolic inflammation. These can influence individuals’ risk for metabolic dysfunctions, with an urgent need to perform sex-based analysis to improve prevention, treatment, and rehabilitation programs. This research work is aimed at characterizing the transcriptional differences present in subcutaneous adipose tissue (SAT) of five obesity affected men versus five obesity affected women, with an additional focus on the role of long non-coding RNAs. Through RNA-sequencing, we highlighted the presence of both coding and non-coding differentially expressed RNAs, and with numerous computational analyses we identified the processes in which these genes are implicated, along with their role in co-morbidities development. We report 51 differentially expressed transcripts, 32 of which were coding genes and 19 were non-coding. Using the WGCNA R package (Weighted Correlation Network Analysis, version 1.70-3), we describe the interactions between coding and non-coding RNAs, and the non-coding RNAs association with the insurgence of specific diseases, such as cancer development, neurodegenerative diseases, and schizophrenia. In conclusion, our work highlights a specific gender sex-related transcriptional signature in the SAT of obesity affected patients.


Introduction
Obesity is a multifactorial disease that, alongside overweight, affects over 1.9 billion adults worldwide [1]. The prevalence of obesity is on the rise in many countries throughout the world; it has tripled since the 1960s, and obesity can thus currently be considered an epidemic and a major health issues [2]. Epidemiologic studies demonstrate that overweight and obesity in both men and women are associated with various aberrations including increased propensity to diabetes, cancer, and lipid profile [3,4]. Even so, there are sex-related differences in obesity including adipocyte functions, sex hormone effects, genetics, and metabolic inflammation issue, which should be considered when studying an individuals' risk for metabolic dysfunctions [3,5,6]. According to the World Health Organization, sex "refers to the biological and physiological characteristics that define men and women," whereas gender "refers to the socially constructed roles, behaviors, activities, and attributes that a given society considers appropriate for men and women" [7]. This knowledge is valuable in devising preventative measures and treatments for both men and women. Sex differences have frequently underestimated the consequences on the daily practice of medicine and in the developing specific therapeutic approaches. In fact, recent pharmacovigilance studies on large populations indicate that the risk of drug side effects in women is about 40% higher than in men during clinical trials [8]. Despite this, relatively few studies have investigated these aspects in obesity research. There is thus an urgent need to investigate sex-based analyses to improve prevention and treatment of this pathology [8][9][10].
Although men and women can both become obese, the incidence and health consequences differ between the sexes, starting from early childhood [10,11]. Moreover, they differ in genetics and in many biological processes such as fat deposition, its mobilization, and its use as a metabolic fuel and the consequences of both excess and insufficient fat stores [12][13][14]. Indeed, women usually have more subcutaneous adipose tissue (SAT) in the gluteal-femoral body area while men have more visceral adipose tissue (VAT) around the abdominal organs, resulting in two different body shapes. Despite the extensive knowledge of adipose tissue distribution, only a few studies have focused on how sex contributes to differences in adipose metabolism [15]. One work including a large cohort of 359,387 participants from nine countries confirmed that there is an association of BMI and waist-to-hip ratio (WHR), with slight differences between the sexes [16]. To support genetic differences between sexes in obesity, a meta-analysis of genome-wide association data from 38,580 individuals, followed by large-scale replication (in up to 70,689 individuals) has identified two loci (TFAP2B and MSRA) and a further locus, near LYPLAL1, suggesting gender-specific relationships with WHR and overall adiposity [17]. An additional metaanalysis of 32 genome-wide association studies for WHR adjusted for BMI identified 13 new loci in or near RSPO3, VEGFA, TBX15-WARS2, NFE2L3, GRB14, DNM3-PIGC, ITPR2-SSPN, LY86, HOXC13, ADAMTS9, ZNRF3-KREMEN1, NISCH-STAB1 and CPEB4 that exhibited a stronger effect on WHR in women than men, providing further evidence for multiple loci that modulate body fat distribution and reveal strong gene-by-sex interactions [18].
Even with an increasing number of genetic differences between women, these do not account for all observed phenotypes. Linkage analysis, especially in twins, confirmed that genetics explains 40-80% of the variance in body mass index (BMI) and in risk of obesity [19]. Indeed, there is a need to increase current knowledge, extending research efforts not only on the role of genetics, but also of epigenetics and transcriptional regulation. To investigate whether robust sex differences exist in adipose tissue gene expression, one very recent work revealed 162 differentially expressed genes, with those with an increased expression in men being associated with sex-hormone receptors, while genes that were elevated in women being implicated in oxidative phosphorylation and adipogenesis [20]. Moreover, the role of the non-coding component of the genome could play a significant role in determining differences in BMI and disease development [21]. Previously classified as "junk DNA", non-coding RNAs (ncRNAs) are a class of molecules with multiple functions in both health and disease [22][23][24][25][26]. They can be classified into subfamilies according to their size, the main ones being small-size microRNAs (miRNA/miR, <200 nucleotides) and long non-coding RNAs (lncRNAs, >200 nucleotides) [24]. The role of miRNAs in genetic sex differences in adipose tissue has recently been studied by performing high-throughput miRNA sequencing in gonadal fat depots of the Four Core Genotypes mouse model, which allowed us to assess the independent effects of gonadal type (male vs. female) and sex chromosome complement (XX vs. XY) on miRNA expression profiles [27]. Moreover, gender and obesity specific microRNAs have been found to be altered in the adipose tissue from lean and obese pigs [28].
As transcriptional regulators, lncRNAs represent emerging targets to study adipocytes development in obese patients [29,30]. Recent studies have addressed the role of lncRNAs in the network governing adipogenesis and adipocytes biology. Some lncRNAs play a key role in these phenomena by interacting with different transcription factors, such as PGC-1alpha, EBF2, ZBTB7B, ZFP516, and PRDM16 [31]. Moreover, lncRNAs could play a remarkable role in explaining gender differences. Indeed, the most notable contribution of lncRNA to sex differences is provided by the lncRNA XIST, encoded on the chromosome X [32][33][34]. Sexual dysmorphisms correlated with lncRNAs have been found for the brain, in male sex determination, and even cardiovascular diseases, but no specific evidence has to this day been reported for the adipose tissue and obesity [35][36][37]. To this end, this research work is aimed at the characterization the transcriptional differences present in the SAT of men and women, with an additional focus on the role of lncRNAs.

Adult Human Adipose Tissue Collection, Isolation and Differentiation
The present study is in accordance with the Declaration of Helsinki, and it was approved by the Ethical Committee of IRCCS Istituto Auxologico Italiano (Ethical Committee approval code #2020_10_20_04). A signed informed consent was obtained from each enrolled patient for tissue sampling. Biopsies of SAT were collected from 5 obese women (OBF; age 41 ± 12.5 years, BMI 38.2 ± 4.6 kg/m 2 ) and 5 obese men (OBM; age 42.4 ± 6.58 years, BMI 36.9 ± 3.5 kg/m 2 ), and their clinical features are reported in Table S1. Surgical biopsies of whole abdominal SAT were collected pre-operatively from obese patients during bariatric surgery procedures and from normal weight patients before aesthetic plastic surgery or abdominal surgery for non-inflammatory diseases. Each collected biopsy was weighed and stored in 1 mL of DMEM (Invitrogen Corporation, Jefferson City, MO, USA) supplemented with 2.5% Bovine Serum Albumin (Sigma, St. Louis, MO, USA) per gram of collected tissue. The biopsy was immediately transferred to the laboratory and processed. A fragment of the whole adipose tissue biopsy was immediately frozen in liquid nitrogen for RNA extraction.

SAT RNA Extraction
Approximately 500 mg of frozen SAT was homogenized in RLT buffer (Qiagen, Hilden, Germany). RNA from SAT was extracted using the RNeasy Mini Kit (Qiagen) according to the manufacturer protocol, and samples were then treated with the RNase-Free DNase Set (Qiagen). Concentration and quality of the extracted RNA were determined by the NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA) and RNA integrity verified by gel-electrophoresis.

Libraries Preparation for RNA-seq and Bioinformatic Data Analysis
RNA-seq libraries were prepared with the CORALL Total RNA-Seq Library Prep Kit (Lexogen, Vienna, Austria) using 150 ng total RNAs. The RiboCop rRNA Depletion Kit (Lexogen, Vienna, Austria) was used to remove rRNA. Qualities of sequencing libraries were assessed with D1000 ScreenTape Assay using the 4200 TapeStation System (Agilent, Santa Clara, CA, USA) and quantified with Qubit™ dsDNA HS Assay Kit (Invitrogen, Carlsbad, CA, USA). RNA processing was carried out using Illumina NextSeq 500 Sequencing. FastQ files were generated via llumina bcl2fastq2 (Version 2.17.1.14-https://support.illumina.com/downloads/bcl2fastq-conversion-software-v2-20.html, last accessed on 15 February 2021), starting from raw sequencing reads produced by the Illumina NextSeq sequencer. The quality of individual sequences were evaluated using FastQC software (version 0.11.9; https://www.bioinformatics.babraham.ac.uk/projects/ download.html#fastqc accessed on 26 March 2021) after adapter trimming with cutadapt software. Gene and transcript intensities were computed using STAR/RSEM software (version 2.7; https://code.google.com/archive/p/rna-star/ accessed on 26 March 2021) [38] using Gencode Release 27 (GRCh38) as a reference, using the "-strandness forward" option. Transcript abundance was obtained using the BlueBee ® Genomics Platform (Lexogen). Differential expression analysis for mRNA was performed using R package DESeq2 [39]. Genes were considered differentially expressed and retained for further analysis with |log2(men SAT/women SAT)| ≥ 1 and an FDR ≤ 0.1. We imposed minimum |log 2 FC| of 1 and an FDR lower than 0.1 as thresholds to differentially expressed genes. The raw data obtained from the RNA-seq analysis was deposited on the Gene Expression Omnibus repository with the accession number GSE166047.

Pathway Analysis
Gene enrichment analysis was performed on coding genes using g:Profiler, ranking terms according to their fold change and using a Bonferroni-Hochberg FDR of 0.05 as the threshold [40]. The R software was used to generate heatmaps (heatmap.2 function from the R ggplots package), PCA plots (prcomp function from the R ggplots package), Volcano plots [41], Dotplot graphs (ggplot2 library), and Pathview graphs (Pathview library [42]). The NDEx plugin [43] was used to perform analysis regarding the tissue expression, subcellular localization, and protein interaction of the differentially expressed genes [44]. The iRegulon plugin was used to identify the transcription factors possibly interacting with the DE RNAs and visualized using Cytoscape [45].

Coding and ncRNAs Co-Expression Analysis
Coding RNAs' co-expression with ncRNAs was performed using the weighted gene coexpression network analysis (WGCNA) R package (https://CRAN.R-project.org/package= WGCNA, last accessed on 25 March 2021) [46]. The soft thresholding power was chosen considering the criterion of approximate scale-free topology. Network nodes represent gene expression profiles, while undirected edges values are the pairwise correlations between gene expressions. Cytoscape software was used for network import and visualization.

Correlation Analyses
For correlation analyses, the DE RNAs were correlated with anthropometrical parameters corresponding to specific patients. For each gene, the raw counts were normalized on the raw counts of EEF2 (Eukaryotic Translation Elongation Factor 2), identified as stable housekeeping gene from the Housekeeping and Reference Transcript Atlas [47]. Prism 8 software (GraphPad Software Inc., La Jolla, CA, USA) was used for statistical analysis, assuming a p value less than 0.05 as the limit of significance.

RNA Extraction and Real-Time PCR
Real-Time PCR was performed with the CFX Connect Real-Time RT-PCR System (Bio-Rad, Hercules, CA, USA) using the SsoAdvancedTM Universal SYBR ® Green Supermix (Bio-Rad, Hercules, CA, USA) as dye. Primers were designed with NCBI's Primer-BLAST tool, and they are reported in Table S2. Gene expression was calculated using the 2 −∆∆Ct method, and 18S was used as an endogenous control. Data were expressed as mean ± SEM. The statistical analysis was performed with Student's t-test. Prism 7 software (GraphPad Software Inc., La Jolla, CA, USA) was used, assuming a p-value less than 0.05 as the limit of significance.

Deep Sequencing of RNAs Expression Profiles in SAT from Obese Women and Men Reveals Transcriptional Differences
This work was aimed at characterizing the implication of sex differences in the subcutaneous adipose tissue (SAT) of obesity affected patients, fully characterizing the underlying molecular signature that leads to a different disease profile in obese men or women. Heatmap analysis of the DE RNAs showed different expression profiles (Figure 1a), and the same can be concluded when looking at the PCA representation ( Figure 1b). Indeed, the two categories tend to group together, revealing a transcriptional dysregulation. The volcano plot reports the DE RNAs (Figure 1c), and specifically 51 DE RNAs were identified, 70% of which were coding genes and 30% of which were non-coding DE RNAs (Table 1)  The list of coding DE RNAs, ranked by their FC, is reported in Table 2. Specifically, it is possible to notice genes involved in transcriptional control, including ribosomal proteins, translation initiation factors and activators, demethylases, regulators of proteins turnover, and even splicing modulators, suggesting a profound transcriptional dysregulation. Via Real-Time PCR, we validated the expression of 3 genes with increased expression in males (TTTY15, UTY, and KDM5D, Figure 1d) and 2 genes with decreased expression (XIST, ROR2, Figure 1e) in an independent cohort of SAT samples obtained from 11 females affected by obesity and 7 males affected by obesity, supporting the evidence of the RNA-seq analysis. May function as a ubiquitin-protein or polyubiquitin hydrolase involved both in the processing of ubiquitin precursors and of ubiquitinated proteins. May therefore play an important regulatory role at the level of protein turnover by preventing degradation of proteins through the removal of conjugated ubiquitin. Triadin; contributes to the regulation of luminal Ca 2+ release via the sarcoplasmic reticulum calcium release channels RYR1 and RYR2, a key step in triggering skeletal and heart muscle contraction. Plays a role in excitation-contraction coupling in the heart and in regulating the rate of heartbeats.

ENSG00000099725
PRKY 3.67 6.11 × 10 -16 Could be the product of a pseudogene. Highly similar to PRKX in the pseudo autosomal region of the X chromosome; the transcripts specific of that gene are potential candidates for nonsense-mediated decay.  The DE RNAs transcripts for OB-M vs. OB-F (|Log 2 FC| ≥ 1) were subjected to pathways analysis using g:Profiler [40], ranking the genes in term of their |FC|. A significant dysregulation was observed in all categories analyzed, which included gene ontology (GO), with 101 significant pathways for molecular function (MF), 530 significant pathways for biological processes (BP), and 73 significant pathways for cellular component (CC), pathway analyses, which highlighted 2 significant pathways for KEGG, 82 significant pathways for Reactome, and 20 significant pathways for WikiPathways, along with 2432 regulatory motifs in DNA (TF), 24 genes targeted by miRNAs (MIRNA), and implications in protein databases such as the human protein atlas (HPA, 143), CORUM (5), and implications for human phenotypes (HP, 289) ( Figure 1d, Table S3). The iRegulon plugin was used to identify the transcription factors possibly interacting with the DE RNAs ( Figure 1e). Specifically, 47 transcription factors were identified as interacting with 25 DE RNAs. Moreover, the 5 CORUM protein binding complexes formed by the DE RNAs were represented, and these were the cytoplasmic 40S ribosomal subunit, the cytoplasmic dynein-2 complex, the Sam-68-p85 P13K-IRS-1-IR signaling complex, and 18S U11/U12 snRNP (Figure 1f).

GO Terms Enrichment Shows Numerous Differences in Genes Pertaining Specific MF, BP, and CC
The GO terms analysis in MF highlighted 101 significant pathways, and the top 10 are displayed as a GO Chord graph in Figure 2a. Interestingly, the highest perturbation seems to be due to the increased expression of two demethylating enzymes: KDM5D and UTY. Indeed, amongst the perturbated processes it is possible to notice "protein demethylase activity", "histone demethylase activity", and "demethylase activity". The pathways involving the highest number of genes appears to be "binding", further subdivided into specific binding of nucleic acid, Wnt-protein, signaling receptor, heterocyclic compound, and protein binding. ClueGO MF analysis highlights as the top significant pathways the "regulation of MAPK activity" along with "cation binding" and "RNA binding and catalytic activity" (Figure 2b). For GO BP analysis, 530 significant terms emerged, and again the top 10 are displayed as a GO Chord graph in Figure 2c. Most pathways involve metabolic processes, interesting in the obesogenic context, including terms such as "protein metabolic process", "cellular protein metabolic process", "organonitrogen compound metabolic process", "cellular macromolecule metabolic process", "macromolecule metabolic process", "nitrogen compound metabolic process", and "primary metabolic process". Protein modifications are also represented, with the rest of the top significant terms pertaining "protein dealkylation", "protein demethylation", and "histone demethylation". ClueGO analysis highlights implication for genes in chordae embryonic development, along with processes already displayed as GO Chord and, interestingly, chromatin organization and specific pathways such as that involving BMP or the MAPK cascade (Figure 2d). Lastly, 73 significant terms emerged in the GO CC analysis, and the top 10 are displayed as a GO Chord graph in Figure 2e. Interestingly, a strong implication is given to the nucleoplasm along with chromatin remodeling complexes such as the "methyltransferase complex" and the "histone methyltransferase complex". Another compartment seems to be involved in the membrane, with terms such as "membrane-enclosed lumen" and "membrane-bounded organelle". ClueGO also reports a strong implication for the "endocytic vesicles and the organelles lumen" (Figure 2f). displayed as a GO Chord graph in Figure 2e. Interestingly, a strong implication is given to the nucleoplasm along with chromatin remodeling complexes such as the "methyltransferase complex" and the "histone methyltransferase complex". Another compartment seems to be involved in the membrane, with terms such as "membrane-enclosed lumen" and "membrane-bounded organelle". ClueGO also reports a strong implication for the "endocytic vesicles and the organelles lumen" (Figure 2f).

Pathway Analysis of DE RNAs Highlights Specific Processes-Involvement
DE RNAs were subjected to pathway analysis in three well renowned databases: Reactome, WikiPathways, and KEGG; g:Profiler analysis highlighted 82 significant Reactome pathways, reported in Table S3. The top 20 pathways, ranked by their significance, are reported in Figure 3a. Interestingly, again epigenetic and transcriptional regulation seems to be the most represented terms. Indeed, chromatin remodeling is represented by "HDMs demethylase histones", "chromatin organization", and "chromatin modifying enzymes", whereas numerous steps of protein synthesis, including all the steps of translation, are involved. ClueGO analysis also highlights the implication of multiple demethylase activities, with more pathways pertaining to specific molecules' involvement and even metabolic processes pertaining, for example, to glucidic metabolism (Figure 3a). Twenty pathways were found after WikiPathways analysis, and these included, again, transcription and translational pathways, along with ribosomal subunits and even an involvement for non-coding RNAs (Figure 3b). The Wnt signaling pathway is represented multiple times, along with insulin signaling and type 2 diabetes. The non-coding and translational dysregulations are also represented after ClueGO analysis, and a strong involvement of the PI3K-Akt signaling pathway is also reported (Figure 3b). KEGG analysis with g:Profiler highlighted two significant pathways, also previously mentioned, which are the Wnt signaling pathway and ribosome involvement (Figure 3c). ClueGO KEGG analysis again highlighted the highest involvement of the PI3K-Akt signaling pathway but also implicated ribosome, RNA transport, and spliceosome as processes that could impact gene expression (Figure 3c).

Characteristics of DE RNAs in SAT of Men versus Women: Interaction, Tissue Expression and Cellular Localization
The coding DE RNAs were comprehensively characterized for their interaction, tissue expression, and subcellular localization (Figure 4). To do so, the STRING database was used to construct an interaction network of DE RNAs, where the nodes are proteins and the edges represent the predicted functional associations. It is possible to see that the proteins belong to three main networks, with MUC20 being an inter-player amongst two of them (Figure 4a). With the NDEx database, it was also possible to investigate the genes expression and cellular localization. The DE RNAs in the SAT of obese men are also expressed in a high number of other tissues. These also include neural specific tissues such as the hippocampus, the cerebellum, and the cerebral cortex; tissues implicated in the genitourinary apparatus such as the fallopian tube, the endometrium, the urinary bladder, and the seminal vesicle etc.; and also the muscle tissue (e.g., skeletal muscle, heart muscle) (Figure 4b). Moreover, in Figure 4c it is possible to observe the known subcellular localization of the DE RNAs. A high number of variable organelles emerge, and it is possible to suggest that the cells of the SAT present with ubiquitary differential expression, in the nucleus as well as the cytoplasm, the mitochondria, and the cytoskeleton, indicating a profound alteration in all cellular functions. , and KEGG database (c), represented as dotplot and ClueGO graphs. For dotplots, the yaxis represents the name of the pathway, the x-axis represents the gene ratio, dot size represents the number of different genes, and the color indicates the adjusted p-value. For ClueGO graphs, each pie segment refers to the % of terms present per group and * p < 0.05.

Characteristics of DE RNAs in SAT of Men Versus Women: Interaction, Tissue Expression and Cellular Localization
The coding DE RNAs were comprehensively characterized for their interaction, tissue expression, and subcellular localization (Figure 4). To do so, the STRING database was used to construct an interaction network of DE RNAs, where the nodes are proteins

Analysis of Differential Expression in Relation to Diseases-Development: Implications for Sex Differences in Secondary Co-Morbidities Development
In order to investigate the potential secondary diseases that could arise in either men or women, the ClinVar database was used. This highlighted, as it is to be expected, diseases associated with the X or Y chromosome (spermatogenic failure, Y-linked or X inactivation, familial skewed), but, interestingly, also metabolic diseases such as diabetes mellitus type 2 and phosphoglycerate dehydrogenase deficiency (Figure 5a). This suggests that the genetic signature found here could be responsible for the development of secondary metabolic complications, and indeed, heart failure also emerges (ventricular tachycardia). The DisGENET database was used to highlight the specific genes' curated evidence, linking them to disease conditions, subdivided by disease category. Specifically, IRS1 and PHGDH emerge as the DE RNAs causative of nutritional and metabolic associated diseases, as they are extremely relevant and characterized in diabetes, lipidemia, obesity (IRS1), and phosphoglycerate dehydrogenase deficiency (PHGDH) (Figure 5b). IRS1 is also involved in the development of coronary heart diseases, whilst the TRDN gene is involved in other forms of heart failure, ventricular tachycardia, and Romano-Ward Syndrome (Figure 5c). The PHGDH is also responsible for the development of neurological genetic Neu-Laxova Syndrome (Figure 5d). Lastly, some genes (HMGN5, IRS1, PHGDH, MXRA5, GOLIM4, and SFRP2) are associated with the development of multiple cancer types (Figure 5e). Interestingly, all the discussed genes present a decreased expression in men for susceptibility genes, and it could thus be speculated that women patients could be more prone to developing associated pathologies, but of course each gene should be considered for its specific correlation with the analyzed condition. As gender medicine is becoming of crucial relevance, a focus was also given to the potential drugs selectively targeting DE RNAs (Figure 5f). Indeed, four specific genes that could be targeted emerged, and these are PHGDH, KDM5D, IRS1, and ZFY.

Role of Non-Coding RNAs in Gender Imbalace in Obesity: Interaction with DE RNAs and Correlated Functions
As lncRNAs are becoming more and more relevant in the context of sex-biology and obesity, a special focus was given to the non-coding DE RNAs, as these could bridge the void in the present knowledge concerning the role of lncRNAs in sex-specific obesity. Indeed, non-coding DE RNAs also present diversifies expression profiles, as shown in Figure 6a with Heatmap visualization and in Figure 6b with PCA clustering. The volcano plot reports the non-coding DE RNAs (Figure 6c), and indeed, 19 non-coding RNAs emerge as differentially expressed, and they are thus reported in Table 3. The possible lncRNAs' implication in diseases was analyzed using the LncRNADiseasev2.0 database [48] to investigate their implication in possible comorbidities development. AC10889.1, TTTY15, XIST, AC006157.1, PAX8-AS1, SNHG25, and JPX were mainly involved in cancer development. Interestingly, apart from AC10889.1 and AC006157.1, the other lncRNAs were implicated, not only in cancer, but also in neurodegenerative diseases (Alzheimer's disease and Huntington's disease) and schizophrenia, suggesting they could contribute to co-morbidity development.     A co-interaction network was built, highlighting all the ncRNAs that interact with coding genes (Figure 6d). The network highlights how 10 out of the 15 ncRNAs interact with 22 coding DE RNAs. This suggest that lncRNAs could exert their functions through a modulation of these interacting coding DE RNAs. In order to see in which processes these lncRNAs could be implicated, the functional enrichment of the DE RNAs specifically associated with lncRNAs was analyzed to identify the specific ontology and pathway implication. Indeed, significant terms were found in all categories analyzed, which include gene ontology (GO), with 49 significant pathways for molecular function (MF), 808 significant pathways for biological processes (BP), and 57 significant pathways for cellular component (CC) pathway analyses, which highlighted 1 significant pathways for KEGG, 53 significant pathways for Reactome, and 30 significant pathways for WikiPathways, along with 1209 regulatory motifs in DNA (TF), 34 genes targeted by miRNAs (MIRNA), and implications in protein databases such as the human protein atlas (HPA, 146), CORUM (3), and implications for human phenotypes (HP, 593) ( Figure 6e, Table S4).
We next specifically investigated these categories and found that the 3 CORUM terms concerned the Sam-68-p85 P13K-IRS-1-IR signaling complex, the cytoplasmic dynein-2 complex, and 18S U11/U12 snRNP (Figure 7a). The GO terms analysis in MF highlighted 49 significant pathways, and the results are displayed as pie graph in Figure 7b. ClueGO MF analysis highlights MAPK activity with 74.07% terms per group, along with catalytic activity with 11.11% terms per group and cation binding and RNA binding with 7.41% terms per group (Figure 2b). Interestingly, the highest perturbation seems to be due to the activity of MAP kinases, which are involved in the regulation of both normal and pathological adipogenesis [49]. For GO BP analysis, 808 significant terms emerged, and again the results are displayed as a pie graph in Figure 7c. ClueGO analysis highlights, as the top significant term, an implication for genes in chordae embryonic development (63.64% terms per group). However, numerous terms involve metabolic processes, which is interesting in the obesogenic context, including terms implicated in the glycoprotein biosynthetic process (2.12% terms per group) and the nucleobase-containing compound metabolic process (4.55% terms per group) (Figure 7c). Lastly, 57 significant terms emerged in the GO CC analysis, and the results are displayed as a pie graph in Figure 7d. Interestingly, a strong implication concerns the nucleus (27.27% terms per group), along with endocytic vesicles (24.24% terms per group). Another compartment that seems to be highly involved is the Golgi apparatus, with 48.48% terms per group (Figure 7d). esting in the obesogenic context, including terms implicated in the glycoprotein biosynthetic process (2.12% terms per group) and the nucleobase-containing compound metabolic process (4.55% terms per group) (Figure 7c). Lastly, 57 significant terms emerged in the GO CC analysis, and the results are displayed as a pie graph in Figure 7d. Interestingly, a strong implication concerns the nucleus (27.27% terms per group), along with endocytic vesicles (24.24% terms per group). Another compartment that seems to be highly involved is the Golgi apparatus, with 48.48% terms per group (Figure 7d).  Each pie segment refers to the % of terms present per group; * p < 0.05 and ** p < 0.01. Pathways analysis with Reactome (e) and Wikipathways (f) databases, represented as a dotplot, where the y-axis represents the name of the pathway, the x-axis represents the gene ratio, dot size represents the number of genes, and the color indicates the adjusted p-value.

Pathway Analysis of Coding and Non-Coding Transcripts Associations Highlights Specific Processes-Involvement
DE RNAs were subjected to pathway analysis in three well renowned databases: Reactome, WikiPathways, and KEGG. Reactome analysis highlighted 53 significant pathways, as reported in Table S4. The top 20 pathways, ranked by their significance, are reported in Figure 7e. Interestingly, Wnt regulation and glycosylation seem to be the most represented terms. Indeed, Wnt signaling is represented by "signaling by Wnt", "negative regulation of TCF-dependent signaling by Wnt ligand antagonists", and "Wnt5A-dependant internalization of FZD2, FZD5, and ROR2", whereas glycosylation and glycan biosynthesis are involved.
Thirty pathways emerged from the WikiPathways database, as reported in Table S4. The first 20 pathways, ranked by their significance, are reported in Figure 7e. The Wnt signaling pathway is represented multiple times, along with in insulin signaling (normal and diabetic condition) and type 2 diabetes. KEGG analysis with g:Profiler highlighted only 1 significant pathway, also previously mentioned, which is the Wnt signaling pathway, as reported in Table S4.

Correlation between DE RNAs and Anhtropometric Parameters
Given the importance of lncRNAs in the biology of sex and obesity, particular attention was given to the correlation between DE RNA and anthropometric parameters, to shed new light on the role of lncRNAs in male/female obese SAT. For each patient, different parameters were considered, and their correlation with DE RNAs was displayed as a scatter plot. The full list of genes correlated with each parameter is reported in Supplementary Table S5 (Figure 8c,d). None of the lncRNAs-DE RNAs were related to fasting glucose, BMI, or age (Figure 8e-g). No significant correlations in DE RNAs were found when comparing gene expression to Cholesterol, LDL, AST, ALT, and white blood cell count (Table S5).

Discussion
The molecular mechanisms underlying the influence of sex-related differences on the development of obesity are still unknown. Moreover, the role of non-coding RNAs in this context has never been investigated. Transcriptional characterization of specific tis-

Discussion
The molecular mechanisms underlying the influence of sex-related differences on the development of obesity are still unknown. Moreover, the role of non-coding RNAs in this context has never been investigated. Transcriptional characterization of specific tissues in obese patients is of crucial relevance in highlighting the role of sex differences in obesity, fully characterizing the underlying molecular signature that leads to a different disease profile in obese men or women. To this end, the results hereby presented are an in-depth analysis of transcriptional differences occurring in the SAT of a total of 10 subjects: 5 obese women (OBF, age 41 ± 12.5 years, BMI 38.2 ± 4.6 kg/m 2 ), and 5 obese men (OBM, OBM, age 42.4 ± 6.58 years, BMI 36.9 ± 3.5 kg/m 2 ). Although the total data size can be considered limited, this study is still relevant in providing novel potential targets to be validated in a wider cohort. Indeed, the results highlight a subset of 51 DE RNAs, 70% of which were coding genes and 30% of which were non-coding DE RNAs. We imposed minimum |log 2 FC| of 1 and an FDR lower than 0.1 as cut-offs to differentially expressed genes in order to have a stringent significance threshold.
Through KEGG, Reactome, WikiPathways, GO enrichment analysis, and CORUM analysis, we were able to deeply investigate pathways, molecular functions, biological processes, cellular components, and protein complexes to gain more in-depth insights over pathways in which coding and non-coding genes were involved. The GO analysis for molecular function highlighted terms such as "protein demethylase activity", "histone demethylase activity", and "demethylase activity" as a result of the increased expression of KDM5D and UTY, and this confirms the existence of differences in epigenetic regulation in males with respect to females [50]. Moreover, the ClueGO MF analysis also highlighted MAPK activity, along with cation binding, RNA binding, and catalytic activity. Furthermore, the GO analysis for biological process highlighted metabolic processes, interesting in the obesogenic context, and protein modifications as also supported by ClueGO analysis. Pathways related to epigenetic and transcriptional regulation as well as Wnt signaling pathways, highlighted by the three renowned databases Reactome, WikiPathways, and KEGG, gives a new perspective on the mechanisms through which sex differences are present in obesity development. Indeed, although the role of Wnt signaling remains largely unexplored, evidences highlights how Wnt signaling contributes to obesity-associated metabolic dysfunction by increasing adipose tissue inflammation [51]. In our work, the two differentially expressed genes involved in this pathway are ROR2 and SFRP2, both with decreased expression in males, suggesting a decreased activation of Wnt signaling mediated by these molecules.
Using the ClinVar database, we were able to evaluate the onset of potential secondary diseases in both men and women. Along with sex-associated diseases, we also highlighted the possible onset of metabolic diseases such as diabetes mellitus type 2. Moreover, the DisGENET database allowed us to highlight specific genes related to disease conditions, divided by category. Specifically, IRS1 and PHGDH emerge as the DE RNAs causative of nutritional and metabolic associated diseases, and these are two regulators extremely relevant and characterized in diabetes, lipidemia, obesity (IRS1), and phosphoglycerate dehydrogenase deficiency (PHGDH). Furthermore, genes presented with a decreased expression in men, suggesting that female patients could be more prone to developing associated pathologies.
Given the importance of lncRNAs in sex-biology and obesity, we also focused on non-coding DE RNAs, as these could bridge the void in present day knowledge concerning the role of lncRNAs in sex-specific obesity. Using the WGCNA R package we were able to evaluate the interactions ncRNAs and coding genes, highlighting how 10 out of the 19 ncRNAs interact with 22 coding DE RNAs. These lncRNAs could thus exert their function through a modulation of the coding genes, although functional in vitro analysis is needed to validate this mechanism. Indeed, lncRNAs can act through several distinct modes, influencing gene expression at multiple levels, starting from chromatin re-arrangements and transcriptional and translational modulation [23,24,52]. Biological characterization will shed light on which processes these lncRNAs affect, which could include adipogenesis, insulin resistance, inflammatory response, or other key mechanisms present in obesity. The lncRNAs hereby reported can also be associated with the insurgence of specific diseases, such as cancer development, neurodegenerative diseases (Alzheimer's disease and Huntington's disease), and schizophrenia, suggesting they could also be crucial players in co-morbidities development and thus affect obesity development and prognosis. In order to see in which processes these lncRNAs could be implicated, the functional enrichment of the involved DE RNAs was analyzed to identify the specific ontology and pathway implication. As for all DE RNAs, the functional enrichment analysis highlighted a dysregulation in MAPK activity pathways as well as RNA binding and Wnt signaling, suggesting a role for lncRNAs in altering expression of the coding genes implicated in obesity through these processes. LncRNAs have been recognized as viable biomarkers for a wide number of diseases and could be used in the design of lncRNA-based therapies for obese patients [21,53].
Moreover, correlations with anthropometrical parameters highlighted how the DE RNAs expression could be influenced by specific biological parameters that differ between men and women, providing a new prospective on possible gene modulation.

Conclusions
In conclusion, even if this work is to be considered a pilot study and will need to be validated in a wider cohort, it provides a comprehensive analysis of the implications of both coding and non-coding genes in relation to obesity and sex differences. Moreover, it represents an interesting basis for the development of new patient-specific therapeutic strategies.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/jpm11050352/s1, Table S1: Clinical features of enrolled patients, Table S2: List of Primers used for Real-Time PCR validation, Table S3: Enrichment analysis of DE RNAs as obtained with g:Profiler. Terms were ranked according to their Fold Change, and a Bonferroni-Hochberg FDR of 0.05 was used as the threshold. The full list of significant pathways (adjusted p.value < 0.05) is reported, Table S4: Enrichment analysis of DE RNAs correlated with nc DE RNAs as obtained with g:Profiler. Terms were ranked according to their Fold Change and a Bonferroni-Hochberg FDR of 0.05 was used as the threshold. The full list of significant pathways (adjusted p.value < 0.05) is reported, Table S5: Correlation analysis of anthropometrical parameters with DE RNAs. The full list of genes correlated with each parameter is reported, and a white cell indicates non-significant correlations, a blue cell terms with * p < 0.05, a green cell terms with ** p < 0.01, and an orange cell terms with *** p < 0.001.

Institutional Review Board Statement:
The present study is in accordance with the Declaration of Helsinki and it was approved by the Ethical Committee of IRCCS Istituto Auxologico Italiano (Ethical Committee approval code #2020_10_20_04).

Informed Consent Statement:
A signed informed consent was obtained from each enrolled patient for tissue sampling.

Data Availability Statement:
The raw data obtained from the RNA-seq analysis is deposited on the Gene Expression Omnibus repository with the accession number GSE166047.