Skip to main content

Transcriptional analysis of liver from chickens with fast (meat bird), moderate (F1 layer x meat bird cross) and low (layer bird) growth potential

Abstract

Background

Divergent selection for meat and egg production in poultry has resulted in strains of birds differing widely in traits related to these products. Modern strains of meat birds can reach live weights of 2 kg in 35 d, while layer strains are now capable of producing more than 300 eggs per annum but grow slowly. In this study, RNA-Seq was used to investigate hepatic gene expression between three groups of birds with large differences in growth potential; meat bird, layer strain as well as an F1 layer x meat bird. The objective was to identify differentially expressed (DE) genes between all three strains to elucidate biological factors underpinning variations in growth performance.

Results

RNA-Seq analysis was carried out on total RNA extracted from the liver of meat bird (n = 6), F1 layer x meat bird cross (n = 6) and layer strain (n = 6), males. Differential expression of genes were considered significant at P < 0.05, and a false discovery rate of < 0.05, with any fold change considered. In total, 6278 genes were found to be DE with 5832 DE between meat birds and layers (19%), 2935 DE between meat birds and the cross (9.6%) and 493 DE between the cross and layers (1.6%). Comparisons between the three groups identified 155 significant DE genes. Gene ontology (GO) enrichment and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis of the 155 DE genes showed the FoxO signalling pathway was most enriched (P = 0.001), including genes related to cell cycle regulation and insulin signalling. Significant GO terms included ‘positive regulation of glucose import’ and ‘cellular response to oxidative stress’, which is also consistent with FoxOs regulation of glucose metabolism. There were high correlations between FoxO pathway genes and bodyweight, as well as genes related to glycolysis and bodyweight.

Conclusions

This study revealed large transcriptome differences between meat and layer birds. There was significant evidence implicating the FoxO signalling pathway (via cell cycle regulation and altered metabolism) as an active driver of growth variations in chicken. Functional analysis of the FoxO genes is required to understand how they regulate growth and egg production.

Background

Advancement in livestock production through selective breeding is perhaps best demonstrated in the poultry industry, where genetic selection, combined with advances in nutrition and improved management, have resulted in increases in meat bird growth in excess of 400% over the past 50 years [1]. Despite intense selection, there is still a significant amount of performance variation observed in commercial meat bird flocks, for both growth and feed efficiency [2, 3]. Feed costs account for ~ 70% of the variable costs of production in chicken meat enterprises [4], therefore optimizing performance is of economic importance to the producer and industry alike. Increased growth has not been achieved without unfavourable consequences as modern meat strains now predisposed to; excessive fat deposition [5, 6], increased leg deformities/lameness [7], metabolic disorders including; pulmonary hypertension, ascites, and sudden death syndrome [8, 9], and altered immune function [10], particularly when compared to slower growing lines such as layers and heritage line meat birds used in these studies.

Studies of different lines of chicken have explored physiological and/or anatomical growth constraints due to differential selection pressure. For example, comparison of heritage line meat birds unselected for growth, and commercial meat birds, demonstrates gross increases in breast muscle mass in modern meat birds [11]. A major difference was identified at d 14 post hatch where breast muscle growth of the heritage line plateaued at ~ 9% of total bodyweight, while breast muscle continued to increase in the commercial strain (at 14% of total bodyweight at d 14 to ~ 18% by d 28). Conversely, organs such as the heart, lungs and digestive system [1, 11, 12] have been shown to decrease as a percentage of bodyweight compared to heritage strains.

Similarly, comparative studies of strains allow for identification of physiological constraints. Experimental models of meat birds identified differential fatty acid metabolism in birds selected for either high or low abdominal fat [13, 14], or very low density lipoprotein (VLDL) plasma concentrations [15]. Comparisons of these lines, regardless of nutritional status, shows that total plasma lipids and lipoprotein levels are higher in the fat line, suggesting a higher rate of hepatic lipogenesis in fat-line birds [16]. Transcriptional analyses of genetically lean and fat chickens [17, 18] as well as juvenile and mature laying hens [19] also reveals differential expression (DE) and regulation of lipogenic genes. Additionally, fat-line birds have been shown to have significant activation of the early steps of insulin signalling at 9 wks post hatch, which may partially account for the increased lipogenesis in the liver [20]. Comparisons of domestic meat birds with the ancestral red jungle fowl, identified an intestinal glucose uptake ‘surge’ by means of increased brush border glucose transporter activity in meat birds at 2 wks post hatch, not seen in the red jungle fowl [21]. The general finding was that the meat birds had decreased glucose transporter activity (with the exception of wk. 2), but had higher glucose transporter capacity, due to an overall increase in small intestinal mass. Furthermore, modern meat birds have been shown to be less immunologically responsive to immune challenges in comparison to heritage lines [10] and more recent studies have associated gut microbes with improved feed conversion ratio (FCR; feed intake per unit of bodyweight gain) in meat birds [22, 23]. These examples are far from exhaustive, but highlight the value of comparing phenotypically different breeds and/or lines with different trait selection histories, to identify key biological pathways involved.

Meat and layer strain chickens have undergone differential genetic selection, with meat strains for high carcass yield and feed efficiency (reduced feed conversion ratio), and layers for high egg production and also reduced feed conversion ratio [24], but also lower bodyweight. Selection pressure on different traits has resulted in meat and layer stains with vastly divergent growth potential, with the bodyweights of meat birds being five times that of layers by d 42 post hatch [25]. This divergent growth rate makes meat birds and layers an excellent phenotypic model to study the underlying biological mechanisms contributing to growth and performance (i.e. FCR). However, negative consequences associated with rapid growth rates of meat birds can complicate comparisons, particularly metabolic disturbances, which in themselves may be associated with dramatic shifts in gene expression. In order to bridge the phenotypic gap between meat birds and layers, we used an intermediate growth phenotype for comparison by crossing layer ISA Brown roosters with a line of commercial meat bird breeder hens, producing an F1 layer x meat bird cross.

RNA-Sequencing (RNA-Seq) has recently been used to explore gene expression in livers of juvenile and laying hens to assess differences in the transcriptome at the different developmental stages [19] and also to study differences in the transcriptome of abdominal fat between genetically lean and fat strains of meat birds [18]. In the current study, we hypothesised that genes driving growth and performance variation in poultry could be discovered in genes DE between three groups of birds with differing growth potentials. We utilised our previous differential growth phenotypes [26] to compare the liver transcriptomes of meat birds, F1 layer x meat bird crosses and layer line males at d 14 post hatch. Day 14 was selected as the primary sampling time point due to the rapid increase in growth observed in meat strains, of which the aim was to capture transcriptional changes at the beginning of this phase for comparison with slower growing birds. The objective was to identify genes and biological pathways contributing to growth and performance differences.

Methods

Birds and management

The University of Adelaide Animal Ethics committee (approval #S-2015-171) and the PIRSA Animal Ethics committee (approval # 24/15) approved all procedures. In total, 150 newly hatched male chicks were obtained from the HiChick Breeding Company Pty Ltd., Bethel, South Australia; n = 50 meat birds (commercial line), n = 50 F1 layer (Isa Brown cockerels) x meat bird (commercial line) crosses and n = 50 layers (Isa Brown). Chicks were placed in a 6 unit rearing pen (n = 25 birds/pen), separated in breed groups (n = 2 pens/breed) in a temperature controlled room at the SARDI PPPI Poultry Research Unit, Roseworthy Campus, The University of Adelaide. All birds were fed a standard commercial meat bird starter diet ad libitum with no added in-feed antimicrobials or coccidiostats, and had unrestricted access to water via nipple drinker lines. The three experimental groups of males were chosen for their growth potential: fast growing (meat bird), moderate (F1 layer x meat bird) and slow growing (layer strain). Feed conversion ratios were recorded weekly as was bodyweight and bodyweight gain. On d 14 post hatch, 36 birds (n = 12 birds/breed) were randomly selected and euthanised by cervical dislocation. Liver tissue samples were rapidly collected, snap-frozen in liquid nitrogen and stored at − 80 °C for RNA extraction and RNA-sequencing.

RNA extraction

Samples were randomly selected for total RNA extraction (n = 6/strain) using an RNeasy Plus Mini Kit (Qiagen, Hilden, Germany). Approximately 80 mg of frozen (− 80 °C) liver tissue was homogenised in 2 mL of Trizol reagent (Invitrogen, Carlsbad, CA). 1 mL aliquots of the Trizol homogenate were combined with 200 μL of chloroform and centrifuged for 15 mins at 4 °C. The upper aqueous phase (350 μL) was transferred to a gDNA eliminator spin column and centrifuged at > 8000 g (14,000 rpm) for 30 s. The flow through (300 μL) was collected and mixed with an equal volume of 70% ethanol and transferred onto RNeasy columns. The remaining collection and wash steps were performed according to the manufacturer’s instructions. RNA was eluted in 200 μL of RNA-free water. Purity and concentration was determined using UV spectrophotometry (Nanodrop 1000; Thermo Scienfic, Wilmington, DE).

RNA-Seq library construction and sequencing

RNA-Seq was carried out by the ACRF Cancer Genomics Facility, Adelaide, SA. The sample quality was analysed on an Agilent Bio-analyser (minimum RIN requirement of 7) and sequencing libraries were made using 2 μL of total RNA. PolyA mRNA isolation was performed using oligo dT beads. Libraries were prepared using KAPA Library Quantification Kits for Illumina platforms (KAPABiosystems, Massachusetts, USA). 2 × 100 nt sequencing was carried out on an Illumin HiSeq 2500 Sequencing System to generate a minimum depth of 25 million reads.

RNA-Seq analysis

Reads were returned in fastq format. FastQC and adaptor sequences were trimmed from the 3′ end of reads with Cutadapt [27]. Hisat2 [28] was used to map reads to the reference genome Galgal5.0 (ftp://ftp.ncbi.nlm.nih.gov/genomes/Gallus_gallus). Duplicate reads were then removed. Stringtie [28] was used to define the transcripts from the read mappings for each sample, and to merge the transcript definitions for all samples. Transcripts were cleaned up using in-house scripts. The number of raw read counts were calculated for each transcript and sample using the function feature Counts of the R package Rsubread [29]. Another R package, edgeR [30] was used to analyse differential gene expression using normalised counts per million transcripts (CPM) to correct for varying depth of sequence among samples. Transcript data were aggregated by gene. Genes where the maximum CPM was < 1 were removed. Gross transcriptome relationships between the three types of bird were analysed by multidimensional scaling of the CPMs.

Functional annotation analysis and statistical analysis

Functional enrichment of the DE genes between meat bird vs layer, meat bird vs cross and layer vs cross and DE between all three groups was conducted for gene ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using the web based tools in DAVID [31, 32]. Only GO terms and KEGG pathways with P < 0.05 were taken into account as significantly enriched among the DE genes. Phenotypic data, including bodyweight, bodyweight gain and liver weights (normalised and actual), were analysed by a one-way ANOVA using SPSS (IBM SPSS Statistics 22). Gene expression levels were correlated with individual bodyweights (all three groups combined) using a Pearson’s correlation in SPSS (IBM SPSS Statistics 22).

Results

Phenotypic data

Bodyweight, bodyweight gain, and liver phenotypic data are presented in Fig. 1. Starting bodyweights (mean ± SEM) at hatch were significantly different between meat bird (44.4 ± 0.4 g); cross (42.5 ± .04 g; P = 0.008) and layer birds (38.5 ± 0.4 g; P < 0.001). At d 14 post hatch, the time of RNA-Seq analysis, bodyweight was significantly different (P < 0.001) between all three groups; meat bird (560 ± 8 g); cross (311 ± 8 g) and layer birds (159 ± 2 g). Bodyweight remained different (P < 0.001) between the three groups for the remainder of the growth period to d 28, with final bodyweights (mean ± SEM) for meat birds (2102 ± 35 g); cross (1037 ± 31 g) and layers (403 ± 6 g).

Fig. 1
figure 1

Growth of liver in meat bird, cross and layer strains; a Bodyweight (g) vs days post hatch; b Bodyweight gain (g) vs weekly intervals; c Normalised liver weight (liver weight/total bodyweight) × 100 vs days post hatch; d Liver weight vs days post hatch. Error bars are ± SEM

Day 0 liver weights (mean ± SEM) did not differ between the meat birds (1.30 ± 0.04 g) and cross (1.19 ± 0.04 g; Fig. 1c). The layer livers (0.99 ± 0.02 g) were however significantly lighter than the meat bird (P < 0.001) and cross (P = 0.002) livers. From d 7 onwards, liver weights were significantly different (at P < 0.001) between all three groups for d 7,-14 and − 28. Normalised liver weights (liver weight/ bodyweight × 100; Fig. 1d) reached maximum weight in layers and crossed birds at d 7 post hatch and declined thereafter. Meat birds had a higher relative ratio and reached maximum relative liver weight later at d 14 post hatch, which was significantly different from cross (P < 0.001) and layer birds (P < 0.001). The meat birds had a more pronounced decline in relative liver weight compared to cross and layer birds between d 14-d 28. By d 28, there was no difference in normalised liver weight between any of the groups (P > 0.05).

Identification of expressed transcripts and gross transcriptional relationships

RNA-Seq generated from 27,010,839 to 52,131,987 raw 2 × 100 paired end reads per sample with the average number being: meat bird (44,346,591), cross (40,568,610) and layer (35,862,746). After filtering the low quality reads, the average number of clean reads and percent retained were; meat bird (43,887,348; 99.0%), cross (40,146,845; 99.0%) and layer (35,447,280; 98.8%). Reads were mapped to the reference genome Galgal5.0. A total of 30,586 genes were identified among the chicken liver libraries, both known and novel. After removal of genes with no or low counts in all samples (< 1 CPM), 16,968 genes remained for analysis. Gross transcriptional analysis was undertaken using multidimensional scaling to determine how similar the transcriptomes were between the three groups. The results showed separate non-overlapping clusters of type; meat bird, layer and their F1 cross, indicating that each has a distinct transcriptome (Fig. 2).

Fig. 2
figure 2

Principal component 1 vs principal component 2 analysis of meat bird (n = 6), cross (n = 6), and layer (n = 6) transcriptomes. Clusters can be seen for meat birds, layers and their F1 cross indicating gross transcriptome differences between the three groups of birds

Identification of differential gene expression

Of the 16,968 genes expressed in at least 1 sample, 6278 genes were found to be DE for at least one of the comparisons of meat birds vs layers, meat birds vs crosses or layers vs crosses (Fig. 3). Of these 6278 genes identified as DE, 5832 were DE between meat birds and layers (19%), 2935 DE between meat birds and crosses (9.6%) and 493 DE between the layers and crosses (1.6%), highlighting that the transcriptome difference was greater in the meat birds than the layer or the cross. Percentages represent; the number of DE genes/ total genes (30,586), identified in the chicken libraries. Consideration of the transcriptome difference relative to body weight increases showed that a 1.8 fold increase in body weight from layer to cross was associated with a 1.6% transcriptome difference. The 2.0 fold increase in body weight between the cross and meat birds was associated with a 9.6% transcriptome difference, while the 3.5 fold bodyweight difference between meat birds and layers was associated with a 19% transcriptome difference. Comparisons between meat birds, crossed and layer birds identified 155 genes that are DE between all three groups. Of these 155 genes, 60% were found to be progressively upregulated in the direction meat bird > cross > layer, 38.1% down regulated meat bird < cross < layer, and 1.9% did not follow any directional pattern associated with growth rate.

Fig. 3
figure 3

Venn diagram illustrating differential gene expression between meat birds (n = 6) vs layers (n = 6) (M v L); meat birds (n = 6) vs cross (n = 6) (M v C), and layers (n = 6) v cross (n = 6) (L v C)

Characterisation of the 155 DE genes

The 155 DE genes were characterised in terms of abundance and fold change. Additionally, correlations were tested between the 155 DE genes with individual bodyweight. The top 10 most abundantly expressed genes were; alpha 2-HS glycoprotein (AHSG), fibrinogen alpha chain (FGA), fibrinogen gamma chain (FGB), fibrinogen beta chain (FGG), ferratin heavy polypeptide 1 (FTH1), compliment C4 (C4), acetyl-CoA transferase 2 (ACAA2), Dihydrolipoamide S-acetyltransferase (DLAT), saccharopine dehydrogenase (SCCPDH) and one unknown (NA) (Table 1). Of these top 10 most abundantly expressed, 6 were down regulated in meat birds (AHSG, FGA, FGG, FGB, FTH1 and C4; meat bird < cross < layer), and 4 upregulated (NA, ACCA2, DLAT and SCCPDH; meat bird > cross > layer).

Table 1 The top 10 most abundantly expressed genes (mean ± SEM) presented as counts per million for male; meat (n = 6), cross (n = 6) and layer (n = 6) birds at d 14 post hatch

Fold changes were calculated using the mean CPMs for meat bird vs layer, meat bird vs cross and layer vs cross (Table 2). The largest fold change detected within the 155 DE genes was a 227-fold upregulation (CPM mean ± SEM) in meat birds (1602.90 ± 83.31) compared to layers (7.06 ± 0.89) for an uncharacterised gene (Un_24875). The second highest fold change was a 147.7 fold upregulation between the meat birds (1.12 ± 0.38) compared to layers (0.01 ± 0.01) for bacterial/permeability-increasing protein-like 3 (BPIL3). One individual bird with significantly increased expression in the meat bird group influenced the magnitude of the BPIL3 fold change. Six of the top 10 highest fold change (all genes) were novel and uncharacterised, highlighting gaps within the chicken genome. Among the top 10 characterised genes were BPIL3, LOC107055086 and LOC107057467 genes which have both been characterised as sperm-associated antigen 4 protein-like, Histamine N-methyltransferase-like (LOC771456), cyclin dependant kinase inhibitor 2B (CDKN2B), platelet glycoprotein VI-like (LOC10087809, leucine protein zipper 2 (LUZP2), butyrophilin subfamily 3 member A2-like (LOC107049070), ubiquitin C-terminal hydrolase L1 (UCHL1) and prostaglandin D2 synthase (PTGDS) (see Table 2).

Table 2 Top 10 fold changes of the 155 DE genes between meat birds, crossed and layer birds

The 155 DE genes were correlated with individual bodyweight. Of the top ten correlated (Table 3), the highest correlation was between dihydrolipoamide S-acetyltransferase (DLAT), which is the E2 component on the pyruvate dehydrogenase complex, linking glycolysis to the citric acid cycle. DLAT was also among the top 10 most abundant of the 155 DE genes. Three of the top 10 genes correlated with bodyweight are novel and uncharacterised, e.g. LOC770248 and two unknown. Other genes highly correlated with bodyweight included quiescin Q6 sulfhydryl oxidase 1 (QSOX1), receptor accessory protein 5 (REEP5), myosin VI (MYO6), transmembrane protein 246 (TMEM246), cyclin G2 (CCNG2) and WW domain binding protein (WBP2).

Table 3 Top 10 genes with highest correlation with individual bodyweight

Functional analysis of DE genes

All 6278 DE genes were analysed for GO terms and KEGG pathways using both edgeR and the web based tools in DAVID [31, 32]. There were 38 biological GO terms (GO: BP) identified for 5832 DE genes (P < 0.05) between meat bird and layer groups, 28 GO terms for 2935 DE genes between meat bird and the cross, and 19 GO terms for 493 DE genes between the layer and cross groups.

To understand the biological differences contributing to growth between the two strains and the F1 cross, we focused on the 155 DE genes among the meat birds, crossed and layer birds. For these 155 DE genes, 27 GO terms were identified (Table 4). Many of the GO terms were found to be significant due to the expression levels of FGA, FGB and FGG, which were among the most abundantly expressed genes. These three genes dominated 20 of the 27 GO terms identified, ranging from fibrinolysis, blood clot formation, fibrin clot formation, plasminogen activation, positive regulation of exocytosis, response to calcium ion and platelet aggregation. However, despite their high abundance, these genes had lower correlations with bodyweight than other DE genes (mentioned above), although still significant at P < 0.01. Of the 155 DE genes, ranked in order of correlation strength with bodyweight, FGA was 89th (r = − 0.874), FGG was 111th (r = − 0.085) and FGB was 124th (r = − 0.836). GO BP terms not largely dominated by FGA, FGB and FGG included positive regulation of glucose import, cellular response to oxidative stress and regulation of cell death. GO CC terms included chromatin and extracellular exosome. The extracellular exosome GO CC term (GO: 0070062) included 21 genes, 6 of which are in the top 10 most abundant (AHSG, FGA, FGB, FGG, FTH1 and ACCA2), 2 in the top 10 fold changes (PTGDS and UCHL1), and 3 in the top 10 correlated with individual bodyweight (QSOX, REEP5 and MYO6).

Table 4 Gene Ontology (GO) terms for the 155 DE genes

KEGG analysis of the 5832 DE genes between meat birds and layers revealed 13 pathways significantly enriched (P < 0.05). Pathways included; metabolic pathway (singular KEGG term), PPAR signalling pathway, biosynthesis of antibiotics, FoxO signalling pathway, cell cycle, drug metabolism, peroxisome, steroid biosynthesis, nicotinate and nicotinamide metabolism, glycine, serine and threonine metabolism, pentose phosphate pathway, glutathione metabolism and fatty acid metabolism. KEGG analysis of meat birds vs cross (2935 DE genes) identified 15 significantly enriched pathways, 10 pathways overlapped with those significant for meat birds vs layers including; metabolic pathway, PPAR signalling pathway, biosynthesis of antibiotics, FoxO signalling pathway, cell cycle, peroxisome, steroid biosynthesis, glycine, serine and threonine metabolism and the pentose phosphate pathway. Three pathways were significantly enriched for the layers vs cross; metabolic pathway, folate biosynthesis and FoxO signalling pathway. The metabolic and FoxO signalling pathway were the only two common pathways between the three types of bird identified (P < 0.05).

KEGG pathway analysis of the 155 genes DE between all three types of birds identified two enriched pathways at P < 0.05 (Table 5). Three genes were enriched for fructose and mannose metabolism (P = 0.024); 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 2 (PFKFB2), phosphofructokinase liver (PFKL), tissue specific transplantation antigen P35B (TSTA3). Five genes were associated with the FoxO signalling pathway (P = 0.001); cell cycle regulators; cyclin D2 (CCND2), cyclin G2 (CCNG2), cyclin-dependent kinase inhibitor 1B (CDKN1B), cyclin-dependent kinase inhibitor 2B (CDKN2B) as well as insulin receptor (INSR). Just falling out of significance at P = 0.053 was the glycolysis/gluconeogenesis pathway, involving three genes: PFKL (overlapping with fructose/mannose metabolism), which is rate limiting in glycolysis, catalysing the transformation of fructose-6-phospate to fructose-1,6-diphosphate [33]; glutamate transporter (SLC1A2), which was upregulated in the direction of meat birds (5.02 ± 0.70) crosses (2.16 ± 0.18) and layers (0.89 ± 0.14); and alcohol dehydrogenase (AKR1A1 or NADP+), which was upregulated in the direction of meat birds (278.63 ± 12.84), crosses (230.01 ± 7.51) and layers (179.24 ± 6.45). INSR occurs in the FoxO pathway, and SLC1A2 also overlaps with the GO term, GO: 0046326, positive regulation of glucose import.

Table 5 Pathways and associated genes identified as enriched by KEGG of the 155 DE genes between meat birds, crossed and layer birds

Discussion

Liver transcriptomes of males of meat birds, F1 layer x meat bird crosses and layer birds were compared to identify DE genes between all three groups. Selection of the groups were based on their fast, moderate and slow growth potential, respectively. Day 14 post hatch was selected as the primary sampling time due to the rapid increase in growth seen in meat birds from 2 to 3 wks of age compared to other strains. By sampling at this time point, it was hoped to capture transcriptional changes at the beginning of rapid growth phase to further understand the biological factors associated with the high growth rates seen in meat birds.

The results of this study revealed that selection for growth or egg laying is associated with altered transcriptomes between meat and layer birds. Bodyweight at d 14 post hatch was 1.8 fold higher for crosses vs layers, and also 2.0 fold higher for meat birds vs crosses (meat bird > cross). The difference in transcriptomes associated with birds of differing bodyweights was quite remarkable. Of the total genes analysed, 1.6% were DE between crosses and layers; 9.6% DE between meat birds and crosses; and 19% DE between meat birds and layers. The differences in gene expression observed between the meat birds, layers and their F1 cross are not all driving the increases in bird size, particularly given the confounding effect of the many metabolic disturbances modern meat birds exhibit. These include; excessive fat deposition [5, 6], increased skeletal defects [7], pulmonary hypertension, sudden death syndrome [8, 9] and altered immune function [10]. However, it is likely that the drivers of growth are represented in the DE genes, particularly those that differ between all three groups.

GO and KEGG analyses of DE genes for meat birds vs layers, meat birds vs crosses and layers vs crosses identified overlapping biological functions that were affected and may contribute to the differential growth between types of birds. Two affected KEGG pathways were identified between all three comparisons; metabolic pathway (singular KEGG term) and the FoxO signalling pathway. The Forkhead box O (FoxO) genes central to this pathway are a family of transcription factors that regulate gene expression related to cell cycle regulation, cell survival, and metabolism, including glucose and lipid metabolism [34]. KEGG analysis of the 155 genes DE between each types of birds again identified the FoxO signalling pathway, enriched at P < 0.05. The fructose and mannose metabolism pathway was also enriched, with an overlap of genes involved in glycolysis, as well as the GO term ‘positive regulation of glucose import’. Among the functions of the FoxO signalling pathway is maintenance of homeostasis, particularly in response to stress [34].

FoxOs have previously been identified as potential candidate genes for growth in chickens. A genome-wide association study using a reciprocal cross between White Recessive Rock (WRR) and Xinghua (XH) chickens, identified a 1.5 Mb region on chromosome 1 containing 5 SNPs, including a SNP 8.9 kb upstream of FoxO1 for bodyweight at 22–24 d and 70 d post hatch [35]. FoxO1 contained two SNPs in the intron region of the gene; however, these two SNPs were not significantly associated with growth traits. The authors questioned whether a regulatory mechanism was involved in the significant SNP effects associated with growth traits located up and downstream of FoxO1. The most significant SNP for average daily gain at d 42 was in a region containing gene LOC770248, which is uncharacterised. Comparatively, LOC770248 was amongst the top 10 genes correlated (r = 0.934) with individual d 14 bodyweight in the present study. The identification of LOC770248 as a potential regulator of growth traits suggests further investigation is warranted to characterise the function of the encoded protein. More recently, RNA-Seq of the breast muscle of WWR and XH chickens at 7 wks post hatch identified FoxO3 as a candidate gene (supported by siRNA analysis and association analysis) for further investigation into breast muscle growth in the chicken [36]. The significant enrichment of the FoxO signalling pathway in all comparisons in the current study strongly supports the contribution of this pathway to the growth differences between meat birds, layers and their F1 cross.

Of the 155 DE genes identified between the three types of bird, five genes associated with the FoxO signalling pathway were upregulated (meat birds > crosses > layers). These were insulin receptor (INSR), as well as genes essential for cell cycle regulation, cyclins CCND2, CCNG2 and cyclin-dependent kinase inhibitors CDKN2B. Down regulation of CDKN1B was seen in meat birds compared with crossed and layer birds (meat birds < crosses < layers). Cyclins, such as CCND2, activate cyclin-dependent protein kinases (CDKs) which form complexes to transition the cell from one cell cycle state to another [37], for example; activation of cyclin-D dependent kinases initiates progression of the cell cycle through the G1 phase [38]. CCND2 binds to several types of CDKs, with the main partners CDK4 and CDK6 [37]. We did not find CDK4 in this gene set (of the 30,586), however, found abundant levels of CDK6, although not DE expressed. CDKs are normally present in the cell in excess of their cyclin partner [37], which was the case at the RNA expression level of CDK6:CCND2 for meat birds, crosses and layers, with ratios of 1.3, 3.0 and 100.3 respectively. Interestingly, we found high DE of CDKN2B between all three groups (meat birds > crosses > layers) which inhibits the activity of CDK4 and CDK6. CDKN2B is known to weaken the binding of D-type cyclins and as well as interact with the catalytic domains of CDK4 and CDK6 as a potent inhibitor of kinase activity [39]. Meat birds were the only group that had higher levels of CDKNB2 relative to either CDK6 or CCDN2, and the ratios for CDKNB2:CCND2 and CDKNB:CDK6 decreased (meat birds > crosses > layers) in both instances. CCKNB2 was also amongst the top genes categorized by fold difference, being 28 fold higher in meat birds compared with layers.

Cyclin CCNG2 was also upregulated in meat birds compared to cross and layer birds. Unlike ‘conventional’ cyclins that promote cell cycle progression, CCNG2 upregulation in murine B cells is associated with cell cycle arrest or apoptosis in response to inhibitory stimuli, and conversely, CCND2 is down regulated during G1 phase growth arrest [40, 41]. There is limited information of CCNG2 activity in birds or in the liver for comparison. One study however compared Arbor Acres meat birds divergently selected for lean and fat lines, and identified CCNG2 with a 0.209 and 0.249 lean/fat fold change at 2 and 4 wks respectively in liver tissue, which is similar to the fold change we observed between layers/meat birds (0.296) [42]. These studies, together with CCNG2 being amongst the top 10 DE genes correlated with bodyweight in the present study, supports the differences between meat birds, crosses and layers being a result of differential cell cycle progression between the three types of birds.

Here we report that the liver (as a percentage of total bodyweight) reaches maximum size at d 14 post hatch in meat birds compared with the crossed and layer birds, where the ratio between liver and body size is lower, and the relative liver weight maximum is reached earlier, at d 7 post hatch. By d 28 there was no difference in relative liver weight (~ 3%) between any of the groups. In many plants and animals, organ scaling is controlled at the level of cell number [43]. However, for meat birds, although liver weight continued to increase at the same rate as the cross and layer birds from d 14-d 28, the expression studied, combined with the higher deceleration in relative liver weight from d 14 onwards in meat birds suggests that either; a) the total cell cycle time is increased or b) there is a decrease in growth fraction. Growth fraction has been defined as the number of cells remaining in the cell cycle vs the total organ cell number [44]. Therefore a decrease in growth fraction is likely due to fewer dividing cells as more remain in the G0 cell phase [44]. This would be supported by increased expression of CCGN2 and CDKNB2. Thus, increased growth in the meat birds compared to cross and layer birds, likely results from hypotrophy (increased cell size) rather than hypoplasia (increased cell number). Hypotrophy via cellular polyploidy in the liver is not uncommon, with polyploid cells appearing late in fetal development, coinciding with terminal differentiation [45]. Polyploidy is associated with rapid growth by facilitating an increase in cell volume without division, which may permit cells to be more metabolically active [43, 46]. Without histological analysis on hepatocytes, increased polyploidy is speculative, however, there is evidence in this study to suggest the meat birds are more metabolically active.

The insulin receptor (INSR) was down regulated in broilers compared with cross and layer birds and was amongst the 155 genes with DE between the three groups. INSR was also enriched to the FoxO signalling pathway. The INSR pathway is conserved from flies to humans, and is a key sensor of nutrient availability, playing an important role in the control of cellular proliferation, cellular size and response to nutrient availability [47, 48]. Insulin regulates not only glucose metabolism, but also lipid homeostasis by increasing lipogenesis in the case of nutrient excess. In hepatocytes, activation of FoxO promotes the expression of key gluconeogenetic and glycogenolytic enzymes in the fasted state, resulting in increased hepatic glucose production [49]. In the fed state, high insulin blocks FoxO activity through the PI3-kinase (PI3K)-Atk pathway [49]. Atk phosphorylates the FoxO protein, retaining it in the cytoplasm in its inactive state [48,49,50]. This would favour glucose uptake and glycolysis. Furthermore, FoxOs have been shown to directly regulate the insulin signalling response to nutrients in C2C12 lines [48]. Upregulated insulin mRNA levels were associated with dephosphorylation of FoxO1, conversely down regulated insulin mRNA levels were associated with phosphorylation of FoxO1 [48]. As phosphorylation of FoxO1 results in decreased activation of FoxO1, it would be anticipated that this direct effect would result in decrease gluconeogenesis and increased glucose uptake and glycolysis. A major limitation in this study is that without functional analysis of the FoxO genes themselves, we cannot determine their activation status.

The lower expression levels of INSR in meat birds compared with crosses and layers however, is consistent with increased levels of phosphofuctokinase (PFKL; upregulated meat birds > crosses > layers), glutamate transporter SLC1A2, and AK1A1, which would be expected with increased levels of glycolysis, particularly as PFKL is a rate limiting enzyme in glycolysis. PFKL catalyses the transformation of fructose-6-phospate to fructose-1,6-diphosphate [33]. Furthermore, the pyruvate dehydrogenase complex (PDC) links glycolysis to the citric acid cycle. Therefore it is significant that dihydrolipoamide S-acetyltransferase (DLAT) was in the top 10 most abundantly expressed genes, upregulated in meat birds, and showed the highest correlation with bodyweight (r = 0.968). DLAT is the E2 component of the PDC, catalysing the oxidative reaction of pyruvate (end product of glycolysis) to acetyl-CoA in the mitochondria. Interestingly, chickens have been shown not to accumulate pyruvate in the liver, so the increase in DLAT is also consistent with the conversion to, and utilisation of acetyl-CoA in the mitochondria as soon as pyruvate is formed [51].

Conclusion

In this study, we used RNA-Seq to show that the transcriptomes of meat birds, layers (and the F1 cross between them) are highly divergent, particularly between meat and layer type birds. Metabolic pathway (singular KEGG term) and the FoxO signalling pathway were identified as significantly enriched in comparisons between the three types of birds, with trends between meat, crossed and layer birds. Functional analysis of the 155 genes DE between all three strains also identified enrichment of the FoxO signalling pathway, particularly genes related to cell cycle regulation and the insulin receptor. These data suggest that differences in cell cycle regulation and glucose metabolism are associated with differences in growth rate, and provide evidence that meat birds have a higher rate of glycolysis. Functional analysis of the chicken hepatic FoxO genes and associated pathway targets warrants further investigation to determine the role of this pathway in regulating the growth of meat birds.

References

  1. Zuidhof MJ, Schneider BL, Carney VL, Korver DR, Robinson FE. Growth, efficiency, and yield of commercial broilers from 1957, 1978 and 2005. Poult Sci. 2014;93(12):2970–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Emmerson DA. Commercial approaches to genetic selection for growth and feed conversion in domestic poultry. Poult Sci. 1997;76(8):1121–5.

    Article  CAS  PubMed  Google Scholar 

  3. Tallentire CW, Leinonen I, Kyriazakis I. Breeding for efficiency in the broiler chicken: a review. Agron Sustain Dev. 2016;36:66.

    Article  Google Scholar 

  4. Aggrey SE, Karnuah AB, Sebastian B, Anthony NB. Genetic properties of feed efficiency parameters in meat type chickens. Genet Sel Evol. 2010;42(1):25.

    Article  PubMed  PubMed Central  Google Scholar 

  5. Ji B, Ernest B, Gooding JR, Das S, Saxton AM, Simon J, et al. Transcriptomic and metabolic profiling of chicken adipose tissue in response to insulin neutralization and fasting. BMC Genomics. 2012;13:441. https://0-doi-org.brum.beds.ac.uk/10.1186/1471-2164-13-441.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Fouad AM, El-Senousey HK. Nutritional factors affecting abdominal fat deposition in poultry: a review. Asian-Australas J Anim Sci. 2014;27(7):1057–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Bessei W. Welfare of broilers. Worlds Poult Sci J. 2006;62(3):455–66.

    Article  Google Scholar 

  8. Julian RT. Production and growth related disorders and other metabolic diseases of poultry. Vet J. 2005;169(3):350–69.

    Article  CAS  PubMed  Google Scholar 

  9. Olkowski AA, Wojnarowicz C, Nain S, Ling BB, Alcorn JM, Laarveld B. A study on pathogenesis of sudden death syndrome in broiler chickens. Res Vet Sci. 2007;85(1):131–40.

    Article  PubMed  Google Scholar 

  10. Cheema MA, Qureshi MA, Havenstein GB. A comparison of the immune response of a 2001 commercial broiler with a 1957 randombred broiler strain when fed representative 1957 and 2001 broiler diets. Poult Sci. 2003;82(10):1519–29.

    Article  CAS  PubMed  Google Scholar 

  11. Schmidt CJ, Persia ME, Feierstein E, Kingham B, Saylor WW. Comparison of a modern broiler line and a heritage line unselected since the 1950s. Poult Sci. 2009;88(12):2610–9.

    Article  CAS  PubMed  Google Scholar 

  12. Havenstein GB, Ferket PR, Qureshi MA. Carcass composition and yield of 1957 versus 2001 broilers when fed representative 1957 and 2001 broiler diets. Poult Sci. 2003;82(10):1509–18.

    Article  CAS  PubMed  Google Scholar 

  13. Leclercq B, Blum JC, Boyer JP. Selecting broilers for low or high abdominal fat-initial observations. Br Poult Sci. 1980;21(2):107–13.

    Article  Google Scholar 

  14. Leclercq B, Simon J. Selecting broilers for low or high abdominal fat: observations on the hens during the breeding period. Ann Zootech. 1982;31(2):161–70.

    Article  Google Scholar 

  15. Whitehead CC, Griffin HD. Development of divergent lines of lean and fat broilers using plasma very low density lipoprotein concentration as selection criterion: the first three generations. Br Poult Sci. 1984;25(4):573–82.

    Article  CAS  PubMed  Google Scholar 

  16. Hermier D, Chapman JM, Leclercq B. Plasma lipoprotein profile in fasted and refed chickens of two strains selected for high or low adiposity. J Nutr. 1984;114(6):1112–21.

    Article  CAS  PubMed  Google Scholar 

  17. Resnyk CW, Carré W, Wang X, Porter TE, Simon J, Le Bihan-Duval E, et al. Transcriptional analysis of abdominal fat in genetically fat and lean chickens reveals adipokines, lipogenic genes and a link between hemostasis and leanness. BMC Genomics. 2013;14:557.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Resnyk CW, Chen C, Huang H, Wu CH, Simon J, Le Bihan-Duval E, et al. RNA-Seq analysis of abdominal fat in genetically fat and lean chickens highlights a divergence in expression of genes controlling adiposity, hemostasis, and lipid metabolism. PLoS One. 2015;10(10):e0139549. https://0-doi-org.brum.beds.ac.uk/10.1371/journal.pone.0139549.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Li H, Wang T, Xu C, Wang D, Ren J, Li Y, et al. Transcriptome profile of liver at different physiological stages reveal potential mode for lipid metabolism in laying hens. BMC Genomics. 2015;16:763.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Dupont J, Chen J, Derouet M, Simon J, Leclercq B, Taouis M. Metabolic differences between genetically lean and fat chickens are partly attributed to the alteration of insulin signalling in the liver. J Nutr. 1999;129(11):1937–44.

    Article  CAS  PubMed  Google Scholar 

  21. Jackson S, Diamond J. Metabolic and digestive responses to artificial selection in chickens. Evolution. 1996;50(4):1638–50.

    Article  PubMed  Google Scholar 

  22. Stanley D, Denman SE, Hughes RJ, Geier MS, Crowley TM, Chen H, et al. Intestinal microbiota associated with differential feed conversion efficiency in chickens. Appl Microbiol Biotechnol. 2012;96(5):1361–9.

    Article  CAS  PubMed  Google Scholar 

  23. Stanley D, Hughes RJ, Geier MS, Moore RJ. Bacteria within the gastrointestinal tract microbiota correlated with improved growth and feed conversion: challenges presented for the identification of performance enhancing probiotic bacteria. Front Microbiol. 2016;7:187. https://0-doi-org.brum.beds.ac.uk/10.3389/fmicb.2016.00187.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Druyan S. The effects of genetic line (broilers vs. layers) on embryo development. Poult Sci. 2010;89(7):1457–67.

    Article  CAS  PubMed  Google Scholar 

  25. Zhao R, Muehlbauer E, Decuypere E, Grossman R. Effect of genotype-nutrition interaction on growth and somatotropic gene expression in the chicken. Gen Comp Endocrinol. 2004;136(1):2–11.

    Article  CAS  PubMed  Google Scholar 

  26. Willson N-L, Forder REA, Tearle RG, Nattrass GS, Hughes RJ, Hynd PI. Evaluation of fatty acid metabolism and innate immunity interactions between commercial broiler, F1 layer x broiler cross and commercial layer strains selected for different growth potentials. J Anim Sci Biotechnol. 2017;8:70. https://0-doi-org.brum.beds.ac.uk/10.1186/s40104-017-0202-4.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Martin M. Cutadapt removes adaptor sequences from high-throughput sequencing reads. EBMnet J. 2011;17(1):10.

    Google Scholar 

  28. Pertea M, Kim D, Pertea GM, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Protoc. 2016;11(9):1650–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Liao Y, Smith GK, Shi W. The subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 2013;41(10):e108.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.

    Article  CAS  PubMed  Google Scholar 

  31. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    Article  CAS  Google Scholar 

  32. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths towards the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37(1):1–13.

    Article  Google Scholar 

  33. Uyeda K. Phosphofuctokinase. Adv Enzymol Relat Areas Mol Biol. 1979;48:193–244.

    CAS  PubMed  Google Scholar 

  34. Eijkelenboom A, Burgering BMT. FOXOs: signalling integrators for homeostasis maintenance. Mol Cell Biol. 2013;14(2):83–97.

    CAS  Google Scholar 

  35. 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(2):e30910.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Chen B, Xu J, He X, Xu H, Li G, Du H, et al. A genome-wide mRNA screen and functional analysis reveal FOXO3 as a candidate gene for chicken growth. PLoS One. 2015;10(9):e0137087.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Pines J. Cyclins and cyclin-dependent kinases: a biochemical view. Biochem J. 1995;308(Pt 3):697–711.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Sherr CJ. D-type cyclins. Trends Biochem Sci. 1995;20(5):187–90.

    Article  CAS  PubMed  Google Scholar 

  39. Asghar U, Witkiewicz AK, Turner NC, Knudsen ES. The history and future of targeting cyclin-dependent kinases in cancer therapy. Nat Rev Drug Discov. 2013;14(2):130–46.

    Article  Google Scholar 

  40. Bennin DA, Arachchige Don AS, Brake T, McKenzie JL, Rosenbaum H, Oritz L, et al. Cyclin G2 associates with protein phosphatase 2A catalytic and regulatory B’ subunits in active complexes and induces nuclear aberrations and a G1/s phase cell cycle arrest. J Biol Chem. 2002;277(30):27449–67.

    Article  CAS  PubMed  Google Scholar 

  41. Horne MC, Donaldson KL, Lynn Goolsby L, Tran D, Mulheisen M, Hell JW, et al. Cyclin G2 is upregulated during growth inhibition and B cell antigen receptor-mediated cell cycle arrest. J Biol Chem. 1997;272(19):12650–61.

    Article  CAS  PubMed  Google Scholar 

  42. He Q, Wang SZ, Leng L, MNa W, Wang QC, Li H. Differentially expressed genes in the liver of lean and fat chickens. Genet Mol Res. 2014;13(4):10823–8.

    Article  CAS  PubMed  Google Scholar 

  43. Orr-Weaver TL. When bigger is better: the role of ploidy in organogenesis. Trends Genet. 2015;31(6):307–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Lui JC, Baron J. Mechanisms limiting body growth in mammals. Endocr Rev. 2011;32(3):422–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Gupta S. Hepatic polyploidy and liver growth control. Semin Cancer Biol. 2000;10(3):161–71.

    Article  CAS  PubMed  Google Scholar 

  46. Lee HO, Davidson JM, Duronio RJ. Endoreplication: polyploidy with purpose. Genes Dev. 2009;23(21):2461–77.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. Marr MT II, D'Alessio JA, Puig O, Tjian R. IRES-mediated functional coupling of transcription and translation amplifies insulin receptor feedback. Genes Dev. 2007;21(2):175–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Puig O, Tjian R. Transcriptional feedback control of insulin receptor by dFOXO/FOXO1. Genes Dev. 2005;19(20):2435–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Kousteni S. FoxO1, the transcriptional chief of staff of energy metabolism. Bone. 2012;50(2):437–43.

    Article  CAS  PubMed  Google Scholar 

  50. Puig O, Marr MT, Laure Ruhf M, Tjian R. Control of cell number by Drosophila FOXO: downstream and feedback regulation of the insulin receptor pathway. Genes Dev. 2003;17(16):2006–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Ochs RS, Harris RA. Studies on the relationship between glycolysis, lipogenesis, gluconeogenesis, and pyruvate kinase activity of rat and chicken hepatocytes. Arch of Biochem Biophys. 1978;190(1):193–201.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We would like to thank D. Schultz, M. Bowling and N. Heberle for assistance during the trial and the HiChick Breeding Company Pty Ltd. for providing the birds.

Funding

This work was financially supported by the Australian Poultry Cooperative Research Centre. The funder had no role in study design, data collection, and analysis. The funder read and approved the manuscript for publication. N-LW received a stipend from the University of Adelaide (F J Sandoz Scholarship) and the Australian Poultry Cooperative Research Centre for PhD studies.

Availability of data and materials

The datasets generated and/or analysed during the current study are available in the NCBI Sequence Read Archive (SRA) under BioProject PRJNA392351 [https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/bioproject/PRJNA392351/].

Author information

Authors and Affiliations

Authors

Contributions

N-LW, PH, RF, RH and GN designed the study. N-LW and RF were involved in performing the experiment. N-LW conducted the laboratory analysis. RT performed RNA-seq data analysis. N-LW, PH, RF, GN, RT and JLW were involved in interpreting the results. N-LW performed GO, KEGG and data analysis and wrote the manuscript, PH, RF, RT, JLW, RH and GN revised the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Nicky-Lee Willson.

Ethics declarations

Ethics approval

The experimental protocol used in this study, including animal management, housing, and slaughter procedures were by the University of Adelaide Animal Ethics Committee (approval #S-2015-171) and the PIRSA Animal Ethics committee (approval #24/15).

Competing interests

The authors declare they have no competing interests.

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Willson, NL., Forder, R.E.A., Tearle, R. et al. Transcriptional analysis of liver from chickens with fast (meat bird), moderate (F1 layer x meat bird cross) and low (layer bird) growth potential. BMC Genomics 19, 309 (2018). https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-018-4723-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/s12864-018-4723-9

Keywords