FOSL2 Is Involved in the Regulation of Glycogen Content in Chicken Breast Muscle Tissue

The glycogen content in muscle of livestock and poultry animals affects the homeostasis of their body, growth performance, and meat quality after slaughter. FOS-like 2, AP-1 transcription factor subunit (FOSL2) was identified as a candidate gene related to muscle glycogen (MG) content in chicken in our previous study, but the role of FOSL2 in the regulation of MG content remains to be elucidated. Differential gene expression analysis and weighted gene coexpression network analysis (WGCNA) were performed on differentially expressed genes (DEGs) in breast muscle tissues from the high-MG-content (HMG) group and low-MG-content (LMG) group of Jingxing yellow chickens. Analysis of the 1,171 DEGs (LMG vs. HMG) identified, besides FOSL2, some additional genes related to MG metabolism pathway, namely PRKAG3, CEBPB, FOXO1, AMPK, and PIK3CB. Additionally, WGCNA revealed that FOSL2, CEBPB, MAP3K14, SLC2A14, PPP2CA, SLC38A2, PPP2R5E, and other genes related to the classical glycogen metabolism in the same coexpressed module are associated with MG content. Also, besides finding that FOSL2 expression is negatively correlated with MG content, a possible interaction between FOSL2 and CEBPB was predicted using the STRING (Search Tool for the Retrieval of Interacting Genes) database. Furthermore, we investigated the effects of lentiviral overexpression of FOSL2 on the regulation of the glycogen content in vitro, and the result indicated that FOSL2 decreases the glycogen content in DF1 cells. Collectively, our results confirm that FOSL2 has a key role in the regulation of the MG content in chicken. This finding is helpful to understand the mechanism of MG metabolism regulation in chicken and provides a new perspective for the production of high-quality broiler and the development of a comprehensive nutritional control strategy.

Introduction

Muscle glycogen (MG), as the main source of glucose for muscle glycolysis (Godfrey and Quinlivan, 2016), can supply adenosine triphosphate to muscles through its own degradation, which is closely related to muscle growth and development (Puolanne and Immonen, 2014). Under special physiological conditions (such as intense exercise), it can maintain the stability of blood sugar and regulate the balance of glucose metabolism. In addition, MG content affects various meat quality indicators, such as meat color, tenderness, and pH of livestock and poultry meat (Henckel, 1996; Berri et al., 2007; Le Bihan-Duval et al., 2008). The degradation of MG to lactic acid changes the pH, which affects meat quality; thus, pH is an important indicator of meat quality (Alnahhas et al., 2015).

A previous genome-wide association study (GWAS) analysis of the MG phenotype conducted using the whole-genome resequencing data of 474 chickens identified FOS-like 2, AP-1 transcription factor subunit (FOSL2) transcription factor subunit as a candidate gene whose product affects MG content (Liu et al., 2020). FOSL2 belongs to the activator protein 1 (AP1) transcription factor family, which is expressed in various tissues and is involved in fat metabolism, bone development, and the pathogenesis of various diseases, such as cancer (Foletta et al., 1994; Chinenov and Kerppola, 2001; Eferl and Wagner, 2003). A study on the involvement of FOSL2 in type 2 diabetes found that the epigenetic downregulation of FOSL2 mRNA and protein expression levels by DNA hypermethylation may induce the occurrence of type 2 diabetes in patients (Li et al., 2016). FOSL2 plays a key role not only in human physiology and disease (Eferl et al., 2008), but also in the differentiation of chondrocytes by controlling osteoclast-specific transcription (Karreth et al., 2004). Additionally, FOSL2 regulates bone formation by osteoblasts in mouse and human mainly by regulating the expression of collagen isoforms collagen, type I, alpha 2 chain (COL1A2) and bone gamma-carboxyglutamate protein (BGLAP) (Bozec et al., 2010). Previous studies have reported that the FOSL2 gene protein product activates CEBPB transcription in prostaglandin E₂ (PGE₂)–activated osteoblasts. Although the regulatory effect of transcription factor CEBPB on chicken MG content has also been reported (Sibut et al., 2011), the role of FOSL2 expression in tissues other that bone tissue, such as its role in MG metabolism, has not been described to date. An interesting experiment conducted in mice found that osteoblasts specifically express FOSL2 to regulate Adipoq and Bglap levels, which in turn influences systemic glucose and insulin metabolism (Bozec et al., 2013). Although the study did not clearly reveal that FOSL2 affects the content of MG, it suggested that FOSL2 is closely related to glucose metabolism. MG is synthesized from monomers of uridine diphosphate glucose and serves as a storage form of energy in the body. Thus, MG and glucose are inseparably related.

In this study, breast muscle tissues from the high-MG-content (HMG) group and the low-MG-content (LMG) group of Jingxing yellow (JX-Y) chickens were subjected to transcriptome sequencing. Combined with the weighted gene coexpression network analysis (WGCNA), the protein products of the identified candidate genes were further analyzed to study the effect of gene transcription on chicken MG content from the perspective of gene expression and its relationship with phenotype. Considering the complex function of FOSL2 in chickens of different breeds and sex, we studied the effects of the FOSL2 expression level in breast muscle of Beijing You (BJY) chicken. Furthermore, we also investigated the effects of FOSL2 expression on the MG content in chicken fibroblast DF1 cells. This study found that FOSL2 reduces MG content likely by downregulating the activity of CEBPB.

Materials and Methods

Ethics Statement and Animals

All animals and experimental protocols used in this study were approved by the Beijing Institute of Animal Science, Chinese Academy of Agricultural Sciences (the scientific research department responsible for animal welfare issues) (no.: IAS2019-21).

In this study, the experimental samples were derived from 474 female 98-day-old JX-Y chickens, a synthetic breed by the Institute of Animal Science of Chinese Academy of Agricultural Sciences, with high fat content and unique flavor, which is a better material for meat quality research. Additionally, 30 male and 30 female 98-day old BJY chickens, a unique Chinese local breed, with high flavor characteristics, were selected to evaluate the influence of the breed and sex.

All experimental birds were raised in an environmentally controlled room, in three-story step cages. Basal diets were formulated based on the National Resource Council (1994) requirements and the Feeding Standards of Chickens established by the Ministry of Agriculture, Beijing, China (2004).x1 All birds were provided ad libitum access to feed and water up to 98 days of age. All chickens were individually euthanized using carbon dioxide anesthesia and exsanguination by severing the carotid artery at 98 days of age after 12-h fasting (no additional feed was supplied, and the feed trough was not emptied). After slaughtering, the skin was cut open, and internal tissue samples were collected along the direction of muscle fibers. The pectoralis major muscle samples were weighed, snap-frozen in liquid nitrogen, and stored at −80°C for determination of MG, RNA sequencing, and fluorescence quantitative real-time polymerase chain reaction (qRT-PCR).

Measurement of MG Content

Muscle Glycogen Assay Kit from Nanjing JianCheng Bioengineering Institute (Nanjing, China) was selected. Tissue samples of approximately 100 mg were first cut, with a fluctuation of no more than 0.5 mg among samples. Then, the breast muscle tissue samples were hydrolyzed with 300 μL of lye in a boiling water bath for 20 min and cooled afterward with running water. Subsequently, 0.1 mL of a 5% MG detection solution and 0.9 mL distilled water were added to the hydrolyzed samples and boiled for 5 min. Then, 2 mL of color reagent was added, and the optical density (OD) value was measured at a wavelength of 620 nm using a spectrophotometer. Finally, the MG content of the tested sample was calculated.

Samples Collection for RNA Sequencing and qRT-PCR Analysis

The MG content of the JX-Y chicken female population ranged from 0.2 to 9.0 mg/g (Supplementary Material 1). We included individuals with an MG content of less than 1 mg/g in the LMG group and individuals with MG content of more than 4 mg/g in the HMG group. Transcriptome sequencing was performed on samples from 16 chickens (including three individuals in the HMG group, three individuals in the LMG group, and 10 individuals with random phenotypes). A total of six chickens in the HMG and LMG groups were analyzed by differential gene expression analysis to identify the differentially expressed genes (DEGs) between the two groups. In order to exclude the differences in FOSL2 gene expression caused by extreme phenotypes, we included 10 individuals with random phenotypes on the basis of six individuals for WGCNA. The MG content of BJY chicken was measured in 30 male and 30 female individuals. The MG content in male BJY chickens ranged from 0.8 to 5.1 mg/g, and in female BJY chickens, it ranged from 0.8 to 4.2 mg/g (Supplementary Material 1). We considered male and female individuals with MG content of less than 1 mg/g as the LMG group, and individuals with MG content of more than 2 mg/g as the HMG group. Fluorescence qRT-PCR analysis was performed on six male individuals and six female individuals (three individuals in the HMG group and three individuals in the LMG group).

RNA Extraction and Sequencing

Total RNA was extracted from the breast muscle tissue samples using TRIzol reagent (Invitrogen, Carlsbad, CA, United States). The quality of the RNA was assessed after separation by electrophoresis on a 1% agarose gel, and RNA concentration was determined by spectrophotometry using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, United States). The OD 260/280 values of all samples were within the range of 1.8 to 2.0 for RNA sequencing and qRT-PCR analysis.

We used a cDNA library construction method previously described by Chen et al. (2019). The mRNA was enriched by binding of the mRNA poly-A tail to oligo (dT)–coated magnetic beads and fragmented into small pieces. Single- and double-stranded cDNAs were synthesized using mRNA as a template. The double-stranded cDNA was purified using the QIAQuick PCR purification kit (QIAGEN, Valencia, CA, United States). After purification, end repair, and ligation to sequencing adapters, agarose gel electrophoresis was used for fragment separation and fragment size selection. Finally, PCR enrichment was performed to obtain the final cDNA library. RNA sequencing was performed on an Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, United States) by Frasergen (Wuhan, China), and 150-bp paired-end reads were generated (Supplementary Material 2).

After quality control, the HISAT2 software (Kim et al., 2015) was used to align the second-generation sequence of each sample to the reference genome sequence galGal6.0. The number of reads compared to the transcript of the sample was counted by the comparison result of bowtie2 called by RSEM (Li and Dewey, 2011) and converted by the FPKM (fragments per kilobase of transcript per million bases) method (Trapnell et al., 2010).

For biological duplication, we used the DEseq2 software to analyze the significance of differences in gene and transcript expression (Love et al., 2014) and used the edgeR software for differential expression analysis without biological duplication. The standard for screening the total DEGs is | log2FC| > 1.

qRT-PCR Analysis

All PCR primers were designed at or just outside exon/exon junctions to avoid the amplification of residual genomic DNA using the Primer-BLAST on the National Center for Biotechnology Information (Bethesda, MD, United States) website, and specificity was determined using BLASTN (Supplementary Material 3). First-strand cDNA was synthesized using 2 μg of total RNA, random primers, and oligo(dT) primers, with the Transcriptor First Strand cDNA synthesis Kit (Takara, Dalian, China) following the manufacturer instructions. Target mRNAs were quantified by qRT-PCR analysis using SYBR Green Master Mix (Takara). The qRT-PCR reaction of each sample was performed on the QuantStudio 7 Flex system (Applied Biosystems, Shanghai, China) using 40 cycles (95°C for 3 min, 95°C for 3 s, and 60°C for 34 s). The amplification reaction for each sample was performed in triplicate. The collected data were analyzed using the 2^–ΔΔCT method (Livak and Schmittgen, 2001), and all the results were normalized to the β-actin rRNA gene (Livak and Schmittgen, 2001; Liu et al., 2019).

WGCNA

Network analysis was performed using WGCNA R package version 1.66 (Langfelder and Horvath, 2008). The WGCNA methods have been successfully applied to gene expression data from microarrays (Kadarmideen et al., 2011) and RNA sequencing platforms in animal studies (Kogelman et al., 2014). Before performing WGCNA analysis, the genes obtained through transcriptome sequencing were filtered to remove genes whose average expression level was less than 1. Next, we used the WGCNA program package to construct a weighted gene coexpression network. We first calculated the Pearson correlation between paired genes and constructed a Pearson correlation matrix based on all paired genes. Subsequently, the optimal β value (the β value of JX-Y chicken population = 14) that makes the gene distribution conform to the scale-free network was selected, and the weighted adjacency matrix was constructed. Then, based on the correlation expression value, the adjacency relationship matrix was converted into a topological overlap matrix, and each topological overlap matrix was used as a hierarchical clustering analysis (Zhang and Horvath, 2005). Finally, based on the topological overlap matrix, we used the dissimilarity between genes to cluster the genes and used the dynamic shearing method to cut the tree into different gene modules (i.e., gene clusters with high topological overlap).

After constructing the modules, we calculated the module eigengene (ME) (Zhao et al., 2010), which is defined as the first principal component of the standardized expression profiles. Module–trait relationships (Kogelman et al., 2014), assessed by Pearson correlation between the ME and phenotypic values, were used to select potential modules related to phenotypes. In addition, to assess the relationship between genes and traits, Pearson correlation between gene expression and phenotypic values was calculated and used as gene significance (GS) (Liu et al., 2016), and the average of the absolute GS values within each module was used to determine the module GS (MS).

Comprehensive Analysis of PPI Network

We used the Search Tool for the Retrieval of Interacting Genes (STRING) database¹ to assess protein–protein interaction (PPI) data (Szklarczyk et al., 2015). In addition, in order to determine the relationship between FOSL2 and CEBPB, we used the STRING database and converted the results visually by using the Cytoscape software.

DF1 Cell Culture, Transfection, and Glycogen Detection

The chicken fibroblast DF1 cells were obtained from the cell bank of the Chinese Academy of Sciences. The DF1 cells were cultured in Dulbecco modified Eagle medium supplemented with 10% fetal bovine serum (Gibco, Grand Island, NY, United States) and 1% penicillin-streptomycin (Gibco). The cell lines were cultured at 37°C in a humidified incubator with 5% CO₂ (Li et al., 2020).

The FOSL2 (Gene ID: 421416) overexpressing lentivirus and empty vectors were constructed by Ubigene (Guangzhou, China) and denoted as YOE-LV001-FOSL2 and YOE-LV001-Ctrl, respectively. Cells were plated in six-well plates at a density of 1 × 10⁶ cells per well, and YOE-LV001-FOSL2 and YOE-LV001-Ctrl were transfected into DF1 cells by the Polybrene, a transfection-promoting reagent provided by Ubigene. The formula for virus dosage is V (μL) = 1,000 × MOI × N/T, where MOI is multiplicity of infection, N is number of cells, and T is lentiviral infectious titer.

The transfected cells were selected with 2 μg/mL puromycin for 2 days. Using an inverted microscope, the imaging system is Olympus BX41, and the images are taken at a 100× field of view. Then cells were collected and tested for glycogen content. The amount of glycogen in DF1 cells was quantified using a glycogen content assay kit (Solarbio, BC0340, China) according to the manufacturer’s instructions (Zhao et al., 2017). The cell experiment was repeated three times with samples in triplicate.

Statistical Analyses

The significance of the differences between groups was tested by the Student t test using the SPSS software version 22.0 (IBM Corp., Armonk, NY, United States). Confidence limits were set at 95% and p < 0.05 (^∗) or p < 0.01 (^∗∗) was considered significant. Data are presented as the mean ± standard error of the mean (SEM).

The GLM (general linear model) procedure in the SAS 9.4 software was used to analyze the differential expression of target genes of the HMG and LMG groups in BJY chickens, and the model was as follows:

$y_{ij}=u+P_j+G_i+e_{ij}$

where y_ij is gene expression quantity, u is population mean gene expression, P_j is sex effect, G_i is MG group effect, and e_ij is random error.

Results

Comparison of MG Content in the HMG and LMG Chicken Groups

To study MG metabolism in breast muscle tissue from HMG and LMG groups of chickens, the MG content in breast muscle tissue samples was measured. The results revealed significant differences in the MG content between chickens from the HMG and LMG groups, as shown in Figure 1A. The MG content in the HMG group chickens was significantly (p < 0.01) higher than that in the LMG group chickens.

FIGURE 1

Figure 1. The phenotypic analysis and the RNA sequencing analysis results. (A) The MG content in chickens from the HMG and LMG groups (P = 0.0003). (B) Hierarchical clustering analysis. The hierarchical clustering analysis was performed based on DEGs; the heatmaps of all six samples revealed that the gene expression profiles in the same group were very similar. (C) Changes in gene expression profiles of Jingxing yellow (JX-Y) chicken from the HMG and LMG groups. Red dots represent upregulated genes, blue dots represent downregulated genes, and black dots (No) represent insignificantly DEGs. (D) Expression levels (fold change) of representative genes involved in MG metabolism pathways according to transcriptome analysis data of chickens from the HMG and LMG.

Differential Expression Analysis of Gene Expression Data From HMG and LMG Chicken Samples

The hierarchical clustering analysis showed that only individuals within the same group clustered more closely (Figure 1B). The comparison of the results of the transcriptome sequencing analysis of individuals from the LMG group and HMG group identified a total of 1,171 DEGs, of which 376 were upregulated and 795 were downregulated (Figure 1C and Supplementary Material 4).

The expression of classic glucose metabolism genes in the HMG and LMG groups is shown in Figure 1D. In the HMG group, AMPK and PRKAG expression is higher than that in the LMG group (p < 0.01). The expression of other genes, such as FOSL2 and CEBPB, was higher in the LMG group than that in the HMG group (p < 0.01).

Construction of Weighted Gene Coexpression Network and Module Identification

A total of 15,934 genes were used to construct a weighted gene coexpression network. The gene modules shown in Supplementary Figure 1A were obtained based on the difference in hierarchical clustering. Based on the height value that quantifies the coexpressed similarity between modules, the modules with a value of less than 0.35 are merged. After merging, a total of 16 modules were identified. These coexpression modules are all represented in different colors in Supplementary Figure 1B, and the modules are named by the colors.

FOSL2 Modules Related to MG Traits

We searched for the gene modules most related to MG traits by correlating 16 ME and MG traits (Figure 2A). The ME midnight blue module shows the most significant negative correlation with MG traits in breast muscle tissue (r = −0.68; P = 0.004), which suggests that the genes in this module may play a critical role in MG metabolism.

FIGURE 2

Figure 2. (A) Relationship between gene modules and MG traits in chicken breast muscle tissue. Each row in the table corresponds to a module, and each column to a trait. Each cell contains the corresponding correlation and P value. The module name is shown on the left side of the heatmap. The correlation is shown on the right side of the heatmap, in accordance with the color legend (red means positive correlation, and green means negative correlation). (B) Scatterplots of gene significance versus module membership for MG modules, along with the indicated correlation coefficients and P values. There is a highly significant correlation between GS and MM in the midnight blue module (cor. = 0.82, P = 5.9e-56). (C) Protein–protein interaction network. Each node represents a kind of protein, and the line between the nodes indicates that the pair of proteins has an interaction relationship.

In the midnight blue module, 621 genes that are highly correlated with the corresponding modules and traits were clustered together (Figure 2B). Besides the candidate genes FOSL2 and CEBPB, some genes related to the classical sugar metabolism pathway, such as MAP3K14, SLC2A14, SLC12A4, PPP2CA, SLC38A2, PPP2R5E, and IGBP1, were also found in the midnight blue module. This indicates that the protein products of the FOSL2 and CEBPB genes and those of the sugar metabolism pathway-related genes may coparticipate in similar biological pathways and coregulate MG metabolism. The data in Supplementary Material 5 show the correlation between the aforementioned classic glucose metabolism–related genes and MG traits. The results reveal that the expression of the FOSL2 and CEBPB genes is negatively correlated with MG content. The network analysis of all the genes of the midnight blue module in the JX-Y chicken WGCNA using the STRING software to predict the structure of the protein network revealed that there is an interaction between FOSL2 and CEBPB (Figure 2C).

In vivo Validation of FOSL2 Gene Expression in Different Populations

In order to verify whether FOSL2 and CEBPB can play the same role in the regulation of the MG content in different breeds and different sex groups, BJY chicken from the HMG and LMG groups were used for quantification of their expression.

The HMG and LMG groups included six chickens (three males and three females). The transcript abundance of FOSL2 and CEBPB was verified by qRT-PCR analysis. The data in Table 1 show that the transcript abundance of the FOSL2 and CEBPB was significantly upregulated in the LMG group (p < 0.01). These results were consistent with the results of the transcriptome analysis of JX-Y chicken muscle tissue, indicating that FOSL2 and CEBPB also play a role in reducing MG content in BJY chicken.

TABLE 1

Table 1. The relative quantification (RQ) of mRNA levels of FOSL2 and CEBPB in BJY chicken HMG and LMG groups.

In vitro Validation of FOSL2 Gene

We generated chicken fibroblast DF1 cells overexpressing FOSL2 after lentiviral transfection, and the efficiency of FOSL2 overexpression in DF1 cells was confirmed by fluorescence microscopy of transfected cells (Figures 3A,C). Enhanced green fluorescent protein (EGFP) is a mutant of green fluorescent protein (GFP), which emits more than six times of fluorescence intensity than GFP. EGFP plays a role only after gene expression. In addition, FOSL2-overexpressing DF1 cells had lower glycogen content compared with control group cells (Figure 3B). This is consistent with the results of our previous transcriptome analysis and WGCNA data analysis to mine FOSL2, which predicted their reduced MG content.

FIGURE 3

Figure 3. (A,B) Bright-field images and fluorescence microscopy images of transfected cells of chickens from the control group, respectively. Bright-field images and fluorescence microscopy images of cells from the FOSL2-transfected group, respectively. More than 80% of DF1 cells expressed green fluorescent protein. Scale bar = 100 μm. (B) The glycogen of control and FOSL2 overexpressing DF1 cells (P = 0.046). Control refers to the control group, and FOSL2 refers to the overexpression group of FOSL2. (C) Expression level of FOSL2 gene by qRT-PCR in the control group and FOSL2 overexpression group (P < 0.001). Data are presented as the mean ± SEM. Control refers to the control group, and FOSL2 refers to the overexpression group of FOSL2.RQ: relative quantification.

Discussion

In the past few decades, livestock and poultry genetic breeding programs have focused on economically important traits, and significant progress has been made in genetic improvement. As an important intermediate in the process of sugar metabolism, MG interacts closely with blood sugar through glycogen synthesis. MG is not only involved in maintaining normal life activities of animals, but also affects the production performance and meat quality of livestock and poultry (Le Bihan-Duval et al., 2001; Berri et al., 2007; Alnahhas et al., 2014). However, the key genes that regulate the metabolism of MG in poultry have not yet been identified.

Numerous DEGs were identified in the HMG and LMG chicken groups. Among them, genes coding for proteins involved in the classical glycogen metabolism pathway, such as 5′-AMP–activated protein kinase (AMPK), which is a heterotrimeric protein that is involved in important signaling cascades to regulate cellular energy metabolism (Gowans and Hardie, 2014; Hardie and Ashford, 2014), Although this role of AMPK in glycogen metabolism has been widely established, increasing evidence suggests that the chronic activation of AMPK leads to glycogen accumulation instead of glycogen degradation (Hunter et al., 2011). For instance, it has been shown that the constitutive activation of AMPK via mutations in the γ2 and γ3 subunits leads to the accumulation of glycogen in the skeletal and cardiac muscles of pigs and mice (Milan et al., 2000; Ahmad et al., 2005; Luptak et al., 2007). According to the GeneCards² database, the AMPK gamma subunit (PRKAG3), phosphatidylinositol-4,5-bisphosphate 3-kinase catalytic subunit beta (PIK3CB), protein phosphatase 2 regulatory subunit Bdelta (PPP2R2D), forkhead box O1 (FOXO1), and forkhead box O3 (FOXO3) are involved in the AMPK signaling pathway, thereby regulating MG content.

As expected, it was found that the expression of FOSL2 and CEBPB genes in the breast muscle tissue of the HMG group is different from that of the LMG breast muscle tissue. The information in GeneCards (see text footnote 2) shows that FOSL2 activates CEBPB transcription in PGE₂-activated osteoblasts. In addition, studies in the abdominal fat line broiler breeds bred by the French Ministry of Agriculture have shown that the downregulation of CEBPB affects MG production through a cAMP-dependent signaling pathway (Sibut et al., 2011). Although there is no relevant literature on the involvement of FOSL2 in the regulation of MG content through activation of CEBPB, the DEGs identified by transcriptome analysis of breast muscle tissue from HMG and LMG groups and the existing research results on MG indicate that FOSL2 may regulate the activation of CEBPB and thereby lead to changes in the MG content.

Some studies have reported that the traditional transcriptome sequencing statistical analysis methods have the drawback of not being able to distinguish confounding factors (Barabási and Oltvai, 2004; Ghaemi et al., 2019), but WGCNA can solve some of the analysis issues caused by confounding factors (Camacho et al., 2018; Emilsson et al., 2018; Mangul et al., 2019). Because of the aforementioned problems in transcriptome analysis, 10 individual chickens were randomly added on the basis of the original six individual chickens for transcriptome sequencing and WGCNA. In the WGCNA of JX-Y chickens, the FOSL2 and CEBPB genes present in the ME midnight blue module are highly negatively correlated with MG traits. Besides the candidate genes FOSL2 and CEBPB, this module also comprises many classic genes related to MG and sugar metabolism. For example, MAP3K14, PPP2CA, SLC38A2, PPP2R5E, IGBP1, etc. The above results indicate that the FOSL2 and CEBPB genes are coexpressed and their protein products may have a coregulatory role. Furthermore, network analysis of the MG-related genes in the ME midnight blue module performed using the STRING software to predict the possible relationship between their proteins revealed a regulatory network between FOSL2 and CEBPB. Based on the above findings, we hypothesize that low expression of FOSL2 will inhibit the transcription of CEBPB, thereby decreasing MG production.

Differences in the selection process and selection intensity will lead to differences in traits of different breeds. We selected BJY chickens of different sex to quantitatively analyze the expression of the FOSL2 and CEBPB genes. The results showed that FOSL2 and CEBPB are also differentially expressed in the HMG and LMG chicken groups. These results not only showed that sex did not affect the effect of FOSL2 and CEBPB on the reduction of MG content, but also that both proteins had an effect in reducing MG content in BJY chickens. Currently, there is no report on whether sex differences affect MG content in poultry, but studies in cattle suggest that sex differences do not affect MG content (Komatsu et al., 2015).

Members of the Fos gene family are attracting increasing attention as they may have multifunctional roles in a variety of physiological processes, including fat metabolism, bone development, and the pathogenesis of diseases, such as cancer. In this study, transcriptome analysis and WGCNA revealed that FOSL2 was related to MG content. Overexpression of FOSL2 in chicken fibroblast DF1 cells also confirmed its important role in decreasing the MG content.

At present, there are only a few studies on the key genes involved in MG content regulation in chicken, but more studies have been performed in pigs, cattle, and other animals. For instance, a major MG-related gene, namely, PHKG1, has been discovered in pigs (Ma et al., 2014). In the JX-Y chicken transcriptome analysis in this study, the PHKG1 gene was not found to be differentially expressed. This may be due to the small sample size and differences between species. Although our study has certain limitations, it still advances the existing knowledge of chicken MG traits and provides a theoretical basis for subsequent research.

Conventional transcriptome sequencing analysis has been widely used in the study of various traits of poultry (Luo et al., 2019; Yang et al., 2020). In this study, we used basic differential expression analysis and gene coexpression network analysis to identify molecular network pathways that are involved in the regulation of glycogen content in chicken breast muscle. Furthermore, combined with the results of the whole-genome association analysis in the previous period, a key gene that affects MG content, namely, FOSL2, was identified. The results of the quantitative analysis of FOSL2 expression in BJY chickens showed that its regulatory effect on the MG content is not affected by breed or sex. The results of the analysis in DF1 cells further confirmed the negative regulatory effect of FOSL2 on MG content.

Conclusion

Based on results from previous GWAS and WGCNA, we identified a key gene, namely, FOSL2, which negatively regulates MG content. Our data suggest that FOSL2 may decrease MG content by regulating its downstream gene CEBPB, as both genes exhibit the same pattern of expression. In vivo experiments showed that the effect of FOSL2 on MG was not affected by breed or sex. These results provide new insights and reference for the study of the regulation of MG metabolism.

Data Availability Statement

The raw whole genome sequencing data reported in this paper have been deposited in the Genome Sequence Archive (Wang et al., 2017) in BIG  National Genomics Data Center Members and Partners (2020) under accession number CRA004003 that can be publicly accessed at https://bigd.big.ac.cn/gsa.

Ethics Statement

All animals and experimental protocols used in this study were approved by the Beijing Institute of Animal Science, Chinese Academy of Agricultural Sciences (the scientific research department responsible for animal welfare issues) (No. IAS2019-21).

Author Contributions

XL performed the study, analyzed the data, and drafted the manuscript. LL performed the study and drafted the manuscript. JWa analyzed the data and modified the manuscript. JWe, GZ, and HC contributed to the design of the study. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by grants from the National Natural Science Foundation of China (31872340), the Agricultural Science and Technology Innovation Program (CAAS-ZDRW202005), the Basic Research from Chinese Academy of Agricultural Sciences (Y2019XK06), Basic Research from Institute of Animal Sciences of Chinese Academy of Agricultural Sciences (2019-YWF-YB-07), the Agricultural Science and Technology Innovation Program (ASTIP-IAS04), the Earmarked Fund for Modern Agro-Industry Technology Research System (CARS-41).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2021.682441/full#supplementary-material

Footnotes

1. ^ https://www.string-db.org/
2. ^ https://www.genecards.org/

References

Ahmad, F., Arad, M., Musi, N., He, H., Wolf, C., Branco, D., et al. (2005). Increased alpha2 subunit-associated AMPK activity and PRKAG2 cardiomyopathy. Circulation 112, 3140–3148. doi: 10.1161/circulationaha.105.550806

PubMed Abstract | CrossRef Full Text | Google Scholar

Alnahhas, N., Berri, C., Boulay, M., Baéza, E., Jégo, Y., Baumard, Y., et al. (2014). Selecting broiler chickens for ultimate pH of breast muscle: analysis of divergent selection experiment and phenotypic consequences on meat quality, growth, and body composition traits. J. Anim. Sci. 92, 3816–3824. doi: 10.2527/jas.2014-7597

PubMed Abstract | CrossRef Full Text | Google Scholar

Alnahhas, N., Le Bihan-Duval, E., Baéza, E., Chabault, M., Chartrin, P., Bordeau, T., et al. (2015). Impact of divergent selection for ultimate pH of pectoralis major muscle on biochemical, histological, and sensorial attributes of broiler meat. J. Anim. Sci. 93, 4524–4531. doi: 10.2527/jas.2015-9100

PubMed Abstract | CrossRef Full Text | Google Scholar

Barabási, A. L., and Oltvai, Z. N. (2004). Network biology: understanding the cell’s functional organization. Nat. Rev. Genet. 5, 101–113. doi: 10.1038/nrg1272

PubMed Abstract | CrossRef Full Text | Google Scholar

Berri, C., Le Bihan-Duval, E., Debut, M., Santé-Lhoutellier, V., Baéza, E., Gigaud, V., et al. (2007). Consequence of muscle hypertrophy on characteristics of Pectoralis major muscle and breast meat quality of broiler chickens. J. Anim. Sci. 85, 2005–2011. doi: 10.2527/jas.2006-398

PubMed Abstract | CrossRef Full Text | Google Scholar

Bozec, A., Bakiri, L., Jimenez, M., Rosen, E. D., Catalalehnen, P., Schinke, T., et al. (2013). Osteoblast-specific expression of Fra-2/AP-1 controls adiponectin and osteocalcin expression and affects metabolism. J. Cell Sci. 126, 5432–5440.

Google Scholar

Bozec, A., Bakiri, L., Jimenez, M., Schinke, T., Amling, M., and Wagner, E. F. (2010). Fra-2/AP-1 controls bone formation by regulating osteoblast differentiation and collagen production. J. Cell Biol. 190, 1093–1106. doi: 10.1083/jcb.201002111

PubMed Abstract | CrossRef Full Text | Google Scholar

Camacho, D. M., Collins, K. M., Powers, R. K., Costello, J. C., and Collins, J. J. (2018). Next-generation machine learning for biological networks. Cell 173, 1581–1592. doi: 10.1016/j.cell.2018.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, F., Wu, P., Shen, M., He, M., Chen, L., Qiu, C., et al. (2019). Transcriptome analysis of differentially expressed genes related to the growth and development of the jinghai yellow chicken. Genes (Basel) 10:539. doi: 10.3390/genes10070539

PubMed Abstract | CrossRef Full Text | Google Scholar

Chinenov, Y., and Kerppola, T. K. J. O. (2001). Close encounters of many kinds: fos-jun interactions that mediate transcription regulatory specificity. Oncogene 20, 2438–2452. doi: 10.1038/sj.onc.1204385

PubMed Abstract | CrossRef Full Text | Google Scholar

Eferl, R., Hasselblatt, P., Rath, M., Popper, H. H., Zenz, R., Komnenovic, V., et al. (2008). Development of pulmonary fibrosis through a pathway involving the transcription factor Fra-2/AP-1. Proc. Natl. Acad. Sci. U.S.A. 105, 10525–10530. doi: 10.1073/pnas.0801414105

PubMed Abstract | CrossRef Full Text | Google Scholar

Eferl, R., and Wagner, E. F. (2003). AP-1: a double-edged sword in tumorigenesis. Nat. Rev. Cancer 3, 859–868. doi: 10.1038/nrc1209

PubMed Abstract | CrossRef Full Text | Google Scholar

Emilsson, V., Ilkov, M., Lamb, J. R., Finkel, N., Gudmundsson, E. F., Pitts, R., et al. (2018). Co-regulatory networks of human serum proteins link genetics to disease. Science 361, 769–773. doi: 10.1126/science.aaq1327

PubMed Abstract | CrossRef Full Text | Google Scholar

Foletta, V. C., Sonobe, M. H., Suzuki, T., Endo, T., Iba, H., and Cohen, D. R. (1994). Cloning and characterisation of the mouse fra-2 gene. Oncogene 9, 3305–3311.

Google Scholar

Ghaemi, M. S., DiGiulio, D. B., Contrepois, K., Callahan, B., Ngo, T. T. M., Lee-McMullen, B., et al. (2019). Multiomics modeling of the immunome, transcriptome, microbiome, proteome and metabolome adaptations during human pregnancy. Bioinformatics 35, 95–103. doi: 10.1093/bioinformatics/bty537

PubMed Abstract | CrossRef Full Text | Google Scholar

Godfrey, R., and Quinlivan, R. (2016). Skeletal muscle disorders of glycogenolysis and glycolysis. Nat. Rev. Neurol. 12, 393–402. doi: 10.1038/NRNEUROL.2016.75

PubMed Abstract | CrossRef Full Text | Google Scholar

Gowans, G. J., and Hardie, D. G. (2014). AMPK: a cellular energy sensor primarily regulated by AMP. Biochem. Soc. Trans. 42, 71–75. doi: 10.1042/bst20130244

PubMed Abstract | CrossRef Full Text | Google Scholar

Hardie, D. G., and Ashford, M. L. J. P. (2014). AMPK: regulating energy balance at the cellular and whole body levels. Physiology (Bethesda) 29, 99–107. doi: 10.1152/physiol.00050.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Henckel, P. (1996). “Physiology and biochemistry of muscle fibres in poultry,” in Proceedings from 2nd European Poultry Breeders Roundtable 6-8 September 1995, Foulum.

Google Scholar

Hunter, R. W., Treebak, J. T., Wojtaszewski, J. F., and Sakamoto, K. (2011). Molecular mechanism by which AMP-activated protein kinase activation promotes glycogen accumulation in muscle. Diabetes 60, 766–774. doi: 10.2337/db10-1148

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadarmideen, H. N., Watson-Haigh, N. S., and Andronicos, N. M. (2011). Systems biology of ovine intestinal parasite resistance: disease gene modules and biomarkers. Mol. Biosyst. 7, 235–246. doi: 10.1039/c0mb00190b

PubMed Abstract | CrossRef Full Text | Google Scholar

Karreth, F. A., Hoebertz, A., Scheuch, H., Eferl, R., and Wagner, E. F. (2004). The AP1 transcription factor Fra2 is required for efficient cartilage development. Development 131, 5717–5725. doi: 10.1242/dev.01414

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Langmead, B., and Salzberg, S. L. (2015). HISAT: a fast spliced aligner with low memory requirements. Nat. Methods 12, 357–360. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

Kogelman, L. J., Cirera, S., Zhernakova, D. V., Fredholm, M., Franke, L., and Kadarmideen, H. N. (2014). Identification of co-expression gene networks, regulatory genes and pathways for obesity based on adipose tissue RNA Sequencing in a porcine model. BMC Med. Genomics 7:57. doi: 10.1186/1755-8794-7-57

PubMed Abstract | CrossRef Full Text | Google Scholar

Komatsu, T., Shoji, N., Saito, K., and Suzuki, K. (2015). Effects of genetic and environmental factors on muscle glycogen content in Japanese Black cattle. Anim. Sci. J. 85, 793–798. doi: 10.1111/asj.12201

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Bihan-Duval, E., Berri, C., Baeza, E., Millet, N., and Beaumont, C. (2001). Estimation of the genetic parameters of meat characteristics and of their genetic correlations with growth and body composition in an experimental broiler line. Poult. Sci. 80, 839–843. doi: 10.1093/ps/80.7.839

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Bihan-Duval, E., Debut, M., Berri, C. M., Sellier, N., Santé-Lhoutellier, V., Jégo, Y., et al. (2008). Chicken meat quality: genetic variability and relationship with growth and muscle characteristics. BMC Genet. 9:53. doi: 10.1186/1471-2156-9-53

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., and Dewey, C. N. (2011). RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 12:323. doi: 10.1186/1471-2105-12-323

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Li, S., Hu, Y., Cao, G., Wang, S., Rai, P., et al. (2016). The expression level of mRNA, Protein, and DNA methylation status of FOSL2 of Uyghur in XinJiang in Type 2 diabetes. J. Diabetes Res. 2016:5957404.

Google Scholar

Li, Q., Wang, F., Wang, Q., Zhang, N., Zheng, J., Zheng, M., et al. (2020). SPOP promotes ubiquitination and degradation of MyD88 to suppress the innate immune response. PLoS Pathog. 16:e1008188. doi: 10.1371/journal.ppat.1008188

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Jing, L., and Tu, X. (2016). Weighted gene co-expression network analysis identifies specific modules and hub genes related to coronary artery disease. BMC Cardiovasc. Disord. 16:54. doi: 10.1186/s12872-016-0217-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Liu, X., Cui, H., Liu, R., Zhao, G., and Wen, J. (2019). Transcriptional insights into key genes and pathways controlling muscle lipid metabolism in broiler chickens. BMC Genomics 20:863. doi: 10.1186/s12864-019-6221-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Liu, L., Wang, J., Cui, H., Chu, H., Bi, H., et al. (2020). Genome-Wide Association Study of Muscle Glycogen in Jingxing Yellow Chicken. Genes 11:497. doi: 10.3390/GENES11050497

PubMed Abstract | CrossRef Full Text | Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550–550. doi: 10.1186/S13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, J. J., Zhang, Y., Sun, H., Wei, J. T., Khalil, M. M., Wang, Y. W., et al. (2019). The response of glandular gastric transcriptome to T-2 toxin in chicks. Food Chem. Toxicol. 132:110658. doi: 10.1016/j.fct.2019.110658

PubMed Abstract | CrossRef Full Text | Google Scholar

Luptak, I., Shen, M., He, H., Hirshman, M. F., Musi, N., Goodyear, L. J., et al. (2007). Aberrant activation of AMP-activated protein kinase remodels metabolic network in favor of cardiac glycogen storage. J. Clin. Invest. 117, 1432–1439. doi: 10.1172/jci30658

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, J., Yang, J., Zhou, L., Ren, J., Liu, X., Zhang, H., et al. (2014). A splice mutation in the PHKG1 gene causes high glycogen content and low meat quality in pig skeletal muscle. PLoS Genet. 10:e1004710. doi: 10.1371/journal.pgen.1004710

PubMed Abstract | CrossRef Full Text | Google Scholar

Mangul, S., Martin, L. S., Hill, B. L., Lam, A. K., Distler, M. G., Zelikovsky, A., et al. (2019). Systematic benchmarking of omics computational tools. Nat. Commun. 10:1393. doi: 10.1038/s41467-019-09406-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Milan, D., Jeon, J. T., Looft, C., Amarger, V., Robic, A., Thelander, M., et al. (2000). A mutation in PRKAG3 associated with excess glycogen content in pig skeletal muscle. Science 288, 1248–1251. doi: 10.1126/science.288.5469.1248

PubMed Abstract | CrossRef Full Text | Google Scholar

National Genomics Data Center Members and Partners (2020). Database resources of the national genomics data center in 2020. Nucleic Acids Res. 48, D24–D33. doi: 10.1093/nar/gkz913

PubMed Abstract | CrossRef Full Text | Google Scholar

Puolanne, E., and Immonen, K. (2014). CONVERSION of Muscle to Meat. Glycogen. Oxford: Academic Press.

Google Scholar

Sibut, V., Hennequet-Antier, C., Bihan-Duval, E. L., Marthey, S., Duclos, M. J., and Berri, C. (2011). Identification of differentially expressed genes in chickens differing in muscle glycogen content and meat quality. BMC Genom. 12:112. doi: 10.1186/1471-2164-12-112

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Franceschini, A., Wyder, S., Forslund, K., Heller, D., Huerta-Cepas, J., et al. (2015). STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 43, D447–D452. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Song, F., Zhu, J., Zhang, S., Yang, Y., Chen, T., et al. (2017). GSA: genome sequence archive. Genom. Proteom. Bioinform. 15, 14–18. doi: 10.1016/j.gpb.2017.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, L., Zheng, X., Mo, C., Li, S., Liu, Z., Yang, G., et al. (2020). Transcriptome analysis and identification of genes associated with chicken sperm storage duration. Poult. Sci. 99, 1199–1208. doi: 10.1016/j.psj.2019.10.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, B., and Horvath, S. (2005). A general framework for weighted gene co-expression network analysis. Stat. Appl. Genet. Mol. Biol. 4:17. doi: 10.2202/1544-6115.1128

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, W., Langfelder, P., Fuller, T., Dong, J., and Hovarth, S. (2010). Weighted gene coexpression network analysis: state of the art. J. Biopharm. Stat. 20, 281–300. doi: 10.1080/10543400903572753

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, X., Sun, K., Lan, Z., Song, W., Cheng, L., Chi, W., et al. (2017). Tenofovir and adefovir down-regulate mitochondrial chaperone TRAP1 and succinate dehydrogenase subunit B to metabolically reprogram glucose metabolism and induce nephrotoxicity. Sci. Rep. 7:46344. doi: 10.1038/srep46344

PubMed Abstract | CrossRef Full Text | Google Scholar