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

Abstract

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.

About the Authors:

Shaohua Yang

Contributed equally to this work with: Shaohua Yang, Ying Wang, Lulu Wang

Roles Conceptualization, Data curation, Project administration, Writing – original draft, Writing – review & editing

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

ORCID http://orcid.org/0000-0003-4000-1748

Ying Wang

Contributed equally to this work with: Shaohua Yang, Ying Wang, Lulu Wang

Roles Data curation, Formal analysis, Writing – original draft

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Lulu Wang

Contributed equally to this work with: Shaohua Yang, Ying Wang, Lulu Wang

Roles Data curation, Formal analysis, Methodology, Resources

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Zhaoyuan Shi

Roles Conceptualization, Formal analysis, Software, Validation

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Xiaoqian Ou

Roles Formal analysis, Investigation, Validation, Visualization

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Dan Wu

Roles Resources

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Xinmiao Zhang

Roles Visualization

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Hao Hu

Roles Investigation

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Jia Yuan

Roles Validation

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Wei Wang

Roles Conceptualization, Methodology, Project administration, Supervision

Affiliation: Agricultural Products Quality and Safety Supervision and Management Bureau, Xuancheng, Anhui, P. R. China

Fuhu Cao

Roles Writing – original draft

* E-mail: [email protected] (FC); [email protected] (GL)

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

Guoqing Liu

Roles Conceptualization, Data curation, Funding acquisition, Supervision, Visualization, Writing – original draft, Writing – review & editing

* E-mail: [email protected] (FC); [email protected] (GL)

Affiliation: College of Food Science and Engineering, Hefei University of Technology, Hefei, Anhui, P. R. China

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 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 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 RNase-free 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 BF₃-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 NanoPhotometer^® 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 NEBNext^® Ultra^™ RNA Library Prep Kit for Illumina^® (NEB, USA) according to the manufacturer’s instructions. Using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) following the manufacturer's 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 paired-end 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 p-values) 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 LightCycler^® 480 SYBR Green I Master Kit (Roche Applied Science, Penzberg, Germany) in a 15 μ L reaction on a LightCycler^® 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 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:2n-6, C20:3n-3 and C20:4n-6 were relatively higher (p<0.01) in comparison to the FAL group.

[Figure omitted. See PDF.]

Table 1. Analysis of FAMEs in thigh muscle of Huangshan Black chickens in the group FAH and FAL (Mean±SD).

https://doi.org/10.1371/journal.pone.0195132.t001

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 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 (R²) 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).

[Figure omitted. See PDF.]

Table 2. Basic sequencing data statistics for each sample.

https://doi.org/10.1371/journal.pone.0195132.t002

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.

[Figure omitted. See PDF.]

Fig 1. Volcano plot displaying DEGs within two different comparison groups.

Note: the y-axis shows the mean expression value of log₁₀(q-value), and the x-axis displays the log₂fold 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).

https://doi.org/10.1371/journal.pone.0195132.g001

[Figure omitted. See PDF.]

Table 3. Top 20 expressed genes in thigh muscle tissue with high polyunsaturated fatty acids percentage compared to low polyunsaturated fatty acids percentage.

https://doi.org/10.1371/journal.pone.0195132.t003

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.

[Figure omitted. See PDF.]

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 log₂ (ratio of FAH/FAL) measured by RNA-Seq and qRT-PCR, respectively.

https://doi.org/10.1371/journal.pone.0195132.g002

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.

[Figure omitted. See PDF.]

Table 4. Summary of the GO analysis of 20 differentially expressed genes.

https://doi.org/10.1371/journal.pone.0195132.t004

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, 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 delta-6-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]. 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 water-holding 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.

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

[Figure omitted. See PDF.]

S1 File. Performance traits of AA broilers and Huangshan Black chickens.

https://doi.org/10.1371/journal.pone.0195132.s001

(PDF)

S2 File. Analysis of FAMEs in thigh muscle of Huangshan Black chickens.

https://doi.org/10.1371/journal.pone.0195132.s002

(PDF)

S3 File. PCR primers for qRT-PCR validation of 10 DEGs between the two different comparison groups.

https://doi.org/10.1371/journal.pone.0195132.s003

(PDF)

S4 File. GC-MS analysis of 37-component FAMEs standard mixture.

https://doi.org/10.1371/journal.pone.0195132.s004

(PDF)

S5 File. The basic statistics for RNA-Seq reads generated from thigh muscle tissues with high and low fat acids percentage.

https://doi.org/10.1371/journal.pone.0195132.s005

(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 (R²) between two individuals within each group was calculated based on the FPKM value of each individual.

https://doi.org/10.1371/journal.pone.0195132.s006

(PDF)

Citation: Yang S, Wang Y, Wang L, Shi Z, Ou X, Wu D, et al. (2018) RNA-Seq reveals differentially expressed genes affecting polyunsaturated fatty acids percentage in the Huangshan Black chicken population. PLoS ONE 13(4): e0195132. https://doi.org/10.1371/journal.pone.0195132

References

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. 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. 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. 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. 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. 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. 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. 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. pmid:26853217

10. Miao X, Luo Q, Qin X, Guo Y, Zhao H. Genome-wide mRNA-seq profiling reveals predominant down-regulation of lipid metabolic processes in adipose tissues of Small Tail Han than Dorset sheep. Biochem Biophys Res Commun. 2015; 467: 413–420. 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. 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. 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. 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. 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; 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. 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. 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. pmid:24020486

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 Green-legged Partrigenous breeds and egg production and quality traits in F(2) mapping population. J Appl Genet. 2002; 43: 319–29. pmid:12177521

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, Ozório R, Amaral R, Tocher DR, Monroig Ó, 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. pmid:27693628

25. Solakivi T, Kunnas T, Jaakkola O, Renko J, Lehtimä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. 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. 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. pmid:27482299

30. Chen S, Birk DE. Focus on molecules: decorin, Exp Eye Res 2011; 92: 444–445. 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. pmid:21948712

32. Lorda-Diez CI, García-Porrero JA, Hurlé 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. pmid:24769017

33. Schönherr E, Sunderkötter C, Iozzo RV, Schaefer L. Decorin a novel player in the insulin-like growth factor system. J Biol Chem. 2005; 280: 15767–15772. 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. pmid:20974860

35. Lorda-Diez CI, García-Porrero JA, Hurlé 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. 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 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. 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. 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. 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. 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. pmid:10846039

43. Yang Y, Xiong D, Yao L, Zhao C. An SNP in exon 11 of chicken 5′-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 5′ 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. 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. 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. 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. 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. 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. 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. 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. 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.