Transcriptome-based analysis of key genes and pathways affecting the linoleic acid content in chickens

ZHAO Wen-juan , YUAN Xiao-ya, XIANG Hai MA Zheng CUI Huan-xian, LI Hua#, ZHAO Gui-ping #

1 Guangdong Provincial Key Laboratory of Animal Molecular Design and Precise Breeding, Foshan University, Foshan 528000,P.R.China

2 Institute of Animal Science, Chinese Academy of Agricultural Sciences, Beijing 100193, P.R.China

Abstract Linoleic acid is an essential polyunsaturated fatty acid that cannot be synthesized by humans or animals themselves and can only be obtained externally.The amount of linoleic acid present has an impact on the quality and flavour of meat and indirectly affects consumer preference.However, the molecular mechanisms influencing the deposition of linoleic acid in organisms are not clear.As the molecular mechanisms of linoleic acid deposition are not well understood, to investigate the main effector genes affecting the linoleic acid content, this study aimed to screen for hub genes in slow-type yellow-feathered chickens by transcriptome sequencing (RNA-Seq) and weighted gene coexpression network analysis (WGCNA).We screened for candidate genes associated with the linoleic acid content in slow-type yellow-feathered broilers.A total of 399 Tiannong partridge chickens were slaughtered at 126 days of age, fatty acid levels were measured in pectoral muscle, and pectoral muscle tissue was collected for transcriptome sequencing.Transcriptome sequencing results were combined with phenotypes for WGCNA to screen for candidate genes.KEGG enrichment analysis was also performed on the genes that were significantly enriched in the modules with the highest correlation.A total of 13 310 genes were identified after quality control of transcriptomic data from 399 pectoral muscle tissues.WGCNA was performed, and a total of 26 modules were obtained, eight of which were highly correlated with the linoleic acid content.Four key genes, namely, MDH2,ATP5B, RPL7A and PDGFRA, were screened according to the criteria |GS|>0.2 and |MM|>0.8.The functional enrichment results showed that the genes within the target modules were mainly enriched in metabolic pathways.In this study, a large-sample-size transcriptome analysis revealed that metabolic pathways play an important role in the regulation of the linoleic acid content in Tiannong partridge chickens, and MDH2, ATP5B, RPL7A and PDGFRA were screened as important candidate genes affecting the linoleic acid content.The results of this study provide a theoretical basis for selecting molecular markers and comprehensively understanding the molecular mechanism affecting the linoleic acid content in muscle, providing an important reference for the breeding of slow-type yellowfeathered broiler chickens.

Keywords: chicken, linoleic acid, transcriptome sequencing, weighted gene coexpression network analysis (WGCNA),metabolic pathways

1.Introduction

Chicken is the second most consumed meat in China, and the import volume and production volume are increasing annually.With the transformation of consumer preference from “quantity” to “quality”, the taste, aroma and nutritional quality of chicken affect consumers’ desire to buy.Yellowfeathered broilers are attractive to consumers due to their excellent meat quality and flavour.Among them,Tiannong partridge chicken, as a high-quality local yellow chicken breed, has the quality of both broiler chickens and laying hens, with tender meat and unique soft bone flavour, occupying a certain market in Guangdong, Hong Kong and Macao of China (Guonan 2006).

Poultry products are one of the sources of linoleic acid for humans, and there have been relevant studies on several species, but few studies have been conducted on the Tiannong partridge chicken (Cartoniet al.2022).Linoleic acid is a polyunsaturated essential fatty acid that cannot be synthesized by humans and animals themselves and can only be obtained externally.At the same time, linoleic acid can be used as a source of energy; it undergoes beta-oxidation as a long-chain fatty acid to produce the energy required by the body(Lopaschuket al.2010).As a structural component,linoleic acid can be used as a reference when determining cardiovascular health and as a marker for the uniform distribution of nanovesicles in bionic cell membranes(Letko Khaitet al.2019; Lázaroet al.2021).The production of the biologically active substance arachidonic acid is involved in the regulation of immune function,and when the dietary linoleic acid content is reduced,arachidonic acid levels are reduced, thereby reducing the activation of endotoxic lipopolysaccharides in the body(Martinet al.2016; Innesa and Calder 2018).Linoleic acid is a precursor of volatile organic compounds (VOCs),and changing the concentration of conjugated linoleic acid can change the concentration of VOCs (Niet al.2020); the linoleic acid content is significantly correlated with the volatile content (Yuanet al.2022).Moreover,chicken meat can be used as a source of linoleic acid for humans, having an impact on human diseases such as atherosclerosis and cancer (Whelan and Fritsche 2013; Choqueet al.2014; Calder 2015).The addition of linoleic acid in meat production can improve the quality of meat products and have health benefits (Sunet al.2004).The conversion and master effector genes of linoleic acid in animals are currently unknown.

Only a few studies have shown that linoleic acid is associated with metabolic pathways.Recent reports suggest that phospholipase C delta 1 (PLCD1) may regulate the linoleic acid content and be involved in metabolic pathways, as shown through genome-wide association analysis (GWAS) (Liuet al.2022).The genes of fatty acid synthesis- and degradation-related enzymes such as fatty acid desaturase 1 (FADS1) and fatty acid synthase (FASN) have been screened by weighted gene coexpression network analysis (WGCNA)and by the identification of differentially expressed genes(DEGs) related to fatty acids in pectoral muscles of bulls(Schettiniet al.2022).Meanwhile, some scholars have subjected volatiles to WGCNA, revealing that the modules associated with floral fragrance are related to fatty acid derivatives (Shiet al.2018).A previous study found by WGCNA that genes in the mitogen-activated protein kinase (MAPK) signalling pathway may affect fatty acid composition in intramuscular fat (Zappaterraet al.2021).Therefore, WGCNA can be used to screen for genes that affect the linoleic acid content.To date, studies have only focused on genes affecting the synthesis and degradation of fatty acids; however, few genes affecting the linoleic acid content in chicken meat at the molecular level have been identified.

At the same time, linoleic acid is an essential fatty acid and has a positive effect on nutrition and health, but the molecular mechanisms that influence the deposition of linoleic acid in chicken meat are not well understood.Therefore, we determined fatty acid content in breast muscle tissues of 399 126-day-old Tiannong partridge chickens, sequenced the tissues, and identified candidate genes involved in the linoleic acid content in depth using WGCNA.The accuracy of the findings will be improved by large-scale transcriptome sequencing.It will help forthcoming research on regulation of mechanisms on linoleic acid deposition and metabolism in chicken.

2.Materials and methods

2.1.Experimental animals and sample collection

The experimental material was Tiannong partridge female chickens provided by Guangdong Tinoo’s Foods Group Co., China.The experimental animals were reared in individual cages under the same environmental and nutritional conditions with free access to water and food.The rearing standard was based on the Chicken Feeding Standard (NY/T33-2004) as a reference, with an ambient temperature of 15-27°C, relative humidity of not less than 50% and an average light duration of 16 h.

The animals were reared to 126 days of age, and 399 were randomly selected for slaughter.Samples were collected uniformly from the same location on the left pectoral muscle, rapidly frozen in liquid nitrogen, and stored at -80°C for backup.The fat and fascia were stripped from one side of the pectoral muscle, and the samples were churned on ice using a meat grinder and stored at -80°C for subsequent fatty acid determination.

2.2.Determination of fatty acid levels in the pectoralis muscle by GC

The reagents required for the determination of fatty acid levels were prepared as follows: 1) Chloroacetyl methanol solution: 10 mL of methanol solution was pipetted out, and them 100 mL of chloroacetyl solution was added dropwise, followed by gentle stirring and shaking well.The solution was transferred to a reagent bottle for storage.2) Potassium carbonate solution (7%):7 g of potassium carbonate was dissolved in ultrapure water, gently stirred and shaken well.The solution was then transferred to a 100 mL volumetric flask for volume determination and stored in a reagent bottle.

The fatty acid content in tissues was determined using gas chromatographic (GC) methods (Chiu and Kuo 2020).The 6 g of pectoral muscle was lyophilized at low temperature and ground to powder.Then, 200 mg of lyophilized chicken powder was weighed out and added with 4 mL of chloroacetylmethanol solution in a hydrolysis tube.Then, 1 mL (1.0 mg mL-1) of methyl undecanoate(Sigma, USA) was added as the internal standard,followed by adding 1 mL of hexane (Sigma, USA).The solution was flushed with the inert gas nitrogen for protection, the bottle was capped tightly, and the contents were mixed by shaking.The organic layer was filtered through a 0.2 μm filter membrane and loaded into a sample bottle.The determination was carried out using a gas chromatograph (Agilent, USA).

2.3.RNA extraction and sequencing

RNA was extracted from pectoral muscle samples using the TRIzol method, RNA purity was assessed using a Kaiao K5500 spectrophotometer (Kaiao, China), and RNA integrity and concentration were assessed using an RNA Nano 6000 Assay Kit and a Bioanalyzer 2100 System (Agilent, USA).The A260/A280value ranged from 1.8 to 2.0.RNA samples with RNA integrity>7.5 were used for subsequent sequencing.Sequencing was performed on the Illumina NovaSeq 6000 (Illumina,USA) platform, and 150 bp paired-end reads were generated (Kanget al.2021).The sequenced data were used to check the quality of RNA sequences using Fast QC software (v0.11.9) for raw data, sequencing junction sequences were removed, and reads with low quality and indeterminate base information proportions greater than 5% were filtered out to obtain clean reads.The sequenced reads were matched to the index using STAR software (v2.7.10a) based on the chicken reference genome (Ensembl GRCg6a(GCA_000002315.5)), and the genes and transcripts were quantified using RSEM software (v1.3.3) (Liet al.2011; Ji and Sadreyev 2018).The subsequent series of analyses were performed on this basis.All RNA sequencing data passed quality control prior to analysis(Kanget al.2021).

2.4.Weighted gene coexpression network analysis

WGCNA was performed using the WGCNA package in R.Gene expression data were first filtered to screen for fragments per thousand bases of transcripts per million mapped reads (FPKM) values, mapping read data for genes with a median absolute deviation (MAD)of 75% and >0.01.The samples were then clustered,and outliers were removed using the hcluster function.Based on the principle of optimal scale-free topology, a soft threshold power of 20 was chosen.A new distance matrix was then obtained by converting the adjacency matrix into a topological overlap matrix (TOM) for noise and pseudocorrelation reduction.Then, 1-TOM was calculated as a biologically important measure of the degree of network interoperation and used as a distance metric for gene hierarchical clustering.The dynamic tree cutting algorithm was used to identify coexpressed gene modules, calculate the module eigengene (ME),and merge the modules marked with different colours.Finally, specific modules were identified based on module(ME)-phenotype correlation.Determining module membership (MM) is required to screen for pivotal genes and thus measure the relationship between genes and specific modules (Langfelder and Horvath 2008).In this study, we used the |GS| (representing the correlation between the gene and a given clinical trait)>0.2 and|MM| (representing the correlation between the gene and a given module)>0.8 in the WGCNA package to directly identify hub genes.

2.5.KEGG enrichment pathway analysis

Kyoto Encyclopedia of Genes and Genomes(KEGG) pathway enrichment analysis was performed using genes withP<0.05 in the target module (MEdarkmagenta) from the previous WGCNA screen.The analysis was carried out using the KOBAS database (http://kobas.cbi.pku.edu.cn/).The enriched pathways were screened at a significance level ofP<0.05.

2.6.Statistical analysis of data

The mean and standard deviation (SD) of the linoleic acid content in Tiannong partridge chicken breast muscle were analysed using Microsoft Excel 2019 (Microsoft Corp., Redmond, WA, USA).Statistical significance (P-value) was calculated using at-test, with an apparent significant difference noted whenPwas less than 0.05.

3.Results

3.1.Key module weighted gene coexpression network construction and identification

Transcriptome data were quality-controlled and screened for genes with a median absolute deviation of the top 75%, with a minimum MAD greater than 0.01.A total of 13 310 genes were obtained for subsequent transcriptome analysis.Clustering was performed on 399 samples (cuthighting=28 000)(Fig.1-A).The relationship with soft-thresholding power was calculated using the scale-free topology fit index, settingR2>0.85, resulting in an optimal soft threshold sft=7 (Fig.1-B).Gene expression was clustered according to different samples to obtain 29 coexpressed gene modules (Fig.2), and those with correlation>0.75 were merged to finally obtain 26 correlated modules.

Fig.2 Gene clustering tree (dendrogram).

3.2.Selection and analysis of target modules

To obtain coexpression modules related to the linoleic acid content, correlation analysis was performed on the phenotype matrix to obtain the correlation coefficients andP-values of each module with the linoleic acid content (Fig.3).The modules withP<0.05 were selected as the target modules, and eight target modules were finally identified, namely, the “MEcyan” module,“MEgreen” module, “MEbrown” module, “MEblack”module, “MEmidnightblue” module, “MEdarkgray” module,“MEplum1” module and “MEdarkmagenta” module.The numbers of genes significantly associated with the linoleic acid content in these eight modules were 57, 377, 202,96, 38, 141, 10 and 30, respectively (Appendix A).The genes within the relevant modules were then plotted for identification of significant correlation with the linoleic acid content (Fig.4).

Fig.3 Heatmap of linoleic acid and module correlation.

Fig.4 Correlation chart: Correlation between module membership (MM)and gene significance (GS).The correlation between the linoleic acid content and the module of interest is indicated.A, “MEcyan” module.B,“MEgreen” module.C, “MEbrown”module.D, “MEblack” module.E, “MEmidnightblue” module.F,“MEdarkgray” module.G, “MEplum1”module.H, “MEdarkmagenta”module.

Among them, the “MEcyan” module, “MEgreen”module, “MEbrown” module, “MEblack” module,“MEmidnightblue” module and “MEdarkgray” module were positively correlated (eigengene value=0.16,eigengene value=0.14, eigengene value=0.14,eigengene value=0.14, eigengene value=0.14,eigengene value=0.13).The “MEplum1” module and the“MEdarkmagenta” module were both negatively correlated(eigengene value=-0.14, eigengene value=-0.2).

Four pivotal genes were screened from all the significantly correlated modules using the criteria |GS|>0.2 and |MM|>0.8.Among them, malate dehydrogenase 2(MDH2) and ATP synthase, H+transporting mitochondrial F1complex, beta subunit (ATP5B) were located in the“MEdarkmagenta” module, ribosomal protein L7a (RPL7A)in the “MEcyan” module, and platelet derived growth factor receptor alpha (PDGFRA) in the “MEblack” module.The “MEdarkmagenta” module had the highest absolute value of the correlation coefficient, so the subsequent analysis was based on this module.

3.3.KEGG enrichment analysis of related module genes

KEGG enrichment analysis of the target modules yielded a total of 20 related pathways (Fig.5; Appendix B), and 13 significantly related metabolic pathways were obtained,including oxidative phosphorylation, metabolic pathways,the citrate (TCA) cycle, carbon metabolism, pyruvate metabolism, glyoxylate and dicarboxylate metabolism,cysteine and methionine metabolism, cardiac muscle contraction, glycolysis/gluconeogenesis, phenylalanine,tyrosine and tryptophan biosynthesis, 2-oxocarboxylic acid metabolism, phenylalanine metabolism, and arginine biosynthesis.

Fig.5 Correlation module pathway enrichment analysis.

Most genes were enriched in metabolic pathways,followed by the ribosome pathway, oxidative phosphorylation pathway, TCA cycle pathway, carbon metabolism and pyruvate metabolism pathway.

4.Discussion

Chicken is one of sources of linoleic acid for humans.There are many parameters to assess the quality and flavour of chicken meat, among which fatty acids are among the important parameters used to assess the quality of chicken meat and one of the flavour precursors of chicken meat (Woodet al.2008; Alagawanyet al.2019; Zaunschirmet al.2019; Jinet al.2021).Tiannong partridge chicken is a high-quality local breed in China, so this experiment was conducted at the transcriptome level to mining candidate genes and pathways related to the linoleic acid content in local chickens by using nearly 400 individuals.In this study, we first performed a weighted gene coexpression network analysis of transcripts from pectoral muscle tissue and linoleic acid data from pectoral muscle tissue of Tiannong partridge chicken.A total of 8 significantly related modules were identified, and 4 hub genes,MDH2,ATP5B,RPL7AandPDGFRA, were identified.

Linoleic acid affects metabolism and fat production in animals.Linoleic acid has been reported to modulate fatty acid metabolism and glucose metabolism to regulate glucose homeostasis (Hamilton and Klett 2021).It undergoes beta-oxidation to produce acetyl-CoA for acetyl-CoA transport in the citric-pyruvic acid cycle,which enters the citric acid cycle and is oxidized to CO2and H2O to provide sufficient energy for the organism(Sergielet al.2001; Bennettet al.2020).In myocytes,muscle contraction requires energy, and calcium ions promote TCA cycle (Wanet al.1989).Thus, TCA cycle and the synthesis and catabolism of fatty acids in muscle interact and complement each other.The finding that supplementation with conjugated linoleic acid in cattle provides the body with additional energy is consistent with our finding that linoleic acid is primarily involved in metabolic pathways (Vogelet al.2021).Dietary supplementation with different fatty acids for pigs showed that linoleic acid was most closely involved in metabolism and enriched in phenylalanine metabolism, which is also consistent with the results of this study (Xiet al.2021).It has been found that dietary supplementation with conjugated linoleic acid modulates low-density lipoprotein receptor (LDLR) expression and inhibits hepatic 3-hydroxy-3-methylglutaryl coenzyme A reductase(HMGR) and cholesterol 7 alpha hydroxylase 1 (CYP7A1)expression, thus further affecting lipid deposition (Wang S Het al.2019).An analysis of fatty acid and transcriptome data in duck pectoral muscle tissue revealed that some genes were enriched in the oxidative phosphorylation pathway, which is consistent with our findings (Fanet al.2020).In this study, we identified the pathways oxidative phosphorylation, TCA cycle and carbon metabolism by enrichment analysis of the genes in the target module,which is consistent with the results of previous studies,and we can infer that linoleic acid may play a role in the life activities of animals.

MDH2is involved in metabolic pathways such as TCA cycle, glyoxylate and dicarboxylate metabolism, pyruvate metabolism, cysteine and methionine metabolism, and carbon metabolism.Expression and translation ofMDH2increases under high unsaturated fatty acid conditions(Liuet al.2021) and is higher in the liver than in adipose tissue (Wang G Set al.2019), whereMDH2is highly acetylated, thereby affecting glycolysis (Wanget al.2018).A high-fat diet with palmitic acid reduced gluconeogenesis and sirtuin 3 (SIRT3) gene expression and thus MDH2,phosphoenolpyruvate carboxykinase (PEPCK) and pyruvate carboxylase (PC) enzyme activity (Guoet al.2022).A combined metabolomics and proteomics analysis revealed that arachidonic acid (Aramchol) acts on glucose metabolism and lipid metabolism and activates beta-oxidation and oxidative phosphorylation of fatty acids (Fernandez-Ramoset al.2020).Meanwhile,MDH2has been used as a molecular marker for screening fat shape-related molecules, providing a theoretical basis for molecular breeding (Zhouet al.2012).In summary,MDH2plays an important role in the metabolism of linoleic acid and in regulating homeostasis in animals.

ATP5Bhas a tendency to negatively regulate fatty acid degradation and acts synergistically with oxidative phosphorylation.It plays a role in metabolism in the organism.It has been shown that the oxidation of fatty acids through the addition of soy glycosides is associated withATP5Band can reduce the amount of intramuscular fat (Kitamuraet al.2020).Simultaneous feeding of ractopamine reduced energy metabolism and increased fatty acid degradation, while leading to a negative correlation between ATP5B enzyme expression and acyl-CoA dehydrogenase very long chain (ACADVL) (Zhaiet al.2022).In summary,ATP5Bexhibits a negative regulatory trend with fatty acid degradation and acts synergistically with oxidative phosphorylation.It plays an important role in metabolism, especially in energy metabolism and lipid metabolism.ATP5Bexpression was negatively correlated with the linoleic acid content in the results of this study, in line with previous research (Liuet al.2022).

PDGFRAplays a role in many life activities, such as embryonic development, cell differentiation and proliferation.It was enriched in 6 pathways, including the gap junction pathway, calcium signalling pathway, focal adhesion, regulation of actin cytoskeleton, endocytosis and MAPK signalling pathway.Fibroblasts enriched with this gene appeared as lipid droplets after lipid induction in a model of skeletal muscle injury constructed by injection of glycerol into skeletal muscle and had the potential to produce lipids (Xuet al.2021).Heat stress affects fat deposition and meat quality in birds with a corresponding change inPDGFRAexpression (Guet al.2020; Malilaet al.2022).Excessive fat accumulation in liver tissues of chickens resulted in hepatic steatohepatic haemorrhage,and transcriptome sequencing was performed to identifyPDGFRAas a differentially expressed gene between the control and pathological groups (Zhuet al.2021).

KEGG enrichment analysis of the genes associated with the linoleic acid content in “MEdarkmagenta” module from the WGCNA results also revealed enrichment of oxidative phosphorylation, metabolic pathways, TCA cycle, gluconeogenesis, cardiac muscle contraction and other signalling pathways.There were three overlapping signalling pathways in the KEGG enrichment analysis with the previously screened hub genes, namely, TCA cycle,pyruvate metabolism, and glyoxylate and dicarboxylate metabolism.This result suggests that ingested linoleic acid may be involved in the metabolic and vital activities in animals.Linoleic acid is inseparable from metabolic pathways in animals and is converted to fat, energy or flavourviaTCA cycle, pyruvate metabolism, and other pathways.For example, the heart requires a large amount of energy to beat, and linoleic acid can assist in the beta-oxidation process (Funakoshiet al.2021).Second, polyunsaturated fatty acids can inhibit the active contraction of the heart muscle and thus prevent heart disease (Kang and Leaf 1994).It has been shown in numerous studies by previous researchers that linoleic acid, as a fatty acid, plays an important role in the life activities of animals.

5.Conclusion

In this study, we performed transcriptome sequencing on nearly 400 samples.Eight modules associated with linoleic acid were identified using weighted gene coexpression network analysis.Using the common screening criteria |GS|>0.2 and |MM|>0.8 for the hub genes,MDH2,ATP5B,RPL7AandPDGFRAwere identified.Interestingly, bothMDH2andATP5Bwere found to be located in metabolic pathways.These findings will provide candidate genes for further transformation of linoleic acid at the genomic level.

Acknowledgements

This study was supported by the China Agriculture Research System of MOF and MARA (CARS-41),the Key-Area Research and Development Program of Guangdong Province, China (2020B020222002),the Foshan University High-level Talent Program,China (CGZ07243), the Guangdong Provincial Key Laboratory of Animal Molecular Design and Precise Breeding, China (2019B030301010), the Key Laboratory of Animal Molecular Design and Precise Breeding of Guangdong Higher Education Institutes, China(2019KSYS011), and the Foshan Institute of Science and Technology Postgraduate Free Exploration Fund, China(2021ZYTS36).

Declaration of competing interest

The authors declare that they have no conflicts of interest.

Ethical approval

All animal experiments in this study were approved by the Laboratory Animal Welfare and Animal Experimentation Ethics Monitoring Committee of Foshan University (No.18091801), and were performed in strict accordance with the regulations for the use of laboratory animals.The experimental chickens were housed in ordinary rearing facilities, which conformed to the requirements set out in the national standard “Requirements for Laboratory Animal Environment and Rearing Facilities”(GB14925-2001).The experimental animal care and animal experimentation protocols and conditions were in accordance with the “Management Measures for the Use of Experimental Animals in Guangdong Province”, China.

Appendicesassociated with this paper are available on https://doi.org/10.1016/j.jia.2023.02.019