Transcriptome Analysis of Muscle Reveals Potential Candidate Genes and Pathways Affecting Intramuscular Fat Content in Duroc Pigs

Background Intramuscular fat (IMF) content plays a key role in meat quality of pork. Duroc is a common commercial pig breed, which is extensively used as the terminal male parent of DLY (Duroc × Landrace × Yorkshire) commercial pigs. For identification potential candidate genes and pathways, IMF content of 28 purebred Duroc pigs were measured, and transcriptome differences of muscles from individuals with divergent IMF content were analyzed based on RNA sequencing. Results As a result, a total of 309 DEGs were detected using edgeR and DESeq2 (p < 0.05, |log2(fold change)| ≥ 1) between the the IMF-High and -Low groups. Functional analysis of the DEGs revealed 27 Gene Ontology (GO) terms, two pathways (q < 0.05) and 23 key DEGs, which closely related to lipid metabolism and fat cell differentiation. Among these key DEGs, LEP, SFRP1, CIDEC, PNPLA3 and MT3, are potential candidate genes affecting IMF content due to their biological functions, the significant correlation between these genes’ expression with IMF content values in the 28 Duroc samples and proofs from previous studies. Meanwhile, a weighted gene co-expression network analysis (WGCNA) also revealed seven out of these key DEGs, including ADIPOQ, FABP4, PLIN1, CIDEA, ACAA1, CYP2A19 and ADIRF, were hub genes in turquoise module which may exhibit important functions in influencing IMF content. Furthermore, two pathways, regulation of lipolysis in adipocytes and PPAR signaling pathway, were potential pathways affecting IMF content because of their significant enrichment for the DEGs and co-expression genes of turquoise modules. Conclusions Our study demonstrates that some DEGs and pathways detected have essential role in regulating IMF content, and are supposed to be potential candidate genes and pathways affecting the trait in Duroc pigs. The study provides information for understanding molecular mechanism of IMF content, and would be of help for improving pork quality of Duroc pigs.

3 Background Meat quality is a main economic trait in pig production, which can be evaluated by multiple indicators, such as intramuscular fat (IMF) content, pH, water holding capacity, color and tenderness. Of these, IMF content is arguably one of the most important one, and has closely correlations with some other meat quality traits, including flavor, juiciness and tenderness [1]. IMF corresponds to the amount of fat accumulates between muscle fibers, chemically, which covers the sum of phospholipids, triglycerides and cholesterol.
High level of IMF content or marbling has a positive influence on the eating quality of pork [2]. However, for the past century, thinner backfat is considered as an important parameter for pig production. Although, this selection led to higher muscularity and growth, IMF content, meat juiciness and tenderness obviously decreased [3]. Meanwhile, IMF content is very difficult to be improved by traditional breeding methods. Hence, to satisfy the consumer preference for IMF content, it is worthwhile to understand the molecular basis of IMF content for improving it based on molecular breeding, IMF content is a complex trait, which is regulated by many genes. Quantitative trait loci (QTL) and genome-wide association studies (GWAS) have been performed to determine the underlying genetic factors influencing it in previous studies. To date, 652 QTLs for IMF content have been deposited in pigQTLdb (https://www.animalgenome.org/cgibin/QTLdb/SS/traitmap?trait_ID=16). However, most of QTLs have been distinguished mainly by low-density microsatellite markers using linkage mapping method, and span a large chromosomal region, which always comprise hundreds of genes. So, it is difficult for the further fine-mapping to identify causative genes [4]. The emergence of highthroughput genotyping platforms, SNP arrays, makes it possible find genetic variants and QTLs associated with the traits of interest in a narrower region. Despite these progresses, compared with other traits, the progress of studying the genetic basis of IMF content trait 4 is still slow.
Next generation sequencing of RNA (RNA-seq), a method based on next-generation sequencing (NGS), has the potential to give unprecedented details about the transcriptional landscape of a certain organism. So, it is a new method for studying complex quantitative traits controlled by many interacting genes. Compared with real-time PCR and microarrays, RNA-seq, with low background noise and high dynamic ranges of gene expression levels quantification, integrates advanced molecular biology techniques and bioinformatics [5]. To date, RNA-seq has been carried out in different livestock species, such as cattle [6], sheep [7], and pig [8], to identify genes involved in IMF content and other meat quality traits. Several studies on transcriptome profiling of porcine longissimus dorsi (LD) between breeds with different IMF content [9], or between individuals with extremely high or low IMF content in Berkshire [10] have been reported.
However, studies on transcriptome profiling in Duroc for identify differentially expressed genes (DEGs) and pathways related to IMF content are still not reported.
Duroc pig population is extensively used as the terminal male parent of the DLY (Duroc × Landrace × Yorkshire) commercial pigs due to its high meat quality and large muscle mass. As a coloured breed, Duroc has substantially higher levels of IMF content than white breeds [11]. Moreover, in the study, we found the IMF content of Duroc still varied significantly (1.17% to 4.23%). Accordingly, individuals with extreme IMF content in Duroc population are good materials for transcriptomics study of IMF content. Hence, in the present study, we performed RNA sequencing of LD muscles from Duroc pigs with variant IMF content, analyzed the transcriptome differences between the groups with extremely high and low IMF content, determined key DEGs and biological pathways affecting IMF content. This study provides valuable information for understanding the molecular basis underlying IMF content of pigs.

Results
Phenotypic variation between extreme groups IMF content of LD muscles from a purebred population of 28 Duroc pigs were measured.
The average of IMF content was 2.36% with a range of 1.17% to 4.23%. To select suitable samples used for the transcriptome analysis, we conducted Principal Component Analysis (PCA) using 300 genes with the most variance. The cumulative proportion of the first two principle components (PCs) was 38.19% (PC 1 = 29.90%, PC 2 = 8.29%, Fig.1A). According to the result of PCA, D8, D15, D19, and D21 individuals were assigned to the IMF-High group, and D13, D17, D18, and D26 individuals were assigned to the IMF-Low group. It was consistent with those selected based on IMF content values, suggesting close relationship existed between expressions of these genes and IMF content values. The average of IMF content was 3.70% and 1.47% for IMF-High and -Low group respectively. And the difference between the two groups was of extreme significance (p < 0.01, Table S1).

Summary of sequencing data
In this study, the transcriptome of LD muscle of the 28 Duroc pigs were sequenced using paired-end RNA-seq approach. Totally, 52.43 -86.93 million raw reads per sample were generated. After filtering 2.31% of raw reads, an average of 62.56 million clean reads were used for following analyses (Table S1). of these, 94.36% of clean reads were successfully mapped to Sus scrofa11.1 genome assembly, with 90.34% being uniquely mapped (Table S1). Additionally, out of the reads mapped to the reference genome, 91.58% and 3.86% were located in exon and intron regions, respectively, and the remaining were in intergenic regions (Table S1).
According to alignment result of RNA sequencing data with Sus scrofa genome, a total of 26,918 genes were expressed in eight IMF content extremely individuals. The correlation coefficient of mean gene expression level between IMF-High and -Low groups was 0.96.
The high correlation showed that most of the genes had a similar expression behavior in the two groups. Taking into account those genes with fragments per kilobase of transcript per million fragments mapped (FPKM) value > 1 in more than four pigs, 18,411 and 18,950 genes were found in IMF-High and -Low group, respectively.

Identification of differential expression genes
Using edgeR, 457 genes were identified as DEGs when we compared the genes expression in terms of |log 2 (fold change)| ≥ 1 and p-value < 0.05 between two groups. Of these, 334 and 123 genes were up-and down-regulated, respectively. Whereas, using DESeq2 under the same screening criteria, 316 DEGs were identified. Of these, 241 and 75 genes were up-and down-regulated, respectively. By comparing the DEGs detected by edgeR and DESeq2, a total of 309 overlapped DEGs (Table S2) were discovered by both methods, and were used for the following functional enrichment analyses. Of these, 240 genes had a higher expression and 69 genes had a lower expression in the IMF-High group. The heatmap and hierarchical clustering of the 309 DEGs indicated that the expression patterns were consistent within groups and different between groups (Fig.1B, 1C).
Meanwhile, hierarchical clustering also exhibited the relationships of expression patterns between the eight IMF extreme individuals and other 20 individuals in the 28 Duroc pigs (Fig. 1C).
Furthermore, to confirm the DEGs, correlation analysis was also carried out between FPKM values of DEGs and IMF contents of 28 pigs measured in the study (Table S2). Of these 309 DEGs, 182 DEGs had medium correlation (|r| ≥ 0.3), and 107 were significantly correlated (p < 0.05) between expression and IMF contents, suggesting these DEGs have similar expression pattern in the 28 samples as the eight IMF extreme individuals.
Using DEseq2 and edgeR, we also screened DEGs under stricter conditions (|log 2 (fold change)| ≥ 1 and q-value < 0.05). As a result, 10 and 12 DEGs were identified respectively, and seven genes were overlapped (Table 1). Further gene annotation revealed that two of the seven DEGs (LEP and CIDEC ) directly related with lipid metabolism, three (SPP1, ENSSSCG00000034371 and SFRP1) located in the region of QTLs influencing IMF content and marbling. Furthermore, correlation analysis indicated that the expressions of LEP, CIDEC and SFRP1 were significantly correlated with IMF content (r = 0.55, p= 2.23E-03; r = 0.47, p = 1.20E-02; r = 0.45, p = 1.69E-02) in the 28 Duroc samples. Of these seven DEGs, two have not be annotated withe specific gene name, and the remaining five ones were regarded as key genes affecting IMF content and used for following protein-protein interaction (PPI) network analysis.

Gene Ontology (GO) enrichment analysis of differentially expressed genes
To have a functional view of our DEGs detected between the two groups, we carried out GO term enrichment analysis for them and determined 58 significantly enriched GO terms (q < 0.05, Fig.2A and Table S3). Interestingly, all the significant enriched GO terms were enriched for up-regulated genes. Among the 58 enriched GO terms, most of them belonged to biological process (BP) category, only two terms belonged to molecular function (MF) and cellular component (CC) category, respectively. Importantly, 27 of 58 GO terms were closely associated with lipid metabolism and fat cell differentiation, such as lipid catabolic process, low-density lipoprotein receptor particle metabolic process, white fat cell differentiation, and positive regulation of cholesterol efflux (Fig.2B).
Moreover, other significant GO terms related to some physiological and biological events, such as response to cytokine, response to tumor necrosis factor, alcohol metabolic, and detoxification, were also enriched.
8 Figure 2B illustrates the 18 up-regulated genes involved in the significantly enriched GO terms related to lipid metabolism and fat cell differentiation. Among these genes, APOE, ADIPOQ, PPARG, CEBPA, LIPE, MT3, and PNPLA3, were involved in more than five GO terms, including catabolism of triglyceride-rich lipoprotein constituents (APOE), control of fat metabolism and insulin sensitivity (ADIPOQ), regulation of adipocyte differentiation (PPARG), inducing adipocyte terminal differentiation (CEBPA), and triacylglycerol hydrolysis (LIPE and PNPLA3). Furthermore, the other enriched genes, such as FABP4 [12,13] and FASN [14] are candidate genes affecting IMF content by QTL mapping. By

Pathway enrichment analysis of differentially expressed genes
In addition, we conducted Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for the 309 DEGs. Five pathways were significantly enriched (q < 0.05, Table 2) for the up-regulated DEGs, while no one was detected for the down-regulated genes. Two of the significant pathways, regulation of lipolysis in adipocytes (q = 5.05E-03) and PPAR signaling pathway (q = 5.05E-03), were IMF-related. And these two IMF-related pathways covered eight genes, PLIN1, LIPE, PTGER3, ADRB1, FABP4, ACAA1, PPARG, and ADIPOQ, which were relevant to lipid metabolism and fat cell differentiation. It is noticeable that, except for PTGER3, other genes were all found in enriched GO terms for the DEGs.
Together with the five key genes detected under stricter conditions (|log 2 (fold change)| ≥ 1 and q-value < 0.05), we determined in total number of 23 key genes playing the vital role in regulating IMF content due to one overlapped (Table S4). Among these genes, some 9 (APOE, ADIPOQ, PPARG, CEBPA, LIPE, MT3, and PNPLA3) were involved in more than five GO terms, some (PLIN1, LIPE, ADRB1, FABP4, ACAA1, PPARG, and ADIPOQ) were simultaneously detected by GO and KEGG analyses. The results implied that there could be possible interactions among these genes. To further analyze the interaction relationships of the 23 key genes, we constructed PPI network analysis using Search Tool for the Retrieval of Interacting Genes (STRING).Except FAPB4, CYP2A19 and ADIPOQ, the other genes were successful matched to sus scrofa database of STRING. By the visualization of Cytoscape software, a PPI network consisted of 12 core proteins were constructed (Fig. 3). The top hub genes interacted with more than six of other genes were PPARG, LEP, LIPE and PLIN1.

Co-expressed genes associated with IMF content
The Weighted Gene Co-expression Network Analysis (WGCNA) relies on the assumption that strongly correlated expression levels of gene sets indicate those genes work cooperatively in related pathways and contribute to the resulting phenotype. By WGCNA, we constructed co-expression network with 5175 selected genes and identified 29 modules (Fig. 4A). Six modules, including skyblue module, orange module, turquoise module, pink module, darkgreen module, and darkred, were correlated with IMF content (correlation > 0.4).
For each module, we performed functional enrichment analysis, and found that genes in turquoise module remarkably participated in IMF-related biological processes or pathways.
For the co-expressed genes in turquoise module, 82 significant enriched GO terms were identified (q < 0.05, Table S5) in GO analysis. Among these, positive regulation of lipid localization and lipid particle were two terms related to lipid metabolism (Fig. 4B).
Meanwhile, two pathways were also identified in KEGG analysis (q < 0.05), regulation of lipolysis in adipocytes and PPAR signaling pathway, which were canonical pathways regulating fat deposition (Fig. 4B, Table S6). Notably, the two terms and pathways related to lipid metabolism or regulating fat deposition were all found to be significantly enriched ones in the functional analysis of the DEGs.
Above mentioned results revealed that genes in the turquoise module may exhibit important functions in influencing IMF content. To identify hub genes, we calculated the weight referring to connection strength between two genes in the turquoise module. The weight values ranged from 0.02 to 0.39. Correlation network was constructed with the weight value > 0.36 by Cytoscape (Fig. 4C). The network contained 39 genes, seven of which, ADIPOQ, FABP4, PLIN1, CIDEA, ACAA1, CYP2A19 , and ADIRF, were also among the key genes identified by GO and pathway enrichment analysis of DEGs. Furthermore, except CYP2A19, the expressions of the other six hub genes were all significantly correlated with IMF content (ADIRF, r = 0.42, p = 2.43E-02), or their correlation coefficients were above 0.3, which also proved that these hub genes were potential candidate genes regulating IMF content in Duroc.

Discussion
Duroc, one of the most excellent breeds for its high meat production and meat quality, is usually used as the terminal male parent of the commercial cross-bred pigs such as DLY pigs. Though Duroc is known for its high IMF content [11] compared to the other lean breeds, the IMF content of Duroc population still vary in a wide range. In this study, IMF content of 28 Duroc pigs, which were reared under the same environmental conditions and with the same nutrition levels, varied from 1.17% to 4.23%. Combined with phenotypic data and the result of PCA, four high IMF content individuals and four low IMF content individuals with significant difference (p < 0.0016) were selected from the 28 Duroc pigs.
High throughput RNA-seq technology was employed to identify DEGs between the two groups, and several functional enrichment analyses were followed to identify plausible 11 candidate genes and pathways affecting IMF content. The outcome of the study might be helpful for explaining genetic variations related to IMF content in Duroc.
DESeq2 and edgeR are both very popular and easy-to-use packages for differential expression analysis of RNA-Seq data [15]. In this study, these two methods were used for detecting DEGs, and overlapped genes detected by them were consider as DEGs, which would be helpful to reduce false-positive DEGs detected. As a consequence, 309 DEG were detected (Table S2)  Among the key DEGs detected in the study, LEP, SFRP1, CIDEC, PNPLA3 and MT3, were not only repeatedly discovered as the key genes in the functional enrichment analysis ( Table   1, Fig. 2B), but their expression were significant correlated with IMF content values in the 28 Duroc samples (Table S2). LEP gene encodes a protein, leptin, which is secreted by white adipocytes into the circulation and regulates the body's fat storage system through the inhibition of lipogenesis and the stimulation of lipolysis [16]. However, a study on LEP gene expression in IMF tissue between growing Erhualian and large white pigs, the expression of LEP may contribute to fast sediment of IMF [17]. SFRP1 is a Wnt antagonist in Wnt/β-catenin signaling pathway which plays an important role in regulating adipogenesis and controlling adipocyte differentiation [18]. CIDEC was discovered highly expressed in adipose tissues, and its expression level was significantly higher in obese pigs than in their lean counterparts [19]. The protein encoded by CIDEC, an adipocyte lipid drop protein, negatively regulates lipolysis and promotes triglyceride accumulation [20].
PNPLA3 gene, belongs to the patatin-like phospholipase (PNPLA) family, exhibits lipase and transacylase activity [21]. In human, an association was revealed between PNPLA3 expression and obesity as well as related phenotypes [22]. In pigs, PNPLA3 is mainly expressed in adipose and liver tissue, participating in the restoration of lipid homeostasis upon aberrant intracellular lipid accumulation. Consistent with our results, in Sodhi's study, PNPLA3 was also distinguished as one of important candidate genes affecting IMF content by comparative transcriptome analysis in different pig breeds [23]. MT3 protein, metallothionein 3, is a member of the mammalian metallothionein (MT) family. Wang et al.
found that MT3 was expressed in human adipocytes and played a critical role in protecting adipocytes from hypoxic damage [24]. By the model of MT3 knockout mice, Byun et al.
found MT3 may be involved in central leptin signaling and the consequent increase in peripheral energy expenditure [25]. However, the biological function of porcine MT3 needs to be further clarified. In summary, the above evidences suggested that they are potential candidate genes affecting IMF content. Except the five genes mentioned above, the other key genes, such as FABP4, CIDEA, PPARG, and ADIPOQ, whose correlation coefficients were all above 0.35 (Table S4), should also be paid attentions too.
To identify important pathways affecting IMF content, we conducted KEGG analysis using DEGs (p < 0.05 and |log 2 (fold change)| ≥ 1, Table 2) and co-expression genes in turquoise modules (Table S6). Regulation of lipolysis in adipocytes and PPAR signaling pathway were both identified as enriched pathways for the two gene sets. According to the description of KEGG database (https://www.genome.jp/dbget-bin/www_bget?map04923), regulation of 13 lipolysis in adipocytes is a pathway regulating triacylglycerol to release fatty acids and glycerol for other organs as energy substrates. Some studies have demonstrated that the IMF content in the LD is determined by a balance between fat accumulation and degradation [26]. It is implied that more lipid storage may lead to the increase of lipolysis [27]. Consequently, in the study, regulation of lipolysis in adipocytes pathway enriched in the High-IMF group compared with Low-IMF group. PPAR signaling pathway is a canonical pathway involving in lipid metabolism. The expression pattern of many candidate genes in this pathway were correlated with meat quality traits in LD muscle of pigs [28]. Previous studies showed that the PPAR signaling pathway was significant enriched pathway for DEGs of subcutaneous and intramuscular stromal vascular cells during adipogenic differentiation [29], as well as DEGs of transcriptome analysis between Laiwu and Yorkshire pigs' LD muscles [30]. Therefore, consistent with the previous studies, our results suggested that PPAR signaling pathway was a pivotal pathway affecting IMF content.

Conclusions
In this study, we performed a high throughput RNA sequencing to evaluate the transcriptome profiles differences of eight Duroc pigs with extreme IMF content. A total of 309 DEGs were screened by DESeq2 and edgeR. GO and pathway enrichment analysis, correlation analysis between gene expressions and IMF content in 28 Duoc pigs and WGCNA revealed some potential candidate genes, such as LEP, SFRP1, CIDEC, PNPLA3 and MT3, and pathways, such as regulation of lipolysis in adipocytes and PPAR signaling pathway. These potential candidate genes and pathways play an essential role in affecting IMF content of pigs, and deserved to carry further studies to unravel their specific mechanism on IMF content. This study provides useful information for understanding 14 molecular mechanism of intramuscular fat content in Duroc pigs.

Animals and sample collection
The pig specimens were from a Duroc breeding farm in Shandong province of China. A total of 28 purebred Duroc pigs at 108.29 ± 6.00 kg (mean ± SD) of body weight were used in this study to select pigs representing high-(n=4) and low-IMF (n=4) content groups. These pigs were fed ad libitum with commercial feed and water under the same environmental, and were slaughtered in one batch after euthanized by electric shock. All For the determination of IMF content, the LD muscle samples taken from the last 4 th and 5 th thoracic vertebra, which were removed adipose and connective tissues, were grounded after oven-drying. Then, IMF content were examined using the Soxhlet petroleum-ether method, and expressed as the weight percentage of wet muscle tissue. Each sample was measured in triplicate to ensure its accuracy.

RNA isolation, library construction and sequencing
Total RNA was isolated from 28 porcine LD samples using TRNzol reagent (Invitrogen, Life Technologies). RNA was examined with a NanoDrop 1000 instrument (Thermo Fisher Scientific) and Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA), and RNA integrity was confirmed using a 2100 Bioanalyzer (Agilent Technologies, 15 CA, USA). RNA libraries were constructed using NEBNext® Ultra TM RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations. Briefly, the mRNA was isolated from total RNA using poly-dT oligonucleotides attached to magnetic beads. The mRNA was fragmented, and then reverse-transcribed into first-strand cDNA with random hexamers. After generating double-stranded cDNA, end repair, A-tailing, adaptor ligation and cDNA purification and enrichment were then performed. Paired-end raw sequences (150 bp) were generated using an Illumina HiSeq4000 Sequencing System. Quality assessment of sequencing data and mapping of RNA-seq reads In order to ensure the quality of bioinformatics analysis, raw reads were firstly filtered to obtain clean reads. We filtered raw reads in FASTQ format by removing reads according to the following criteria: containing adapter sequences, ambiguous base-content > 0.1%, 50% of bases whose Q phred quality scores ≤ 20. Clean reads were then mapped to Sus scrofa reference genome (Sscrofa11.1) and the annotation database Ensembl Genes 95 using HISAT2 software [31]. FeatureCounts tool of subread package was used for calculating the number of reads mapped to each gene for estimating gene expression levels [32]. Gene expression levels were also measured in the expected number of fragments per kilobaseof exon per million fragments mapped (FPKM).

Differential gene expression analysis
For analyzing the expression correlation of 28 Duroc pigs and selecting IMF extreme individuals, PCA was performed with R package prcomp using 300 genes with the most variance in 28 Duroc pigs. Combined with phenotype data, eight pigs with extremely high and low IMF content were ultimately determined for the following transcriptome analysis.
Differential expression analysis was carried out by the R package DESeq2 [33] and edgeR [34]. Genes with |log 2 (fold change)| ≥ 1 and p-value < 0.05an adjusted p-value < 0.05 were regarded as DEGs for each methods. To reduce false-positive result, the overlapped DEGs detected by the two methods were regarded as the final DEGs and used for the following functional enrichment analyses. Meanwhile, the resulting p-values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate (FDR). And to detect key genes affecting IMF, we also conducted the differential expression analysis using stricter criteria, |log 2 (fold change)| ≥ 1 and q-value < 0.05. Of particular note is that, in the study, we did not use genders as a sample selection condition. Among the eight individuals, there were only two male pigs which belonged to IMF-High group. Therefore, the DEGs located in Y-chromosome were not taken into account.
To further confirm the DEGs, Person correlation analysis between the expression levels of DEGs and IMF content values of the data set of 28 Duroc pigs were conducted using the R package Hmisc. In addtion, to visualize the DEGs between the two IMF groups and among the 28 Duroc pigs, heatmap was generated using the R package pheatmap based on the FPKM values and hierarchical clustering was analyzed by R function hclust using correlation coefficients of DEGs obtained in Person correlation analysis as the distance between a pair of individuals.

Functional enrichment and PPI network analysis of DEGs
The DEGs detected by DESeq2 and edgeR were subjected to GO and KEGG pathway enrichment analysis using R package clusterProfiler [35]. GO terms and pathways with qvalue < 0.05 were considered to be significantly enriched. Genes, involved in the significantly enriched GO terms and pathways related to lipid metabolism and fat cell differentiation, were treated as key genes affecting IMF.
To evaluate the interactive relationships of the genes and identify hub genes, the key genes, identified by differential gene expression analysis under stricter criteria (q-value < 0.05) and GO and KEGG pathway enrichment analysis of DEGs, were uploaded into STRING (https://string-db.org/cgi/input.pl?sessionId=WcgfuVcCfC0u&input_page_show_search=on).
Active interaction sources of STRING included textmining, experiments, database, coexpression, neighborhood, gene fusion, co-occurrence with a medium confidence score

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.

Supplementary Files
This is a list of supplementary files associated with the primary manuscript. Click to download.