Identication of Key Pathways and Genes in the Subcutaneous Fat Tissue from Different Sexes of Pigs using WGCNA Method

Background: Obesity is a complex disease that will endanger human life and health, and obesity is closely related to subcutaneous fat deposition. Some studies have shown the signicant association between subcutaneous fat deposition and sex. However, the molecular mechanisms for this association are still unclear. The objective of this study is to identify key pathways and genes related to the subcutaneous adipose tissue between different sexes of pigs. Results: A total of 16 coexpression modules were detected, of which two sex-related modules were found. Functional enrichment analyses showed that the genes in the purple module were mainly enriched in some pathways related to lipid metabolism such as Regulation of lipolysis in adipocytes. The genes in the magenta module mainly involved in some pathways related to immunity such as Defense response to virus and Immune response. Furthermore, 7 extracellular matrix (ECM)-related genes and 7 lipid metabolism-related genes were identied in the purple module. In the magenta module, 17 genes and 2 transcription factors participating in Type I interferon response were identied which were associated with immunity. Conclusions: The present study identied key pathways and genes in subcutaneous adipose tissue from different sexes of pigs, which provided some insights into the molecular mechanism involving in fat formation and immunoregulation between pigs with different sexes. These ndings may be helpful for breeding in pig industry and obesity treatment in medicine. matrix was transformed into a topology matrix, and TOM was used to measure the correlation of gene pairs. According to 1-TOM, average linkage hierarchical clustering was performed to classify genes with coherent expression proles into gene modules. Dynamic cutting algorithm was used to identify gene modules from the system cluster tree. Module eigengene (ME) was dened as the rst principal components, and was the representative in module genes. Module membership (MM) was dened as the correlation between ME and gene module. Gene signicance (GS) was indexed by log10 transformation of the P value of T-test. The only requirement was that GS of 0 indicates that the gene was not signicant with regard to the biological question of interest. The GS could take on positive or negative values. Module signicance (MS) was dened as the average GS for all the genes in the module. More detailed description of WGCNA was presented by an original article [13].

The expression matrix containing 5000 genes were used to reconstruct the gene coexpression network using WGCNA method. Pearson correlation matrix among genes were converted into a strengthened adjacency matrix by power β=5 based on scale-free topology criterion with R 2 =0.9 (Fig. 1a). Topological overlap measure (TOM) of each gene pair was calculated, and 16 gene coexpression modules were identi ed by average linkage hierarchical clustering according to TOM-based dissimilarity (1-TOM) (Fig. 1b). There were large differences in the number of genes between the modules, and the smallest midnightblue module contained 59 transcripts, while the largest turquoise module contained 914 transcripts. The purple module contained 227 transcripts, and the magenta module contained 254 transcripts (Table S1).
The eigengene adjacency heatmap depicting the interaction of the identi ed modules was showed in Fig. 1c, while the network heatmap plot of all the genes was showed in Fig. 1d. As shown in Fig. 1e, the correlation analysis between Module eigengene (ME) and sex showed that the purple module was the most signi cantly representative module (cor=-0.73 and P=5e-07), and the magenta module was the third signi cantly representative module (cor=0.5 and P=0.002). Gene signi cance (GS) analysis across modules (P=2.2e-259) showed that the purple module and magenta module had the highest and second highest GS, respectively (Fig. 1f). The module membership (MM) in the purple module possessed the most signi cant correlation across all modules (cor=0. 71 and P<4e-36, Fig. 2a), and the MM in the magenta module possessed the second signi cant correlation across all modules (cor=0.63 and P<1.7e-29, Fig. 2b).
Thus, the purple module and magenta module were selected for further analyses. The expression heatmap of all genes in the purple and magenta module were respectively showed in Fig. 2c and 2d.
GO and KEGG analysis and key genes identi cation GO and KEGG enrichment analysis were performed on all genes in the purple module and magenta module using the Database for Annotation, Visualization and Integrated Discovery (DAVID) database. In the purple module, GO enrichment results showed that 9 biological processes (BPs), 7 cellular components (CCs) and 4 molecular functions (MFs) were signi cantly enriched (P<0.05). KEGG enrichment analysis showed that 6 KEGG pathways were signi cantly enriched (P<0.05). Top 1 BP, CC and MF terms and Top 4 KEGG pathways were showed in Table 1. KEGG pathway analysis showed that the genes in the purple module were signi cant enriched in Regulation of lipolysis in adipocytes, Aldosterone synthesis and secretion, Ras signaling pathway and MAPK signaling pathway, and so on. GO analysis showed that the most signi cant enriched CC term was Extracellular space. The most signi cant enriched BP term was Muscle organ morphogenesis. The most representative MF term was Heparin binding. In addition, the functional analysis for the magenta module was performed, GO and KEGG enrichment results showed that 23 BPs, 3 CCs, 11 MFs and 12 KEGG pathways were signi cantly enriched (P<0.05). Top 2 BPs, Top 1 CC and MF terms and top 4 KEGG pathways were showed in Table 2. In KEGG analysis, genes in the magenta module mainly enriched in pathways related to immunity disease, including In uenza A, Hepatitis C, Measles and Herpes simplex infection, and so on. The most signi cant enriched CC term was External side of plasma membrane, and the signi cant enriched BP terms were related to immunity process, including Defense response to virus and Immune response. The most representative MF term was Double-stranded RNA binding.
In the purple module, the selection criteria of key genes were that genes participated in two terms, and at least involved in one KEGG pathway. Besides, all the genes involved in the pathway of regulation of lipolysis in adipocytes were selected as key genes. So, 10 key genes, LIPE, NPR1, AQP7, NPY1R, PLA2G16, LOC100522721, FGF10, FGFR2, CACNA1G and IGF1 were identi ed and showed in Fig. 3a. In the magenta module, the selection criterion of key genes was that genes involved in at least three terms. 12 genes, FAS, TNFSF10, DDX58, MX1, IFIT1, ADAR, OAS2, OAS1, CXCL9, CXCL10, STAT1 and IRF9 were selected as key genes (Fig. 3b).

PPI network construction and hub genes identi cation
All genes in the purple and magenta modules were used to perform PPI network analysis. At combined score > 0.4, a PPI network including 116 nodes and 169 edges was constructed for the purple module (called the rst PPI network), and was showed in Fig. 4a. For the magenta module, a PPI network containing 153 nodes and 814 edges was established and showed in Fig. 4b (called the second PPI network). Top genes with higher scores were identi ed by centrality analyses in the two modules using the centiScape. The result of centrality analysis for genes with top 10 centrality scores in the rst PPI network and top 20 centrality scores in the second PPI network were respectively showed in Table 3 and Table 4. In the purple module, top 10 genes were COL1A2, LOC100738989,  MMP2, COL1A1, SPARC, POSTN, SERPINH1, IGF1, FMOD and FKBP10, separately. In the magenta module, top 20 genes were USP18, MX1, STAT1, ISG15, IRG6,  MX2, IFIT3, IFIT2, IFIT1, DDX58, OAS1, PARP9, DDX60, IFI44L, LFI35, IFI44, PARP14, IRF9, XAF1 and OAS2, separately. To further excavate the hub gene in the two modules, the cytoHubba was used to screen out hub genes in the whole PPI network. According the Maximal Clique Centrality (MCC) score, top 10 genes were selected as hub genes in the rst PPI network. A sub-network including 10 nodes and 36 edges in the rst PPI network was identi ed and showed in Fig. 4c. These hub genes were COL1A2, MMP2, COL1A1, SPARC, POSTN, SERPINH1, FMOD, ENSSSCG00000017581, ENSSSCG00000029074 and LOC100738989. According to the MCC score, top 20 genes were selected as hub genes in the second PPI network. A sub-network containing 20 nodes and 187 edges in the second PPI network was identi ed and showed in Fig. 4d. These hub genes were USP18, MX1, STAT1, ISG15, IRG6,  ENSSSCG00000027918, MX2, IFIT3, IFIT2, IFIT1, DDX58, OAS1, PARP9, DDX60, IFI35, IFI44, IFI44L, IRF9, XAF1 and OAS2. Ultimately, 8 hub genes, COL1A2, LOC100738989, MMP2, COL1A1, SPARC, POSTN, SERPINH1 and FMOD in the purple module were identi ed simultaneously using centiScape and cytoHubba. Function enrichment analysis showed that the 8 genes mainly involved in some pathways including AGE-RAGE signaling pathway in diabetic complications, Relaxin signaling pathway and Proteoglycans in cancer (Fig. 5d). The signi cant enriched BP terms were Extracellular matrix (ECM) organization and Extracellular structure organization, and so on (Fig. 5a). The signi cant enriched MF terms were Collagen binding and Extracellular matrix binding, and so on (Fig. 5b). The signi cant enriched CC terms were Extracellular matrix and Collagen-containing extracellular matrix, and so on (Fig. 5c).
In the magenta module, 19 hub genes, USP18, MX1, STAT1, ISG15, IRG6, MX2, IFIT3, IFIT2, IFIT1, DDX58, OAS1, PARP9, DDX60, IFI35, IFI44, IFI44L, IRF9, XAF1 and OAS2 were identi ed simultaneously using centiScape and cytoHubba. The functional enrichment analysis was performed for these hub genes. KEGG pathways had Hepatitis C, Coronavirus disease-COVID-19, Measles, In uenza A, Epstein-Barr virus infection, Herpes simplex virus 1 infection, NOD-like receptors signal pathway, RIG-I-like receptors signal pathway and Human papillomavirus infection and so on (Figure 6d). The signi cant enriched BP terms were Defense response to virus, Response to virus, Immune effector process, Defense response to other organism and Response to type I interferon (Fig.  6a).The enriched MF terms were Double-stranded RNA binding and GTP binding, and so on (Figure 6b). The enriched CC terms were Dendritic spine and Neuron spine, and so on (Fig. 6c).

Differential expression analysis
By the functional analysis and PPI network analysis, 18 genes and 24 genes were selected as key genes in the purple and magenta module, respectively. In the purple module, differential expression analysis showed that 17 genes except for NPY1R (P=0.0502) had signi cant differences in the purple module between two sex groups (P<0.05 or 0.01, Fig. 7a and 7b). Among 17 genes with signi cant differences in expression, FMOD, SPARC, POSTN, SERPINH1 and LOC100522721 were signi cantly different (0.01< P<0.05), and other 12 genes were remarkable signi cant different (P<0.01). Expression level of NPR1, LIPE, AQP7, PLA2G16 and LOC100522721 in female pigs were signi cantly higher than that in male pigs, while that of COL1A2, MMP2, COL1A1, SPARC, POSTN, SERPINH1, FMOD, LOC100738989, IGF1, FGF10, FGFR2 and CACNA1G in females were signi cantly lower than that in males. In the magenta module, 19 genes except for USP18, MX1, IRG6, DDX60 and ADAR including DDX58 (P<0. and TNFSF10 (P<0.01) were signi cantly differentially expressed between two sex groups, and all the genes were more highly expressed in female pigs ( Fig. 7c and 7d).

Key pathways
In our study, two modules including the purple and magenta modules closely related to sex in the subcutaneous fat tissue of pigs were identi ed using WGCNA method. Function enrichment analysis showed that the genes in the purple module were mainly enriched in some pathways including Regulation of lipolysis in adipocytes, Aldosterone synthesis and secretion, Ras signaling pathway and MAPK signaling pathway. These pathways are related to lipid metabolism [27,28]. Yao et al. found that male gonadal steroids de ciency caused by castration signi cantly raised the serum lipid levels and increased fat accumulation in the pigs, and DEGs identi ed between the boars and barrows were mainly involved in metabolism of lipid, carbohydrate, amino acid and immune diseases pathways [29]. A study showed that male gonosteroids in boars could signi cantly reduce blood lipid levels and fat accumulation in pigs. Moreover, functional enrichment analysis showed that the downregulated genes in entire males as compared with immunocastrated males were signi cantly enriched in Adipocytokine signaling pathway, AMPK signaling pathway, Insulin signaling pathway and PPAR signaling pathway [30]. However, The genes in the magenta module were mainly involved in some pathways related to disease immunity, such as In uenza A, Hepatitis C, Measles, Herpes simplex infection, Defense response to virus and Response to virus, and so on. A study demonstrated that sex steroids also could regulate the immune system [31]. After hepatitis C virus and in uenza A virus infection, in ammatory immune responses and repair of damaged tissue differed between the sex and impacted the outcome of infection [32,33]. Therefore, it could be inferred that sex differences might affect fat deposition and immunity in subcutaneous fat tissue of pigs by affecting relevant pathways related to lipogenesis, lipid metabolism and disease immune.

Hub genes in the purple module
In the purple module, 8 hub genes including COL1A2, MMP2, COL1A1, SPARC, POSTN, SERPINH1, FMOD and LOC100738989 were identi ed by the PPI analysis. Functional enrichment analysis showed that these genes were signi cantly enriched in AGE-RAGE signaling pathway in diabetic complications, Relaxin signaling pathway, Extracellular matrix (ECM) organization, Extracellular structure organization, ECM, Collagen-containing ECM, Collagen binding and ECM binding. The result was very similar to that from the study of Poklukar et al. [30], and their ndings showed that the upregulated genes in entire males as compared with immunocastrated males and surgical castrates were signi cantly enriched in extracellular region/matrix (ECM) cellular components, ECM receptor interaction and focal adhesion pathways. Some genes responsible for the differences in backfat deposition between the three male sex categories were identi ed including COL1A2, COL6A3, POSTN, P4HA3, DCN, FMOD, MMP2 and MMP27 [30]. In our study, COL1A2 gene was identi ed to have the highest degree in the PPI network, which indicated that the gene was at the core of the network. In the ECM remodeling, COL1A2 and COL1A1 genes involves in the synthesis of collagen, which are the major components of ECM [34]. POSTN gene is crucial for the collagen cross-linking and the ECM maintenance [35,36]. Similarly, FMOD gene is required for proper collagen folding and ECM stabilization [37]. SPARC gene can promote the collagen brils development and posttranslational modi cations [38] and stabilizes the tissue and maintains the balance of lipogenesis and lipolysis [39]. SERPINH1 gene is responsible for the correct folding and secretion of collagen [40], and can regulate the ECM proteins expression [41]. MMP2 gene involves in the ECM degradation [42]. In the ECM remodeling, COL1A2, COL1A1, FMOD, POSTN, SPARC, SERPINH1 and MMP2 play a key role, which in uence adipose tissue remodeling. In our study, the 7 genes were all signi cant differences in expression between different sexes, and expression of the 7 genes in boars were higher compared with sows. Jeong et al. measured the expression levels of ECM-related genes in different adipose tissues from bulls, cows and steers of Korean cattle (Hanwoo), and found that the expressions of ECM-related genes in the omental adipose tissue of cows and steers are decreased, and expression levels of most ECM-related genes were generally similar between cows and steers [43]. Poklukar et al. found that castration of male pigs resulted in the downregulation of genes involved in ECM dynamics [30]. The results of these studies are similar to this study. Thus, it might be inferred that sex differences could affect the expression levels of ECMrelated genes (COL1A2, COL1A1, SPARC, POSTN, FMOD, SERPINH1 and MMP2), progress to in uence adipose tissue remodeling, and the higher expression of these genes in boars might make boars have a greater fat remodeling capacity than sows.
In the purple module, 10 hub genes including LIPE, NPR1, NPY1R, AQP7, PLA2G16, LOC100522721, IGF1, FGF10, FGFR2 and CACNA1G were identi ed by the function enrichment analysis. Except for NPY1R, other genes were signi cant differences in expression between different sexes. Among these DEGs, the expression levels of LIPE, NPR1, AQP7, PLA2G16 and LOC100522721 in sows were signi cantly higher than that in boars, and other four genes including LIPE, NPR1, AQP7 and PLA2G16 excepting for LOC100522721 participated in Regulation of lipolysis in adipocytes. LIPE gene encodes hormone sensitive lipase which cleaves free fatty acids from the intracellular triacylglycerol [44], and LIPE gene is signi cantly associated with intramuscular fat content and fatty acid composition of pigs [45,46]. NPR1 gene promotes the proliferation, differentiation and lipolysis of preadipocytes [47]. AQP7 promotes glycerine excretion in adipose tissue and are highly expressed in adipose tissue [48], and AQP7 de ciency is associated with adipose tissue triglyceride accumulation and adult obesity [49]. Similarly, PLA2G16 is a major phospholipase A2 of adipose tissues that is speci c to adipocytes [50]. PLA2G16 encodes a major negative regulator of adipocyte lipolysis, and PLA2G16 de cient mice has a markedly higher rate of lipolysis, and a signi cantly reduced adipose tissue mass [51,52]. According to the above, LIPE, NPR1, AQP7 and PLA2G16 mainly involve in the regulation of lipolysis. At present, the four genes had not been reported in the study of fat deposition in pigs between different sexes.
Furthermore, among these DEGs, the expression levels of IGF1, FGF10, FGFR2 and CACNA1G were signi cantly lower in sows than in boars. Functional enrichment analysis showed that IGF1 participated in Ras signaling pathway, and FGF10 and FGFR2 participated in Ras signaling pathway and MAPK signaling pathway, while CACNA1G participated in Aldosterone synthesis and secretion and MAPK signaling pathway. A study showed that IGF1 action is inhibited in castrated animals, which affects adipocyte proliferation and differentiation [30]. FGF10, belongs to the broblast growth factor family, which is widely involved in the regulation of cell growth, proliferation, differentiation and regulation of metabolism [53,54]. FGFR2 appears to be the most commonly distributed growth factor receptors and mediates most biological functions of the FGF family. Some studies suggested that FGF10 stimulated preadipocyte proliferation and differentiation through activating FGFR2 [55,56]. Additionally, CACNA1G gene is reported to be associated with RFI, DMI and FCR [57,58], not been reported in adipocyte proliferation and differentiation of pigs at present. As the above, IGF1, FGF10 and FGFR2 promote adipocyte proliferation and differentiation. Currently, the three genes had not been reported in the adipocyte proliferation and differentiation of pigs between different sexes. Because the expressions of the three genes were higher in boars, which indicated that the three genes might more promote adipocyte proliferation and differentiation of boars. In summary, LIPE, NPR1, NPR1, PLA2G16, IGF1, FGF10 and FGFR2 in the purple module were lipid metabolism-related genes, which might play a different regulatory role in the process of fat deposition and adipocyte proliferation and differentiation of pigs with different sexes.

Hub genes in the magenta module
In the magenta module, 19 genes were identi ed as hub genes by the PPI analysis. Functional enrichment analysis showed that these genes were mainly enriched in some diseases and signal pathways, such as Hepatitis C, Coronavirus disease-COVID-19, RIG-I-like receptors (RLRs) signal pathway and NOD-like receptors (NLRs) signal pathway, and so on. RLRs and NLRs can sense viral components, such as double/single-strand RNA and DNA. The RLRs play essential roles in the production of type I interferons (IFNs) and proin ammatory cytokines in cell type-speci c manners [59]. Among these genes, 15 genes, DDX58, STAT1, IRF9, OAS2, OAS1, IFIT1, IFIT2, IFIT3, ISG15, IFI35, IFI44, IFI44L, MX2, XAF1 and PARP9 were signi cant differences in subcutaneous adipose tissue of pigs between two sex groups. It has been reported that DDX58 gene belongs to one of the crucial members of RLRs family, and DDX58 can produce type I interferon (IFN) [60, 61], which activates STAT1 and STAT2 and recruits IRF9 gene. The STAT1-STAT2-IRF9 complex activates IFN-stimulated gene factor 3 (ISGF3), which drives the expression of IFN-stimulated genes (ISGs) [ [70], IFI44L [70], PARP9 [70], XAF1 [71] and MX2 [72] belonged to the Type I ISGs. OAS1 and OAS2 are the 2′-5′ oligoadenylate synthetase family which consists of antiviral enzymes induced by IFN and is responsible for destroying virus-derived dsRNA [65]. IFIT1, IFIT2 and IFIT3 belong to interferon induced proteins with tetratricopeptide repeats family, which have a great importance in antiviral and immune regulation [67] and restrict virus infection through altering cellular protein synthesis [73]. ISG15, IFI35, IFI44 and IFI44L belong to IFN-induced protein family. ISG15 gene is one of the most abundantly induced ISGs and can effectively control certain viral and bacterial infections [74]. A study indicated that IFI35, IFI44 and IFI44L played important roles in interferon signaling, anti-bacterial and immune regulation [75]. Beside, PARP9 gene promotes the expression of pro-in ammatory genes [76]. XAF1 gene involves in immunoregulatory [71]. MX2 gene mediates antiviral innate immunity [77].
In the magenta module, 12 genes were identi ed as key genes by the function enrichment analysis. Among these genes, 10 genes, DDX58, STAT1, IRF9, OAS2, OAS1, IFIT1, CXCL9, CXCL10, TNFSF10/TRAIL and FAS were signi cant differences between two sex groups. The six genes, DDX58, STAT1, IRF9, OAS2, OAS1 and IFIT1 have been discussed before. Some studies showed that CXCL9 [78], CXCL10 [79] and TNFSF10/TRAIL [70] were also Type I ISGs. The proin ammatory chemokines CXCL9 and CXCL10 associate with in ammation and proliferation [80, 81]. TRAIL-system plays a role in regulating autoimmunity [82]. FAS gene involves in activation of the Type I IFN/STAT1 axis [83]. From the above, a total of 19 DEGs in the magenta module participated in Type I IFN response, and it could be inferred that DDX58 gene rst promoted the production of type I IFN, while FAS gene participated in activating the Type I IFN/STAT1 axis. Second, the combination of STAT1, STAT2 with IRF9, which drove the expression of 15 ISGs (OAS2, OAS1, IFIT1, IFIT2, IFIT3, ISG15, IFI35,  IFI44, IFI44L, PARP9, XAF1, MX2, CXCL9, CXCL10 and TNFSF10) which encoded corresponding immune proteins to carry on speci c immune action. Currently, the 17 genes and two TFs (STAT1 and IRF9) had not been reported in the disease immunity of pigs between different sexes. In this study, the expressions of these genes were higher in female pigs. Therefore, the 17 genes, DDX58, FAS, OAS2, OAS1, IFIT1, IFIT2, IFIT3, ISG15, IFI35, IFI44, IFI44L, PARP9, XAF1, MX2, CXCL9, CXCL10 and TNFSF10) and two TFs (STAT1 and IRF9) participated in Type I IFN response which might play a different regulatory role in immunity of pigs between different sexes and might enhance the immunity ability of female pigs though Type I IFN stimulated pathway in cells.
However, there were some limitations in this study. On one hand, small sample size limited the statistical power to identify the hub genes. On the other hand, molecular biological experiments were required to validate the function of these hub genes.

Conclusions
In summary, 16 coexpression modules were identi ed using WGCNA method, and two of these modules were closely related to sex. Functional enrichment analysis showed that the purple module was related to lipid metabolism, while the magenta module was associated with immunity. 7 ECM-related genes (COL1A2, COL1A1, SPARC, POSTN, FMOD, SERPINH1 and MMP2) and 7 genes related to lipid metabolism (LIPE, AQP7, NPR1, PLA2G16, IGF1, FGF10 and FGFR2) were identi ed in the purple module. In the magenta module, 17 genes and two TFs participating in Type I IFN response were identi ed which were associated with immunity. In conclusion, the identi ed key pathways and genes in subcutaneous adipose tissue from different sexes of pigs, which provided some insights into the molecular mechanism involving in fat formation and immunoregulation between male and female pigs. These ndings may be helpful for breeding in pig industry and obesity treatment in medicine.

RNA-seq data and samples information collection
The normalized gene expression dataset of GSE61271 was used in this study, and was retrieved from the public NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo/). The GSE61271 dataset was from a study of Kogelman et al. [84], consisting of 36 pig (crossbred F2 of Duroc and Göttingen miniature pig) samples with 100 kg weight (17 females and 19 males). WGCNA WGCNA was used to construct the gene coexpression network, and identify the coexpression gene modules. The WGCNA package (version 1.13) based on R was used to perform WGCNA [13]. First, the expression matrix was converted into an adjacency matrix, and an unsupervised coexpression relationship was constructed based on the adjacency matrix using Pearson correlation coe cients for gene pairs. The correlation adjacency matrix was strengthened by power β (soft threshold), and the power parameter was selected based on scale-free topology criterion.
Second, the adjacency matrix was transformed into a topology matrix, and TOM was used to measure the correlation of gene pairs. According to 1-TOM, average linkage hierarchical clustering was performed to classify genes with coherent expression pro les into gene modules. Dynamic cutting algorithm was used to identify gene modules from the system cluster tree. Module eigengene (ME) was de ned as the rst principal components, and was the representative in module genes. Module membership (MM) was de ned as the correlation between ME and gene module. Gene signi cance (GS) was indexed by log10 transformation of the P value of T-test. The only requirement was that GS of 0 indicates that the gene was not signi cant with regard to the biological question of interest. The GS could take on positive or negative values. Module signi cance (MS) was de ned as the average GS for all the genes in the module. More detailed description of WGCNA was presented by an original article [13].
Finally, the statistical signi cance of the relationship between modules and sexes was analyzed by calculating the Pearson correlation coe cient. For studying the genes in the module related to sex, the top two modules with the high GS scores among all modules were selected as the modules of signi cant for further analysis.

PPI network construction and analysis
The interactive relationships among genes encoding proteins in the signi cant gene coexpression module were analyzed by constructing a PPI network. The interactive information among genes encoding proteins was retrieved from online the Search Tool for the Retrieval of Interacting Genes (STRING) database (version 11.5, https://string-db.org/). The gene pairs with a combined score ≥ the medium con dence 0.4 were used to construct the PPI network. Cytoscape GO and KEGG analysis GO and KEGG pathway terms of all genes in the signi cant module were analyzed using the online DAVID database (version 6.80, https://david.ncifcrf.gov/).
The cut-off criterion was set as P-value<0.05. Cytoscape (v3.8.0) software was used to construct and visualize the interactive relationships between genes and functional enrichment terms in the whole network. GO and KEGG enrichment analyses for hub genes in PPI network were implemented using the R-package clusterPro ler software [86]. The cut-off criterion was set as q-value<0.5. GO annotation result includes three main bodies: biological process (BP), molecular function (MF) and cellular component (CC).

Hub genes identi cation
Hub genes were identi ed by the cytoHubba algorithm and centiScape 2.2 plugin in the cytoscape software, and key genes were identi ed using GO and KEGG network analysis. CytoHubba algorithm was used to gain top genes in the whole PPI network, and the criterion of selecting hub genes was that top nodes ranked by Maximal Clique Centrality (MCC) [87]. The centiScape algorithm was applied to calculate the centrality of each gene in the whole PPI network, including Degree Centrality, Betweenness Centrality, Closeness Centrality, EigenVector Centrality, Stress Centrality and Radiality Centrality, and so on [88]. The genes with higher Centrality scores were identi ed as hub genes.

Statistical analysis
All statistical analyses were performed using R software. The comparison of gene expression in subcutaneous fat tissues with females and male pigs was analyzed using Mean±Standard Deviation. T-test was used to estimate the statistical difference, and 0.01<P<0.05 was considered as signi cant difference, and P<0.01 was considered as remarkable signi cant different, and P>0.05 was considered as not signi cant difference in expression level between two groups.   Tables   Table 1 The results of GO and KEGG enrichment analysis for the purple module using DAVID tool Table 2 The results of GO and KEGG enrichment analysis for the magenta module using DAVID tool Table 3 Top 10 genes with higher scores identi ed by centrality analyses using centiScape in the purple module  Table 4 Top 20 genes with higher scores identi ed by centrality analyses using centiScape in the magenta module      (gray lines) between nodes indicated the interaction of genes in the network. Size of circles represented degree number of node, and the bigger the size, the more degree the node. Color represented betweenness, and the darker the color, the more betweenness the node. (c) PPI sub-network for the purple module.

Abbreviations
There were 10 nodes and 36 edges in the rst sub-network. (d) PPI sub-network for the magenta module. There were 20 nodes and 187 edges in the second sub-network. Size of circles represented degree number of node, and the bigger the size, the more degree the node. Color represented rank, and the darker the color, the more rank the node.  indicated the degree of enrichment, with smaller P.adjust indicating terms that were more likely to play signi cant functional roles. Top 10 terms ordered by P.adjust for the KEGG enrichment analysis.

Figure 7
The expression boxplot of key genes in the purple and magenta modules. (a) 8 key genes, and (b) 10 key genes in the purple module. (c) 12 key genes, and (d) 12 key genes in the magenta module. 0.01<P<0.05 was considered as signi cant difference, and P<0.01 was considered as remarkable signi cant different, and P>0.01 was considered as not signi cant difference between two groups in expression level. X-axis represented the genes, and Y-axis represented value of gene expression. F (red) represented female pigs, m (blue) represented male pigs.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. TableS1Thenumberofgenesineachofthe16modules.docx