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

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.


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.

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 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 BF 3 -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 1 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 1 Ultra™ RNA Library Prep Kit for Illumina 1 (NEB, USA) according to the manufacturer's instructions.Using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina) following the manufacturer's

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 LightCycler 1 480 SYBR Green I Master Kit (Roche Applied Science, Penzberg, Germany) in a 15 μ L reaction on a LightCycler 1 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).

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.

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

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.

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][25][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 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.

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.

Fig 1 .
Fig 1. Volcano plot displaying DEGs within two different comparison groups.Note: the y-axis shows the mean expression value of log 10 (q-value), and the x-axis displays the log 2 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

Fig 2 .
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 2 (ratio of FAH/FAL) measured by RNA-Seq and qRT-PCR, respectively.https://doi.org/10.1371/journal.pone.0195132.g002

Table 1 . Analysis of FAMEs in thigh muscle of Huangshan Black chickens in the group FAH and FAL (Mean ±SD).
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 2 ) 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).

Table 3 . Top 20 expressed genes in thigh muscle tissue with high polyunsaturated fatty acids percentage compared to low polyunsaturated fatty acids percentage. Symbol CHR No. Reads Log2 fold change Gene name q-value Gene function
72E-07 Fibril organization, circumferential growth and tissue repair.