RNA-Seq reveals differentially expressed genes affecting polyunsaturated fatty acids percentage in the Huangshan Black chicken population

PLOS ONE, Apr 2018

Fatty acids metabolic products determine meat quality in chickens. Identifying genes associated with fatty acids composition could provide valuable information for the complex genetic networks of genes with underlying variations in fatty acids synthesis. RNA sequencing (RNA-Seq) was conducted to explore the chicken transcriptome from the thigh muscle tissue of 6 Huangshan Black Chickens with 3 extremely high and low phenotypic values for percentage of polyunsaturated fatty acids (PUFAs). In total, we obtained 41,139,108–44,901,729 uniquely mapped reads, which covered 74.15% of the current annotated transcripts including 18964 mRNA transcripts, across all the six thigh muscle tissue samples. Of these, we revealed 274 differentially expressed genes (DEGs) with a highly significant correlation with polyunsaturated fatty acids percentage between the comparison groups based on the ratio of PUFA/SFA. Gene ontology and pathway analysis indicated that the DEGs were enriched in particular biological processes affecting fatty acids metabolism, biosynthesis of unsaturated fatty acids (USFAs), and cell junction-related pathways. Integrated interpretation of differential gene expression and formerly reported quantitative trait loci (QTL) demonstrated that FADS2, DCN, FRZB, OGN, PRKAG3, LHFP, CHCHD10, CYTL1, FBLN5, and ADGRD1 are the most promising candidate genes affecting polyunsaturated fatty acids percentage.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:


RNA-Seq reveals differentially expressed genes affecting polyunsaturated fatty acids percentage in the Huangshan Black chicken population

April RNA-Seq reveals differentially expressed genes affecting polyunsaturated fatty acids percentage in the Huangshan Black chicken population Shaohua Yang 0 2 3 4 Ying Wang 0 2 3 4 Lulu Wang 0 2 3 4 Zhaoyuan Shi 0 2 3 4 Xiaoqian Ou 0 2 3 4 Dan Wu 0 2 3 4 Xinmiao Zhang 0 2 3 4 Hao Hu 0 2 3 4 Jia Yuan 0 2 3 4 Wei Wang 0 1 2 4 Fuhu Cao 0 2 3 4 Guoqing Liu 0 2 3 4 0 demonstrated that FADS2 , DCN, FRZB, OGN, PRKAG3, LHFP, CHCHD10, CYTL1 1 Agricultural Products Quality and Safety Supervision and Management Bureau , Xuancheng, Anhui , P. R. China 2 University-Research in Hefei University of Technology XC2016JZBZ03 Dr Shaohua Yang; Fundamental Research Funds for the Central Universities JZ2017HGTA0228 Dr Shaohua Yang; the innovation Project for College Students 3 College of Food Science and Engineering, Hefei University of Technology , Hefei, Anhui , P. R. China 4 Editor: Marinus F.W. te Pas, Wageningen UR Livestock Research , NETHERLANDS Fatty acids metabolic products determine meat quality in chickens. Identifying genes associated with fatty acids composition could provide valuable information for the complex genetic networks of genes with underlying variations in fatty acids synthesis. RNA sequencing (RNA-Seq) was conducted to explore the chicken transcriptome from the thigh muscle tissue of 6 Huangshan Black Chickens with 3 extremely high and low phenotypic values for percentage of polyunsaturated fatty acids (PUFAs). In total, we obtained 41,139,108± 44,901,729 uniquely mapped reads, which covered 74.15% of the current annotated transcripts including 18964 mRNA transcripts, across all the six thigh muscle tissue samples. Of these, we revealed 274 differentially expressed genes (DEGs) with a highly significant correlation with polyunsaturated fatty acids percentage between the comparison groups based on the ratio of PUFA/SFA. Gene ontology and pathway analysis indicated that the DEGs were enriched in particular biological processes affecting fatty acids metabolism, biosynthesis of unsaturated fatty acids (USFAs), and cell junction-related pathways. Integrated interpretation of differential gene expression and formerly reported quantitative trait loci (QTL) FBLN5, and ADGRD1 are the most promising candidate genes affecting polyunsaturated fatty acids percentage. Introduction During recent decades, the breeding of meat type poultry focused on increasing growth performance and improving breast and thigh meat yields. Although the impressive progress made in these meat quality traits, there were some poor performances, such as excessive deposition of abdominal fat, deterioration of taste quality, and decreased sensory acceptability for 2017CXCYS235 Dr Shaohua Yang; Wanjiang Institute of Poultry Technology XC2015AKKG002 Dr Shaohua Yang; the Science and Technology Bureau of AnHui XC2015XKKJ02, XC2015AKKG003 Dr Shaohua Yang. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. consumers. As an indigenous breed in China, the Huangshan Black Chicken has a distinct appearance and quality in meat and egg products. Compared with the Arbor Acres (AA) broiler, the Huangshan Black Chicken is highly popular in China because of its polyunsaturated fatty acids concentration. Therefore, the elucidation of the precise molecular mechanisms underlying fatty acids traits in Huangshan Black chickens will have both economic and biological consequences. In the past several decades, candidate gene analysis, quantitative trait locus (QTL) mapping, and genome-wide association studies (GWASs) have been the main approaches to identify genes or causal mutations for meat quality traits in chickens. A large number of promising genetic associations and genomic regions have been successfully identified. As of December 21, 2017, 8,363 QTLs for 383different traits have been reported in 277 papers in chickens (http://www.animalgenome.org/cgi-bin/QTLdb/GG/index). Of these, 339 QTLs for abdominal fat traits and 144 for breast muscle traits have been detected in the chicken chromosomal regions. Moreover, GWASs can be used to identify the genes and variants underlying important traits more precisely [1±2]. In chickens, GWASs have already been performed for egg production and quality [ 3 ], growth [4±5], meat quality [6±7] and disease resistance traits [8±9]. Although these techniques have contributed significantly to our better understanding of mechanisms underlying meat quality traits [ 6,7 ], several potential limitations still exist. One major limitation is the fine mapping required to identify causative variants. Additionally, some novel genes or biological pathways associated with the target trait may be excluded unintentionally. In recent years, next generation sequencing (NGS) technologies are increasingly being used for identifying differential expression as well as opportunities to explore novel transcripts [ 10 ]. Of these, RNA-Seq has been widely utilized to detect differentially expressed genes (DEGs) between two gene expression patterns, causative variants, and alternative splicing events. In chickens, many studies of RNA-Seq have been conducted using intestinal mucosal [ 11 ], heart [ 12 ], uterine [ 13 ], and ovarian tissues [ 14 ]. However, limited studies on the transcriptome of thigh muscle tissues have been reported. The identification of DEGs in thigh muscle tissue represents the first step toward clarifying the complex biological properties of meat quality traits. Therefore, the regulation of fat deposition in chickens at a genome wide level remains to be elucidated. In the present study, we used RNA-Seq technology to examine the genome-wide transcription profile in thigh muscle tissues between two groups of Huangshan Black Chickens with extremely high and low phenotypic values of polyunsaturated fatty acids. We then proposed key candidate genes affecting polyunsaturated fatty acids percentage by conducting integrated analysis. The identified candidate genes could lead to improved selection of chicken while providing new insights into meat quality traits. Materials and methods Ethics statement All procedures for animal handling were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) of Hefei University of Technology (Permit Number: DK838). Animals diet and sample collection Huangshan Black Chickens were maintained in free-ranging flocks in a standardized farm (Anhui conservation farm for Huangshan Black Chicken, Huangshan, China), using a diet as: maize 64.0%, wheat bran 16.0%, full-fat soybean 10.0%, fish meal 5.0%, feed yeast powder 2.0%, bone meal 1.5%, inorganic additives 0.7%, Lysine 0.3%, Methionine 0.2%, salt 0.3%. Ten male chickens with an average weight of 1.82 kg at 120 days old were selected randomly for 2 / 14 our detection. To keep the environment factors identical, all chickens were raised in the same way which was free access to food and water in natural lighting. Performance traits of Huangshan Black chickens were shown in supporting information (S1 File). 12 h after feed was withheld, the chickens were handled by electroshock, bled, and dismembered. The thigh muscle tissue samples from the left leg of the chickens were removed within 30 min after slaughter. The samples for each chicken were carefully collected into an RNasefree tube for RNA isolation, immediately frozen in liquid nitrogen, and kept at −80 ÊC until required for RNA isolation. Meanwhile, sufficient samples were minced and kept at −20 ÊC for fatty acids analysis. Fatty acids analysis Fatty acids samples from ten different individuals were methylated according to Ren et al [ 15 ] with some modifications. 1 mL of 1 M KOH-methanol was added to the lipids (100 μL) for esterification at 65 ÊC for 30 min. Then, after being cooled to room temperature, a 2 mL mixture of boron fluoride-methanol (140 g BF3-ether per liter of methanol) was added to deal with the fatty acids, and then heated at 65 ÊC for 5 min. 1 mL of saturated NaCl solution and 1 mL n-hexane were then added at room temperature. The liquid was allowed to separate into 2 phases using a centrifuge (Thermo Scientific, Wilmington, DE, USA) at 2000 rpm for 5 min. The upper phase containing fatty acid methyl esters (FAMEs) was collected. The FAMEs samples were subsequently analyzed by Gas Chromatography-Mass Spectrometer (GC-MS) using an Agilent 5975C GC (Santa Clara, CA, USA) equipped with a quadruple mass spectrometer (flame ionization detector) and a capillary polar HP-88 cyanopropyl column (60 m × 0.25 mm ID × 0.20 μm film). With a flow rate of 2 mL / min, Helium was used as the carrier gas. Initial column temperature was maintained at 125 ÊC, 8 ÊC per minute from 125 ÊC to 145 ÊC, then raised to 220 ÊC at 2 ÊC / min and maintained for 67 min. Meanwhile, the temperatures of the injector and FID detector were both set at 250 ÊC. As the internal standard, nonadecanoic acid (C19:0) (Sigma, Saint Louis, MO, USA) was used to quantity the fatty acids. The details of the ten samples were detected as shown in supporting information (S2 File). Of these, six samples (polyunsaturated fatty acid high (FAH): FAH1, FAH2, FAH3; polyunsaturated fatty acid low (FAL): FAL1, FAL2, FAL3) were divided into two groups with extremes of the phenotypic values for PUFA/(SFA+USFA) to detect DEGs for sequential analyses. RNA isolation and quality assessment The thigh muscle tissues of six samples were disrupted with liquid nitrogen and total RNA was extracted with TRIzol reagents (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions [ 16 ]. Using the RNase-free DNase I (Invitrogen, Carlsbad, CA, USA), DNA contamination was removed from the RNA by incubating for 30 min at 37 ÊC. The purity and concentration of the RNA samples were assessed on a NanoPhotometer1 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The integrity of the RNA samples was assessed with the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Library preparation and RNA sequencing As input material, a total of 3 μg RNA from per sample was used for the RNA sample preparations. The transcriptome library for sequencing was constructed using the NEBNext1 Ultra™ RNA Library Prep Kit for Illumina1 (NEB, USA) according to the manufacturer's instructions. Using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) following the manufacturer's 3 / 14 recommendations, the index-coded samples were clustered on a cBot Cluster Generation System. After cluster generation, the library preparations were sequenced using an Illumina HiSeq 2000 platform and 100 bp paired-end reads were generated; this was followed by FASTQ file generation and the failed reads elimination by CASAVA ver.1.8.2 (Illumina). Quality control for paired-end reads Using CASAVA ver.1.8.2 (Illumina), the sequencing-derived raw images were transformed into raw reads by base calling. The raw reads were cleaned by our self-written Perl scripts. In this step, low quality reads (more than half of the reads with a phred base quality score of less than 5) were removed to obtain clean reads. In addition, the description statistics for the clean data, such as Q20, Q30, and GC-content were calculated for high-quality downstream analysis. All downstream analyses were based on the clean reads. Reads mapping to the reference genome The chicken reference genome UMD 4.1 and model annotation files were downloaded directly from the website (ftp://ftp.ensembl.org/pub/release-75/fasta/gallus_gallus/dna/) to be utilized for the assembly. The index of the reference genome was built using Bowtie v2.2.3 and pairedend clean reads for each individual chicken were aligned to the reference genome using TopHat v2.0.12. Additionally, HTSeq v0.6.1 was used to count the reads numbers mapped to each gene. Identification of DEGs Differential expression analysis of different groups (the high and low groups with phenotypic values for percentage of polyunsaturated fatty acids) was performed using the DESeq R package (1.10.1). Using a generalized linear model based on the negative binomial distribution, DESeq2 provides statistical counts for determining DEGs in digital gene expression data. Furthermore, the Hochberg and Benjamini method was used to adjust the p-values to control for the false discovery rate. The fold changes (in log 2 scale), p-values, and q-values (corrected pvalues) of the DEGs were acquired in the output files from DESeq2. An adjusted p-value of 0.05 was assigned as the threshold for significant differential expression. GO and gene functional analysis of DEGs GO and pathway enrichment analysis of DEGs was implemented in the GOseq R package (version 2.12), in which gene length bias was corrected. GO terms and KEGG pathways (http:// www.genome.jp/kegg/) with p-value less than 0.01 were considered significantly enriched among the differential expressed genes. Validation of RNA-Seq results with qRT-PCR To validate the sequencing results, qRT-PCR was carried out on 10 randomly selected DEGs. Using identical samples with RNA-seq, total RNA was converted to cDNA with Superscript III (Invitrogen, CA, USA) following the manufacturer's instructions and was used as PCR templates. Primers were designed via Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/primer3/input. htm) and are shown in S3 File. qRT-PCR was carried out in triplicate with the LightCycler1 480 SYBR Green I Master Kit (Roche Applied Science, Penzberg, Germany) in a 15 μ L reaction on a LightCycler1 480 (Roche), using the following program: 95 ÊC for 8 min; 45 cycles of 95 ÊC for 10 s, 60 ÊC for 15 s, and 72 ÊC for 10 s; 72 ÊC for 10 min. The mRNA levels of the DEGs were normalized by the housekeeping genes GAPDH, β-actin and 18s rRNA, in the 4 / 14 corresponding samples. The relative gene expression values were calculated using the 2−ΔΔCt method. Finally, the correlations between RNA-Seq for 10 genes and the mRNA expression level from qRT-PCR were estimated using R (V3.2). Results Fatty acids profiles determined by GC-MS A typical chromatogram for the analysis of the 37-component FAMEs reference standard was shown in S4 File. Interestingly, Compared with the AA broiler, the Huangshan Black Chicken displayed higher polyunsaturated fatty acids percentage and intramuscular fat ratio (as shown in S1 File). Meanwhile, all the target compounds can be baseline separated by HP-88GC column with excellent peak shape as indicated in the chromatogram. Meanwhile, the fatty acids profiles of different chicken thigh muscle tissues in the comparison group of FAH and FAL was shown in Table 1. The analysis of the fatty acid profile showed a higher level of PUFAs in the group FAH. Particularly, contents of C22:6n-3 was higher (p<0.01), while that of C18:2n6, C20:3n-3 and C20:4n-6 were relatively higher (p<0.01) in comparison to the FAL group. RNA sequencing of thigh muscle tissue We acquired a total of 358.88 million clean reads with an average of 59.81 million (range, 57.31 to 62.63 million) for each sample (Table 2). Alignment of the sequence reads against the 5 / 14 chicken reference genome UMD 4.1 yielded 69.01±74.15% of uniquely aligned reads across the six samples, of which 71.3±80.0% fell in annotated exons, 5.2±8.2% were located in introns, and the remaining 14.8±20.5% were assigned to intergenic regions (S5 File). The data sets analyzed are available in the NCBI GenBank (https://www.ncbi.nlm.nih.gov/genbank) and the BioProject ID is PRJNA412788. Additionally, 18,964 mRNA transcripts were detected. Furthermore, the correlation coefficient (R2) between the six individuals within the FAH and FAL groups was calculated based on the FRPM value of each sample and was shown to be 0.946 and 0.969, respectively, indicating that the similarity of the three biological samples within each group was sufficiently high (S6 File). Top genes expressed in the thigh muscle tissue The differential gene expression profile between FAH and FAL was examined using the RPKM method. In total, 274 DEGs were detected significantly different between the FAH and FAL groups. Of these, 43 genes were up regulated while 231 genes were down regulated. A volcano plot of the two comparison groups that are differentially expressed and illustrate distinct transcriptional profiles is displayed in Fig 1. Moreover, using integrated analysis of RNA-Seq and gene function, the top 20 genes with the highest absolute value of expression in the thigh muscle tissue between FAH and FAL are shown in Table 3. Strikingly, the fat associated genes FADS2, OGN, and CD2 accounted for a significant proportion. Validation of differentially expressed genes To validate the RNA-Seq results, 10 random DEGs including FADS2, ABI3BP, FBLN1, DCN, LUM, FRZB, OGN, CA10, EDA2R, and ZIC4 were selected for qRT-PCR analysis. The correlations between the mRNA expression level from qRT-PCR and RNA-Seq were all consistent (Fig 2), validating the reproducibility and repeatability of gene expression data in this study. Gene ontology enrichment and pathway analysis To further study the functional associations of the detected 20 DEGs, gene ontology (GO) analysis was performed. Several important GO categories were enriched (p < 0.01), including GO processes related to synthesis, transport, and metabolic processing of lipids. For fatty acids traits, the important pathways identified were `fatty acids metabolic process,'`acyl-CoA desaturase activity',`lipid biosynthetic process,' and `unsaturated fatty acids biosynthetic process,º which also involved several candidate genes. The detailed gene function and pathway analysis of genes are shown in Table 4. 6 / 14 Fig 1. Volcano plot displaying DEGs within two different comparison groups. Note: the y-axis shows the mean expression value of log10(q-value), and the x-axis displays the log2fold change value. The blue dots represent the transcripts did not reach statistical significance (q > 0.05); the red dots represent whose expression levels were significantly different (q < 0.05). Candidate genes Integrated analysis of DEGs, GO, and pathway results, QTL databases and gene function allows us to suggest FADS2, DCN, FRZB, OGN, PRKAG3, LHFP, CHCHD10, CYTL1, FBLN5, and ADGRD1 as the 10 promising candidate genes for affecting fatty acids concentration. In addition, 10 genes (ABI3BP, NAV3, LUM, CD2, BMPER, LAMB1, FAP, SOD3, FSTL1, and CD38) were also revealed to be associated with fatty acids traits by the significant expression levels of DEGs and the QTL databases. The details of the above candidate genes identified are listed in Table 2. Discussion Polyunsaturated fatty acids play vital roles in multiple physiological processes, they participate in structural functions as major components of biomembranes [ 17 ], metabolic energy production, ligands for transcription factors, and messengers in cellular pathways [ 18 ]. Additionally, they can regulate the metabolism of lipids and promote the growth and development of animals. For poultry, the content of polyunsaturated fatty acids in adipose tissues directly affects the flavor of the meat. In our study, the Huangshan Black Chicken displayed higher polyunsaturated fatty acids percentage by comparing performance traits with AA broiler. Nonetheless, 7 / 14 OGN 12 CYTL1 4 CHCHD10 15 FRZB ADGRD1 FBLN5 CD2 LHFP BMPER LAMB1 FAP LUM SOD3 Fatty acid desaturase 2 2.13E-05 Regulated unsaturation of fatty acids, fatty acid beta-oxidation) Decorin 5.35E-08 Collagen fibril assembly, tumor suppression ABI family member 3 binding protein 7.36E-08 Heparin binding and glycosaminoglycan binding Neuron navigator 3 7.83E-05 ATPases associated with a variety of cellular activities Protein Kinase AMP-Activated Non- 1.33E-05 cell proliferation, cell differentiation in immune responses Catalytic Subunit Gamma 3 oglycin 6.92E-07 Ectopic bone formation and osteoblast differentiation Cytokine like 1 1.45E-07 Receptor binding bear the CD34 surface marker Coiled-coil-helix-coiled-coil-helix 2.17E-06 Cristae morphology maintenance, oxidative phosphorylation, and domain containing 10 Mitochondrial protein import Frizzled-related protein 3.12E-06 Involved in the regulation of bone development Adhesion G protein-coupled receptor D1 4.27E-04 Transduced extracellular signals through G proteins Fibulin 5 1.01E-04 Promoted adhesion of endothelial cells CD2 molecule 3.12E-02 Immune recognition with LFA3 on antigen presenting cells Lipoma HMGIC fusion partner 1.04E-04 A gene associated with translocation-associated lipoma BMP binding endothelial regulator 7.97E-03 Inhibited osteoblast differentiation of the chondrogenic cells Laminin subunit beta 1 9.75E-05 Cell adhesion, differentiation, signaling and metastasis Fibroblast activation protein 1.43E-03 Fibroblast growth or epithelial-mesenchymal interactions Lumican 9.72E-07 Fibril organization, circumferential growth and tissue repair. Superoxide dismutase 3 9.78E-03 Antioxidant enzyme that catalyzed the conversion of superoxide radicals into hydrogen peroxide and oxygen Follistatin like 1 1.12E-03 An autoantigen associated with rheumatoid arthritis CD38 molecule 4.12E-03 An intracellular calcium ion mobilizing messenger Fig 2. Correlations of mRNA expression level of 10 randomly DEGs between high and low polyunsaturated fatty acids percentage using RNA-Seq and qRT-PCR. Note: the x- and y-axis correspond to the log2 (ratio of FAH/FAL) measured by RNA-Seq and qRT-PCR, respectively. PLOS ONE | https://doi.org/10.1371/journal.pone.0195132 8 / 14 the precise mechanisms of Huangshan Black Chicken contributing to fatty acids composition remain unclear. Compared with traditional cDNA microarray technologies, RNA-Seq has many advantages, such as greater dynamic range, reduced bias, lower frequency of false-positive signals, and higher reproducibility [19±20]. Moreover, the correlations between RNA-Seq and the mRNA expression level from qRT-PCR were relatively high [ 10,20 ]. Three biological replicates were used for each condition to ensure broader application in our study design, the more replicates, the greater the detection power [21]. By comparative analysis and imperative validation, we detected 274 DEGs between Huangshan black chickens with extremely high and low phenotypic values for polyunsaturated fatty acids percentage. Among them, 10 genes were identified to be located within QTL areas [22] that were affirmed to have large genetic effects on fatty acids composition, including FADS2, DCN, FRZB, OGN, PRKAG3, LHFP, CHCHD10, CYTL1, FBLN5, and ADGRD1. Fatty Acid Desaturase 2 (FADS2) is one of the key limiting enzymes in the lipid metabolic pathway, which converts linoleate and alpha-linolenate into PUFAs [ 23 ]. The SNPs in the 3' untranslated regions of the FADS2 gene have significant genetic effects on the composition of fatty acids in gene expression activity in milk and blood [24±26]. Zhu et al [ 27 ] suggested that the SNPs of the FADS2 gene affect the content of essential fatty acids in muscle, and played a role in the early-stage growth rate of chickens. To investigate changes in the muscle transcriptome by increased consumption of omega-6 and omega-3 fatty acids in the pig gluteus medius muscle, Ogøuszka et al [ 28 ] showed that FADS2 may be an important gene for fatty acids metabolism. By liver transcriptome induced by a diet enriched with omega-6 and omega-3 fatty acids, Szostak et al [ 29 ] showed that FADS2 is responsible for coding enzymes delta6-desaturase. FADS2 was reported to negatively regulate fat synthesis. This evidence is consistent with the results in this study. Considering the performance traits, we supposed that FADS2 act mainly in the omega-6 metabolic pathway. Integrated analysis indicated that FADS2 is one of the most important candidate genes for polyunsaturated fatty acids percentage in chicken. As a small leucine-rich proteoglycan, Decorin (DCN) distributed in the extracellular matrix and reported to be associated with the cell membranes in tissues [30±31]. In the chicken, the existence of the core protein influencing two glycosaminoglycan chains has also been described [ 32 ]. DCN always acts as a ligand for the receptor tyrosine kinases including the insulin-like growth factor receptor [ 33 ] and the hepatocyte growth factor receptor [ 34 ]. 9 / 14 Furthermore, the expression of the DCN gene can promote the differentiation of myoblasts, the formation of muscle fibers, and the regeneration of muscle [ 35 ]. Similarly, our study revealed that DCN was near to the peak positions of four QTLs for fat traits. Frizzled motif associated with bone development (FRZB) is a protein-coding gene. As a member of the Wnt signaling pathway, FRZB can influence normal cellular processes through activating frizzled receptors [ 36 ]. In addition, FRZB is a competitor for the cell-surface G-protein receptor Frizzled affecting skeletal development in the embryo and fetus [ 37 ]. Bennett et al [ 38 ] have shown that FRZB1, FRZB2, and FRZB5 are expressed in preadipocytes and mediate the repressive effects of Wnt on adipogenesis. However, Soukas et al [ 39 ] have observed that Frizzled4 is expressed in primary adipocytes but not in 3T3-L1 cells. Wang et al [ 40 ] indicated that FRZB expression has a positive association with fat deposition and a negative association with muscle growth and inferred that FRZB may be a major candidate gene for growth traits in pigs. Combined with the function of FRZB in the metabolism, our data implied that FRZB might be involved in fat metabolism through the Wnt signaling pathway. As a regulatory subunit of the AMPK, protein kinase AMP-activated non-catalytic subunit gamma 3 (PRKAG3) takes part in regulating cellular energy homeostasis in a wide variety of tissues and cells [41±42]. By inactivating ACC oxidase (ACC) and HMG-CoA reductase (HMGCR), PRKAG3 was reported to be associated with meat quality [ 43 ], which was consistent with the previous study [ 44 ]. Likewise, nucleotide variants of PRKAG3 were able to produce significant effects on fat traits, such as final pH, meat color, and waterholding capacity in pigs [ 45 ]. Hence, PRKAG3 was considered as a major gene affecting fat traits. In addition, comprehensive analysis of differential expression patterns, biological functions, and QTL, revealed that six other genes, namely OGN, LHFP, CHCHD10, CYTL1, FBLN5, and ADGRD1, were also associated with fat acids composition traits to some extent. GO and IPA analysis showed that both OGN and ADGRD1 are involved in accumulation of adipocyte, apoptosis, and differentiation. Osteoglycin (OGN) is involved in matrix assembly, cellular growth, and migration [ 46 ]. It is believed that OGN may be a vital humoral bone anabolic factor produced by muscle tissues [ 47 ]. Donati et al indicated that OGN regulated lipid differentiation through the Wnt/β-catenin signaling pathway [ 48 ]. Lipoma HMGIC Fusion Partner (LHFP) belongs to a subset of the super family of tetraspan transmembrane protein encoding genes. Mutations in the LHFP gene result in translocation-associated lipoma. CHCHD10 belongs to a family of mitochondrial proteins and plays a role in the mitochondrial DNA stability maintenance and mitochondrial cristae morphology [ 49 ]. CYTL1 is a cytokine-like protein specifically expressed in the bone marrow and mainly takes part in DNA repair, metabolism, and cell migration. FBLN5, a matricellular protein, plays critical roles in cell proliferation [ 50 ], vascular remodeling, and smooth muscle development [ 51 ]. Mice lacking in FBLN5 exhibit systemic elastic fiber defects [ 52 ]. As a member of the adhesion-GPCR family of receptors, ADGRD1 is involved in the control of fat and adipocyte differentiation [ 53 ]. No previous studies have linked CHCHD10 or CYTL1 with lipid differentiation and further study of these genes seems to be warranted. A total of 274 genes were found to differ significantly between FAL and FAH, while some of the genes with known functions [ 54 ], e.g. FASN, FABP4, and SCD1, for fat acids composition and metabolism did not differ in the present study. It is likely that these genes with strong effects have been fixed by long-term genetic selection and no obvious differences have been observed between the comparison groups. It is also likely that different chicken populations were tested in previous studies. 10 / 14 Conclusions This study provided a global view of the complexity of the transcriptome of thigh muscle tissues of six Huangshan Black Chickens, and revealed 274 DEGs between Huangshan Black Chickens with extremely high and low phenotypic values of polyunsaturated fatty acids percentage. Integrated analysis of differential gene expression, QTL data and biological functions indicated that ten genes, including FADS2, DCN, FRZB, OGN, LUM, LHFP, CHCHD10, CYTL1, FBLN5, and ADGRD1 represent the most promising candidates affecting meat fatty acids percentage of chicken. Supporting information S1 File. Performance traits of AA broilers and Huangshan Black chickens. (PDF) S2 File. Analysis of FAMEs in thigh muscle of Huangshan Black chickens. (PDF) S3 File. PCR primers for qRT-PCR validation of 10 DEGs between the two different comparison groups. (PDF) S4 File. GC-MS analysis of 37-component FAMEs standard mixture. (PDF) S5 File. The basic statistics for RNA-Seq reads generated from thigh muscle tissues with high and low fat acids percentage. (PDF) S6 File. Correlation between biological replicates within three samples. The x- and y-axis correspond to the FPKM value of each sample, respectively. The correlation coefficient (R2) between two individuals within each group was calculated based on the FPKM value of each individual. (PDF) Author Contributions Conceptualization: Shaohua Yang, Zhaoyuan Shi, Wei Wang, Guoqing Liu. Data curation: Shaohua Yang, Ying Wang, Lulu Wang, Guoqing Liu. Formal analysis: Ying Wang, Lulu Wang, Zhaoyuan Shi, Xiaoqian Ou. Funding acquisition: Guoqing Liu. Investigation: Xiaoqian Ou, Hao Hu. Methodology: Lulu Wang, Wei Wang. Project administration: Shaohua Yang, Wei Wang. Resources: Lulu Wang, Dan Wu. Software: Zhaoyuan Shi. Supervision: Wei Wang, Guoqing Liu. Validation: Zhaoyuan Shi, Xiaoqian Ou, Jia Yuan. 11 / 14 Visualization: Xiaoqian Ou, Xinmiao Zhang, Guoqing Liu. Writing ± original draft: Shaohua Yang, Ying Wang, Fuhu Cao, Guoqing Liu. Writing ± review & editing: Shaohua Yang, Guoqing Liu. 12 / 14 22. Wardecka B, Olszewski R, Jaszczak K, Zieba G, Pierzchaøa M, Wicińska K. Relationship between microsatellite marker alleles on chromosomes 1±5 originating from the Rhode Island Red and Greenlegged Partrigenous breeds and egg production and quality traits in F(2) mapping population. J Appl Genet. 2002; 43: 319±29. PMID: 12177521 13 / 14 1. Fan B , Du ZQ , Gorbach DM , Rothschild MF . Development and application of high-density snp arrays in genomic studies of domestic animals . Asian-Australas J Anim Sci . 2010 ; 23 : 833 ± 847 . 2. Groenen MA , Megens HJ , Zare Y , Warren WC , Hillier LW , Crooijmans RP , et al. The development and characterization of a 60 K SNP chip for chicken . BMC Genomics . 2011 ; 12 : 274 . https://doi.org/10. 1186/ 1471 -2164-12-274 PMID: 21627800 3. Liu W , Li D , Liu J , Chen S , Qu L , Zheng J , et al. A genome-wide SNP scan reveals novel loci for egg production and quality traits in white leghorn and brown-egg dwarf layers . Plos One . 2011 ; 6: e28600 . https://doi.org/10.1371/journal.pone. 0028600 PMID: 22174844 4. Gu X , Feng C , Ma L , Song C , Wang Y , Da Y , et al. Genome-wide association study of body weight in chicken F2 resource population . PLoS One . 2011 ; 6: e21872 . https://doi.org/10.1371/journal.pone. 0021872 PMID: 21779344 5. Xie L , Luo C , Zhang C , Zhang R , Tang J , Nie Q , et al. Genome-wide association study identified a narrow chromosome 1 region associated with chicken growth traits . Plos One . 2012 ; 7: e30910 . https:// doi.org/10.1371/journal.pone. 0030910 PMID: 22359555 6. Liu R , Sun Y , Zhao G , Wang F , Wu D , Zheng M , et al. Genome-wide association study identifies loci and candidate genes for body composition and meat quality traits in beijing-you chickens . Plos One . 2013 ; 8: e61172 . https://doi.org/10.1371/journal.pone. 0061172 PMID: 23637794 7. Zhang T , Fan QC , Wang JY , Zhang GX , Gu YP , Tang Y. Genome-wide association study of meat quality traits in chicken . Genet Mol Res . 2015 ; 14 : 10452 ± 10460 . https://doi.org/10.4238/ 2015 .September. 8 .6 PMID: 26400276 8. Fife MS , Howell JS , Salmon N , Hocking PM , Van Diemen PM , Jones MA , et al. Genome-wide SNP analysis identifies major QTL for Salmonella colonization in the chicken . Anim Genet 2011 ; 42 : 134 ± 140 . https://doi.org/10.1111/j.1365- 2052 . 2010 . 02090 .x PMID: 20579012 9. Sun Y , Li Q , Hu Y , Sun Y , Liu R , Zheng M , et al. Genome wide association study of immune traits in chicken F2 resource population . J Anim Breed Genet 2016 ; 133 : 197 ± 206 . https://doi.org/10.1111/jbg. 12186 PMID: 26853217 10. Miao X , Luo Q , Qin X , Guo Y , Zhao H . Genome-wide mRNA-seq profiling reveals predominant downregulation of lipid metabolic processes in adipose tissues of Small Tail Han than Dorset sheep . Biochem Biophys Res Commun . 2015 ; 467 : 413 ± 420 . https://doi.org/10.1016/j.bbrc. 2015 . 09 .129 PMID: 26420224 11. Truong AD , Hong YH , Lillehoj HS . High-throughput sequencing reveals differing immune responses in the intestinal mucosa of two inbred lines afflicted with necrotic enteritis . Vet Immunol Immunopathol . 2014 ; 166 : 116 ± 124 . 12. Blech-Hermoni Y , Ladd AN . Identification of transcripts regulated by CUG-BP, Elav-like family member 1 (CELF1) in primary embryonic cardiomyocytes by RNA-seq . Genom Data . 2015 ; 6 : 74 ± 76 . https:// doi.org/10.1016/j.gdata. 2015 . 08 .014 PMID: 26366374 13. Liu L , Fan Y , Zhang Z , Yang C , Geng T , Gong D , et al. Analysis of gene expression and regulation implicates C2H9orf152 has an important role in calcium metabolism and chicken reproduction . Anim Reprod Sci . 2017 ; 176 : 1± 10 . https://doi.org/10.1016/j.anireprosci. 2016 . 11 .002 PMID: 27889102 14. Wang J , Zhao C , Li J , Feng Y , Gong Y . Transcriptome analysis of the potential roles of FOXL2 in chicken pre-hierarchical and pre-ovulatory granulosa cells . Comp Biochem Physiol Part D Genomics Proteomics 2017 ; 21 : 56 ± 66 . https://doi.org/10.1016/j.cbd. 2016 . 12 .003 PMID: 28076754 15. Ren LJ , Huang H , Xiao AH , Lian M , Jin LJ , Ji XJ . Enhanced docosahexaenoic acid production by reinforcing acetyl-CoA and NADPH supply in Schizochytrium sp . HX-308. Bioprocess and Biosyst Eng . 2009 ; 32 : 837 ± 843 . 16. Chomczynski P , Sacchi N. Single-step method of RNA isolation by acid guanidinium thiocyanate-phenol-chloroform extraction . Anal Biochem . 1987 ; 162 : 156 ± 159 . https://doi.org/10.1006/abio. 1987 .9999 PMID: 2440339 17. Kartikasari LR , Hughes RJ , Geier MS , Makrides M , Gibson RA . Dietary alpha-linolenic acid enhances omega-3 long chain polyunsaturated fatty acid levels in chicken tissues . Prostaglandins Leukotrienes & Essential FattyAcids 2012 ; 87 : 103 ± 109 . 18. Twining CW , Lawrence P , Winkler DW , Flecker AS , Brenna JT . Conversion efficiency of alpha linolenic acid to omega-3 highly unsaturated fatty acids in aerial insectivore chicks . The Journal of Experimental Biology 2017 ; https://doi.org/10.1242/jeb.165373 PMID: 29217628 19. Cui X , Hou Y , Yang S , Xie Y , Zhang S , Zhang Y , et al. Transcriptional profiling of mammary gland in Holstein cows with extremely different milk protein and fat percentage using RNA sequencing . BMC Genomics . 2014 ; 15 : 226 . https://doi.org/10.1186/ 1471 -2164-15-226 PMID: 24655368 20. Li C , Cai W , Zhou C , Yin H , Zhang Z , Loor JJ , et al. RNA-Seq reveals 10 novel promising candidate genes affecting milk protein concentration in the Chinese Holstein population . Sci Rep 2016 6 : 26813 . https://doi.org/10.1038/srep26813 PMID: 27254118 21. Rapaport F , Khanin R , Liang Y , Pirun M , Krek A , Zumbo P , et al. Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data . Genome Biol 2013 14: R95 . https://doi.org/ 10.1186/gb-2013 -14-9-r95 PMID: 24020486 23. Xu C , Wang M , Zhu X , Zhang Y , Xia H , Liu X , et al. Effects of SNPs in the 3' Untranslated Regions of FADS2 on the Composition of Fatty Acids in Milk of Chinese Holstein . Sci Sin Agric . 2016 ; 49 : 2194 ± 2202 . 24. Lopes-Marques M , OzoÂrio R , Amaral R , Tocher DR , Monroig O , Castro LF . Molecular and functional characterization of a fads2 orthologue in the Amazonian teleost, Arapaima gigas . Comp Biochem Physiol B Biochem Mol Biol . 2017 ; 203 : 84 ± 91 . https://doi.org/10.1016/j.cbpb. 2016 . 09 .007 PMID: 27693628 25. Solakivi T , Kunnas T , Jaakkola O , Renko J , LehtimaÈki T , Nikkari ST . ( 2013 ) Delta-6-desaturase gene polymorphism is associated with lipoprotein oxidation in vitro . Lipids Health Dis . 2013 ; 12 : 1 ± 6 . 26. Song Z , Cao H , Qin L , Jiang Y. A case-control study between gene polymorphisms of polyunsaturated fatty acid metabolic rate-limiting enzymes and acute coronary syndrome in Chinese Han population . BioMed Res Int . 2013 : 928178 . https://doi.org/10.1155/ 2013 /928178 PMID: 23555103 27. Zhu SK , Tian YD , Zhang S , Chen QX , Wang QY , Han RL , et al. Adjacent snps in the transcriptional regulatory region of the fads2 gene associated with fatty acid and growth traits in chickens . Genet Molecular Res . 2014 ; 13 : 3329 ± 3336 . 28. Ogøuszka M , Szostak A , Te Pas MFW , Poøawska E , Urbański P , Blicharski T , et al. A porcine gluteus medius muscle genome-wide transcriptome analysis: dietary effects of omega-6 and omega-3 fatty acids on biological mechanisms . Genes Nutr . 2017 ; 12 : 4. https://doi.org/10.1186/s12263-017 -0552-8 PMID: 28163789 29. Szostak A , Ogøuszka M , Te Pas MF , Poøawska E , Urbański P , Juszczuk-Kubiak E , et al. Effect of a diet enriched with omega-6 and omega-3 fatty acids on the pig liver transcriptome . Genes Nutr . 2016 ; 11 : 9. https://doi.org/10.1186/s12263-016 -0517-4 PMID: 27482299 30. Chen S , Birk DE . Focus on molecules: decorin, Exp Eye Res 2011 ; 92 : 444 ± 445 . https://doi.org/10. 1016/j.exer. 2010 . 05 .008 PMID: 20493188 31. Zagris N , Gilipathi K , Soulintzi N , Konstantopoulos K. Decorin developmental expression and function in the early avian embryo . Int J Dev Biol . 2011 ; 55 : 633 ± 639 . https://doi.org/10.1387/ijdb.113321nz PMID: 21948712 32. Lorda-Diez CI , GarcÂõa-Porrero JA , Hurle JM , Montero JA . Decorin gene expression in the differentiation of the skeletal connective tissues of the developing limb . Gene Expr Patterns . 2014 ; 15 : 52 ± 60 . https:// doi.org/10.1016/j.gep. 2014 . 04 .003 PMID: 24769017 33. SchoÈnherr E , SunderkoÈtter C , Iozzo RV , Schaefer L . Decorin a novel player in the insulin-like growth factor system . J Biol Chem . 2005 ; 280 : 15767 ± 15772 . https://doi.org/10.1074/jbc.M500451200 PMID: 15701628 34. Buraschi S , Pal N , Tylerrubinstein N , Owens RT , Neill T , Iozzo RV . Decorin antagonizes Met receptor activity and down-regulates {beta}-catenin and Myc levels . J Biol Chem . 2010 ; 285 : 42075 ± 42085 . https://doi.org/10.1074/jbc. M110.172841 PMID: 20974860 35. Lorda-Diez CI , GarcÂõa-Porrero JA , Hurle JM , Montero JA . Decorin gene expression in the differentiation of the skeletal connective tissues of the developing limb . Gene Expr Patterns 2014 ; 15 : 52 ± 60 . https:// doi.org/10.1016/j.gep. 2014 . 04 .003 PMID: 24769017 36. Moon RT , Kohn AD , De Ferrari GV , Kaykas A. WNT and beta-catenin signaling: diseases and therapies . Nat Rev Genet . 2004 ; 5 : 691 ±701 https://doi.org/10.1038/nrg1427 PMID: 15372092 37. Kongkham PN , Northcott PA , Croul SE , Smith CA , Taylor MD , Rutka JT . The SFRP family of WNT inhibitors function as novel tumor suppressor genes epigenetically silenced in medulloblastoma . Oncogene . 2010 ; 29 : 3017 ± 3024 . https://doi.org/10.1038/onc. 2010 .32 PMID: 20208569 38. Bennett CN , Ross SE , Longo KA , Bajnok L , Hemati N , Johnson KW , et al. Regulation of Wnt signaling during adipogenesis . J Biol Chem . 2002 ; 277 : 30998 ± 31004 . https://doi.org/10.1074/jbc.M204527200 PMID: 12055200 39. Soukas A , Socci ND , Saatkamp BD , Novelli S , Friedman JM . Distinct transcriptional profiles of adipogenesis in vivo and in vitro . J Biol Chem . 2001 ; 276 : 34167 ± 34174 . https://doi.org/10.1074/jbc. M104421200 PMID: 11445576 40. Wang Z , Li Q , Zhang B , Lu Y , Yang Y , Ban D , et al. Single nucleotide polymorphism scanning and expression of the FRZB gene in pig populations . Gene . 2014 ; 543 : 198 ± 203 . https://doi.org/10.1016/j. gene. 2014 . 04 .023 PMID: 24731717 41. Minokoshi Y , Alquier T , Furukawa N , Kim YB , Lee A , Xue B , et al. Amp-kinase regulates food intake by responding to hormonal and nutrient signals in the hypothalamus . Nature . 2013 ; 428 : 569 ± 574 . 42. Winder WW , Holmes BF , Rubink DS , Jensen EB , Chen M , Holloszy JO . Activation of AMP-activated protein kinase increases mitochondrial enzymes in skeletal muscle . J Appl Physiol . 2000 ; 88 : 2219 ± 2226 . https://doi.org/10.1152/jappl. 2000 . 88 .6.2219 PMID: 10846039 43. Yang Y , Xiong D , Yao L , Zhao C . An SNP in exon 11 of chicken 50-amp-activated protein kinase gamma 3 subunit gene was associated with meat water holding capacity . Anim Biotech . 2016 ; 27 : 13 ± 16 . 44. Zhao CJ , Wang CF , Deng XM , Gao Y , Wu C . Identification of single-nucleotide polymorphisms in 50 end and exons of the PRKAG3 gene in Hubbard White broiler, Leghorn layer, and three Chinese indigenous chicken breeds . J Anim Breed Genet . 2006 ; 123 : 349 ± 352 . https://doi.org/10.1111/j.1439- 0388 . 2006 . 00607 . x PMID : 16965409 45. Demeure O , Liaubet L , Riquet J , Milan D. Determination of PRKAG1 coding sequence and mapping of PRKAG1 and PRKAG2 relatively to porcine back fat thickness QTL . Anim Genet . 2004 ; 35 : 123 ± 125 . https://doi.org/10.1111/j.1365- 2052 . 2004 . 01102 . x PMID : 15025572 46. Williamson RE , Darrow KN , Giersch AB , Resendes BL , Huang M , Conrad GW , et al. Expression studies of osteoglycin/mimecan (OGN) in the cochlea and auditory phenotype of Ogn-deficient mice . Hear Res . 2008 ; 237 : 57 . https://doi.org/10.1016/j.heares. 2007 . 12 .006 PMID: 18243607 47. Tanaka K , Matsumoto E , Higashimaki Y , Katagiri T , Sugimoto T. Role of osteoglycin in the linkage between muscle and bone . J Biol Chem . 2012 ; 287 : 11616 ± 11628 . https://doi.org/10.1074/jbc. M111. 292193 PMID: 22351757 48. Donati G , Proserpio V , Lichtenberger BM , Natsug K , Sinclair R , Fujiwara H , et al. Epidermal Wnt/β-catenin signaling regulates adipocyte differentiation via secretion of adipogenic factors . Proc Nat Acad of Sci USA . 2014 ; 111 : 1501 ± 1509 . 49. An J , Shi J , He Q , Lui K , Liu Y , Huang Y , et al. CHM1/CHCHD6, a novel mitochondrial protein linked to regulation of mitofilin and mitochondrial cristae morphology . J Biol Chem . 2012 ; 287 : 7411 ± 7426 . https://doi.org/10.1074/jbc. M111.277103 PMID: 22228767 50. Schiemann WP , Blobe GC , Kalume DE , Pandey A , Lodish HF . Context-specific effects of fibulin-5 (DANCE/EVEC) on cell proliferation, motility, and invasion. Fibulin-5 is induced by transforming growth factor-beta and affects protein kinase cascades . J Biol Chem . 2012 ; 277 : 27367 ± 27377 . 51. Yanagisawa H , Davis EC , Starcher BC , Ouchi T , Yanagisawa M , Richardson JA , et al. Fibulin-5 is an elastin-binding protein essential for elastic fibre development in vivo . Nature . 2002 ; 415 : 168 ± 171 . https://doi.org/10.1038/415168a PMID: 11805834 52. Okuyama T , Shirakawa J , Yanagisawa H , Kyohara M , Yamazaki S , Tajima K , et al. Identification of the matricellular protein Fibulin-5 as a target molecule of glucokinase-mediated calcineurin/NFAT signaling in pancreatic islets . Sci Rep . 2017 ; 7 : 2364 . https://doi.org/10.1038/s41598-017 -02535-0 PMID: 28539593 53. Chan YF , Jones FC , McConnell E , Bryk J , Bunger L , Tautz D. Parallel selection mapping using artificially selected mice reveals body weight control loci . Curr Biol . 2012 ; 22 : 794 ± 800 . https://doi.org/10. 1016/j.cub. 2012 . 03 .011 PMID: 22445301 54. Li WJ , Zhao GP , Chen JL , Zheng MQ , Wen J . Influence of dietary vitamin E supplementation on meat quality traits and gene expression related to lipid metabolism in the Beijing-you chicken . Brit Poul Sci . 2009 ; 50 : 188 ± 198 .

This is a preview of a remote PDF: http://journals.plos.org/plosone/article/file?id=10.1371/journal.pone.0195132&type=printable

Shaohua Yang, Ying Wang, Lulu Wang, Zhaoyuan Shi, Xiaoqian Ou, Dan Wu, Xinmiao Zhang, Hao Hu, Jia Yuan, Wei Wang, Fuhu Cao, Guoqing Liu. RNA-Seq reveals differentially expressed genes affecting polyunsaturated fatty acids percentage in the Huangshan Black chicken population, PLOS ONE, 2018, DOI: 10.1371/journal.pone.0195132