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

Genetic selection for enhanced growth rate in meat-type chickens (Gallus domesticus) is usually accompanied by excessive adiposity, which has negative impacts on both feed efficiency and carcass quality. Enhanced visceral fatness and several unique features of avian metabolism (i.e., fasting hyperglycemia and insulin insensitivity) mimic overt symptoms of obesity and related metabolic disorders in humans. Elucidation of the genetic and endocrine factors that contribute to excessive visceral fatness in chickens could also advance our understanding of human metabolic diseases. Here, RNA sequencing was used to examine differential gene expression in abdominal fat of genetically fat and lean chickens, which exhibit a 2.8-fold divergence in visceral fatness at 7 wk. Ingenuity Pathway Analysis revealed that many of 1687 differentially expressed genes are associated with hemostasis, endocrine function and metabolic syndrome in mammals. Among the highest expressed genes in abdominal fat, across both genotypes, were 25 differentially expressed genes associated with de novo synthesis and metabolism of lipids. Over-expression of numerous adipogenic and lipogenic genes in the FL chickens suggests that in situ lipogenesis in chickens could make a more substantial contribution to expansion of visceral fat mass than previously recognized. Distinguishing features of the abdominal fat transcriptome in lean chickens were high abundance of multiple hemostatic and vasoactive factors, transporters, and ectopic expression of several hormones/receptors, which could control local vasomotor tone and proteolytic processing of adipokines, hemostatic factors and novel endocrine factors. Over-expression of several thrombogenic genes in abdominal fat of lean chickens is quite opposite to the pro-thrombotic state found in obese humans. Clearly, divergent genetic selection for an extreme (2.5–2.8-fold) difference in visceral fatness provokes a number of novel regulatory responses that govern growth and metabolism of visceral fat in this unique avian model of juvenile-onset obesity and glucose-insulin imbalance.


Introduction
For more than three decades, scores of papers have described various aspects of growth and nutrient metabolism in the divergently-selected FL and LL chickens originally developed by Leclerq et al [25]. In general, the FL and LL cockerels have similar growth rates with a 2.5-fold difference in abdominal fatness and higher breast muscle weights in the LL. The FL chickens always exhibit a lower plasma glucose level without overt hyperinsulinemia found in mammals, a peculiar condition which Simon et al. [30,31] described as a "glucose-insulin imbalance". Hypertriglyceridemia of the FL chickens indicates higher rate of hepatic lipogenesis from carbohydrate metabolism, most likely the consequence of a small increase in insulin-sensitivity in the FL chickens [26,31]. The FL and LL chickens are able to maintain their respective fat or lean phenotype independently of altered energy sources [32], emphasizing genetic regulation of phenotypic expression. These metabolic peculiarities in our polygenic model of juvenileonset obesity have been extensively examined by nutritional and metabolic perturbations [26,[33][34][35][36]. Furthermore, transcriptional profiling of multiple tissues [2,18,24,[37][38][39][40][41] and high-throughput surveys of variations in genome sequence and structure, including extensive quantitative trait loci (QTL) and expression (eQTL) analyses [17,19,20,[42][43][44] in the FL and LL chickens have begun to identify causal genes and to unravel the genetic and molecular basis for their divergence in lipid metabolism and visceral fatness.
The present descriptive transcriptomics study, using RNA sequencing (RNA-Seq) analysis, was designed to expand our catalog of expressed adipose genes at 7 wk with a dual goal (1) to determine the most transcriptionally-active biological processes in abdominal fat and (2) to establish major functional differences between the abdominal fat transcriptomes of FL and LL chickens. First, a functional characterization of the 900 highest expressed (HE) adipose genes independent of genotype was provided by a bioinformatics analysis. Second, we identified 1687 differentially-expressed (DE) genes from the comparison of transcripts in abdominal fat of FL and LL chickens. Ingenuity Pathway Analysis has revealed the over-expression of numerous genes involved in hemostasis, lipid catabolism, and endocrine signaling in the LL. In contrast, the up-regulation of several key adipogenic and lipogenic genes in abdominal fat of the FL chickens suggests that in situ lipogenesis could make a more substantial contribution to the expansion of adipose mass in the chicken than previously recognized.

Animals and tissue preparation
The FL and LL chickens were bred and raised at INRA UE1295 Pôle d'Expérimentation Avicole de Tours, F-37380 Nouzilly, France, as described previously [24]. Briefly, 8 birds from each genotype (FL and LL) were randomly selected for tissue sampling at six ages (1, 3, 5, 7, 9, and 11 wk), weighed, bled into heparinized syringes, and killed by cervical dislocation. Abdominal fat was quickly dissected, weighed, a sample was immediately snap frozen in liquid nitrogen, and stored at −75°C for further processing. All animal procedures were performed under the strict supervision of a French government veterinarian and in accordance with protocols approved by the French Agricultural Agency, the Scientific Research Agency, and the Institutional Animal Care and Use Committee at INRA, Nouzilly, France. These procedures were also in compliance with the United States Department of Agriculture guidelines on the use of agricultural animals in research and approved by the University of Delaware Agricultural Animal Care and Use Committee. [45] followed by DNase I treatment. Sample quality was analyzed with an RNA 6000 Nano Assay kit and the Model 2100 Bioanalyzer (Agilent Technologies; Palo Alto, CA). The rRNA ratio (28S/18S) was determined and all samples had an RNA integrity number (RIN) greater than 9.0. Sequencing libraries were made from 1 μg of total adipose RNA with the Illumina RNA Sample Prep Kit v2 following standard Illumina protocols. Individual RNA samples were indexed (bar-coded) to enable multiplexing of libraries within sequencing lanes. Libraries were pooled and sequenced using an Illumina HiSeq 2000 Sequencing System at the Delaware Biotechnology Institute, University of Delaware. Three separate schemes were used for paired-end (101 bp) sequencing of 8 libraries (4 FL and 4 LL) across two sequencing lanes per run. In Scheme A, two sequencing lanes were used for multiplexing of two FL and two LL samples per lane (n = 4/lane). Two libraries (1 FL and 1 LL) in sequencing lane 2 of Scheme A had low quality control (QC) scores and were eliminated from further analyses. Consequently, the two low QC libraries were re-sequenced in individual lanes in Scheme B (n = 1/lane). Finally, all eight (4 FL and 4 LL) libraries were multiplexed and sequenced in two replicate lanes in Scheme C (n = 8/lane). All samples in Schemes 1-3 were merged into one file of 12 samples from the FL and 12 samples from the LL cockerels for further analysis.

RNA sequence (RNA-seq) analysis
All reads generated from the three sequencing schemes (12 FL and 12 LL) described above were included in the RNA-Seq analysis using CLC Genomics Workbench 5.1 software (CLC bio, Cambridge, MA). The data analysis included sequence data filtering, read mapping, transcript and gene identification, analysis of differential gene expression, and functional annotation.
Sequence data filtering. Twenty-four short-read (101 base pairs) sequencing samples (12 FL and 12 LL) from the 3 sequencing schemes were de-multiplexed and imported into CLC Genomics Workbench, separately. Several QC trimming methods were used within the CLC Genomics Workbench software, including quality trimming, ambiguity trimming and adapter trimming with default settings applied before mapping to the reference chicken genome.
Read mapping and transcript/gene identification. The reference genome for the chicken (Gallus gallus, build 2.1) in FASTA format and the corresponding annotation file in GTF format were obtained from Ensembl (ftp.ensembl.org/pub/release-64), which represents 17,934 genes and 22,298 transcripts. Two hundred nucleotides of flanking region upstream and downstream of known genes were also included in the analysis. The short paired-end read sequences (101 bp x 2) were mapped to the reference chicken genome sequence, with mapping parameters that enforced: (1) a maximum of two mismatches and (2) reads must map with 90% of the bases aligned to the reference sequence with 80% similarity. Non-specific matches (reads mapped to multiple places in the reference genome) were excluded from the analysis.
Differential expression analysis. The unique exon reads count (including the exon-exon and exon-intron junctions) for the reads mapped to a gene and its flanking regions were used as the raw expression value for that gene. This raw expression value was normalized to the median of the total mapped reads across the 24 samples to account for variation in original library concentration and multiplexing number. The 24 sequencing samples were divided into two genotypes (FL and LL), resulting in 12 replicates for each genotype. Normalized expression values were analyzed as a beta-binomial model [46] to detect differential expression. The twosided P-value was corrected using the false discovery rate (FDR) adjustment to account for multiple hypothesis testing procedures [47]. Genes with FDR-adjusted P-value (0.05) were considered to be statistically significant. To ensure the biological relevance, a condition of fold change 1.2 (or -1.2) was added on top of FDR-adjusted P-value (0.05); Genes with FDR-adjusted P-value (0.05) and fold change 1.2 (or -1.2) were considered to be differentially expressed (DE) in this study. The fold-change threshold (±1.2-fold) for DE genes is based on our extensive experience in functional genomics and transcriptional profiling of multiple tissues from various chicken models. A recent RNA-Seq study of breast muscle in chickens afflicted with "Wooden Breast" disease [48] used a similar significance level [FDR-adjusted P-value (0.05) and ±1.3 fold-change] to identify DE genes using default settings for the Cuffdiff procedure in the open-source software, Cufflinks (http://cole-trapnell-lab.github.io/ cufflinks/).

Availability of supporting data
The RNA-Seq reads in Sequence Read Archive (SRA) format were deposited into the National Center for Biotechnology Information Gene Expression Omnibus (NCBI GEO) database under the accession # GSE42980. Further data sets supporting the present results are included within the article and in supporting information files.

Quantitative RT-PCR analysis
For verification of expression, quantitative real-time PCR (qRT-PCR) analysis was performed on a subset of 47 DE genes identified by RNA-Seq analysis. First-strand cDNA synthesis was performed by incubation of a 13 μl reaction (containing 1 μg of total DNase-treated RNA, 1 μl of 100 μM oligo dT 20 , 1 μl of 10 mM dNTP mix, and water to 13 μl total volume) for 5 min at 70°C and placed on ice for 2 min. A master mix containing 5 μl of 5x first-strand synthesis buffer, 1 μl of 0.1 M DTT, 1 μl of RNaseOUT, and 200 U of SuperScript III reverse transcriptase (Invitrogen, Carlsbad, CA) was added (final reaction volume of 20 μl). Primers were designed for qRT-PCR using Primer Express v2.0 software (Applied Biosystems, Foster City, CA). Detailed information for each primer pair, including gene name, gene symbol, forward and reverse primer sequences, GenBank accession number and amplicon size, are provided in S1 Table. An ABI Prism Sequence Detection System 7900HT was used to perform the qRT-PCR assays, using 10 ng of total RNA, Power SYBR green PCR master mix (Applied Biosystems, Foster City, CA), and 400 nM of each primer pair (Sigma-Aldrich, St. Louis, MO) in duplicate wells. A disassociation step was used to validate specific amplification and verify absence of primer dimers. PCR products were analyzed using agarose gel electrophoresis to compare product size to the expected amplicon size. The cycle time (Ct) for each sample was normalized to the corresponding sample geometric mean of housekeeping genes [49]. We selected two housekeeping genes [pantothenate kinase 1 (PANK1) and ribosomal protein L14 (RPL14)] based on their invariability in qRT-PCR analysis and identified as the most stably-expressed genes using RefFinder software (http://www.leonxie.com/referencegene.php). The 2 -(ΔΔCt) formula was used to calculate relative abundance of transcripts [50]. The statistical analysis of normalized qRT-PCR data across three ages (3, 7 and 11 wk) was performed using a general linear model (GLM) procedure in Statistical Analysis System (SAS v9.3; Cary, NC). These data were analyzed using a two-factor analysis of variance to determine main effects (P0.05) of genotype (G), age (A), and the interaction of age with genotype (A × G). Where genes were only used for qRT-PCR analysis at one age (7 wk), a student's T-test was used to identify significant (P0.05) differences between the FL and LL genotypes.

Mapped reads and detection of genes and transcripts
Sequence data from the 24 samples were mapped to the reference genome (Gallus gallus, build 2.1). Table 1 presents a summary of the RNA-Seq analysis including the number of mapped reads and detection of corresponding chicken genes and transcripts (see S2 Table for more  details). The original sequencing run (Scheme A) was completed by multiplexing 2 FL and 2 LL samples (N = 4) in two separate sequencing lanes. Two samples (1 FL and 1 LL) in Lane 2 had low quality scores; therefore, the low-quality data was eliminated from further RNA-Seq analysis. Nonetheless, Scheme A provided an average of 32.5 M mapped reads for the FL (N = 3) and 36 M mapped reads for the LL (N = 3), which allowed detection of 73% genes and 65% of transcripts across the FL and LL genotypes. Subsequently, the two low quality score libraries were re-sequenced in separate lanes in another sequencing run (Scheme B), which gave the highest average detection levels of 78% for genes and 70% for transcripts across both genotypes. In Scheme C, the 8 libraries (4 FL and 4 LL samples) were multiplexed and sequenced in duplicate lanes within the same sequencing run. Scheme C represents the balanced block design with two technical replicates as described by Auer and Doerge [51] for proper statistical analysis of RNA-Seq experiments. Comparing sequencing depths (averaged across the FL and LL chickens), Scheme C allowed detection of 71% of genes and 63% of transcripts by multiplexing 8 libraries per lane which were sequenced in duplicate lanes. The average number of reads mapped across Scheme A, B and C was greater for the LL (41.6 M) than the FL (35.4 M), which is reflected in the slightly higher number of expressed genes found in abdominal fat of the LL cockerels. The overall average across genotype (FL and LL) and sequencing schemes (A, B and C) shows that 45% (38.5 M) of the total reads were mapped, which equates to identification of 74% of genes (13,265/17,934) and 66% of transcripts (14,724/22,298) based on the reference chicken genome. Genes from the reference chicken genome were mapped to UniProtKB accession numbers by the Protein Information Resource (PIR) ID mapping service [52]. The assigned fold-change values (i.e., FL/LL expression ratios) were based on the number of normalized reads from the RNA-Seq analysis. A power analysis was performed on this RNA-Seq dataset using the web-based software program "Scotty" (http://euler.bc.edu/marthlab/scotty/scotty.php) [53]. This analysis demonstrates that our sample size of four birds/genotype (n = 4), sequenced across three depths [38.5 million mapped reads per sample averaged across depths (Table 1)], had power to detect 80% of genes expressed in each genotype (P0.01) at 1.5 fold difference and greater than 90% at a fold-change of 2 (S1 Fig). Further, the "Scotty" program performed hierarchical clustering using Spearman correlation as the distance metric. This correlation analysis grouped the two genotypes distinctly, where the individuals within each genotype (FL and LL) were closely linked.

Abdominal fat transcriptome of fed FL and LL chickens
First, the 900 highest-expressed (HE) genes, defined as genes with an average (across both genotypes) of >4289 reads/gene, were identified from RNA-Seq analysis of abdominal fat in FL and LL chickens at 7 wk (S3 Table). Of the HE genes, 164 were expressed higher in the FL (>1.2-fold difference), while 155 HE genes were up-regulated (< -1.2-fold difference) in visceral fat of the LL. Second, we identified 1687 DE genes with a FDR-corrected P-value (P0.05) and fold change 1.2 (or -1.2); and of these, 1182 DE genes were expressed higher in abdominal fat of LL chickens, whereas only 505 DE genes were expressed higher in FL chickens (S4 Table). A working list of 607 functional genes, associated with lipid metabolism (lipogenesis, lipolysis, lipid transport, etc.), was compiled from the RNA-Seq datasets using Ingenuity Pathway Analysis software (S5 Table). The Venn diagram (Fig 1) shows the intersection among HE genes, DE genes, and lipid metabolism genes, with 25 genes in common across all three gene sets. A total of 164 DE genes were shared between HE and DE gene lists. And, 108 DE genes are shared between the 1687 DE genes and the 607 lipid metabolism genes. Further, 87 genes are shared between the 900 HE genes and the 607 lipid metabolism genes.

Ingenuity Pathway Analysis (IPA) of gene expression
Analysis of highest expressed genes in abdominal fat of FL and LL chickens. The 900 highest-expressed (HE) genes in abdominal fat (S3 Table) were submitted to the Ingenuity Knowledge Base (http://www.ingenuity.com/) for functional annotation, mapping of genes to canonical pathways, and identifying gene interaction networks. There were 828 DE genes identified by IPA as "Analysis ready" (i.e., annotated in the Ingenuity Knowledge Base). A summary of the over-represented IPA functional categories found for 828 "Analysis ready" HE genes is presented in Table 2 Table 2) are also of particular interest. Of the 21 HE genes associated with "Hepatic Fibrosis", 17 HE genes are expressed higher in abdominal fat of the LL chickens, including 5 collagen genes (COL3A-COL6A), fibronectin 1 (FN1), spondin 1 and 2 (SPON1, SPON2), fibrillin 1 (FBN1), transforming growth factor, beta receptor II (TGFBR2) and thrombospondin 1 (THBS1). The "PPAR/RXR Activation" category includes 28 genes that are highly expressed in the FL and involved in lipid synthesis [thyroid hormone responsive spot 14 alpha (THRSPA), stearoyl-CoA desaturase (SCD1), fatty acid synthase (FASN), sterol regulatory element binding factor 2 (SREBP2), lipoprotein lipase (LPL) and glycerol-3-phosphate dehydrogenase 1 (GDP1)]. Likewise, the "Mechanism of Gene Regulation by PPARs" category includes 18 HE genes.
Many HE genes found in abdominal fat of the FL and LL chickens are key transcriptional regulators of lipogenesis and adipogenesis (Fig 2). For example, several transcription factors [SREBF2, THRSP, nuclear receptor subfamily 1, group H, member 3 (NR1H3) or liver-activated receptor alpha (LXRA), and peroxisome proliferator-activated receptor gamma (PPARG)] interact with each other and ultimately effect the transcription of several downstream target genes (Fig 2A). Some HE targets of SREBF2 include fatty acid desaturase 2  Table), differentially expressed (DE) genes (S4 Table), and Ingenuity annotated genes known to be involved in lipid metabolism (S5 Table). The numbers in overlapping arcs indicate the number of genes shared between and among these three categories.
Analysis of differential gene expression between FL and LL chickens. The dataset of 1687 DE genes (S4 Table) was submitted to Ingenuity Pathway Analysis (IPA) (http://www. ingenuity.com/) for functional annotation and mapping to canonical metabolic and regulatory pathways. From these, a total of 1322 DE genes were identified as "Analysis ready" by IPA. A summary of the Ingenuity Pathway Analysis of these DE genes is presented in Table 3, which includes the major functional categories: "Top Canonical Pathways, Diseases and Disorders, Molecular and Cellular Functions, Physiological System Development and Function, Top Gene Interaction Networks, Top Tox [Toxicology] Lists, and Top 10 Up-regulated and Down-regulated Genes". Annotated lists are provided for DE genes assigned by IPA to the "Top Canonical Pathways" (S6 Table) and "Top Tox Lists" (S7 Table) functional categories. A total of 1687 differentially-expressed (DE) genes from RNA-Seq analysis were submitted to IPA, which provided 1322 "Analysis Ready" DE genes for functional annotation and mapping to canonical pathways and gene interaction networks. P-values were determined by IPA software using Functional gene interaction networks were identified by Ingenuity Pathway Analysis (IPA). Genes are colored based on fold-change values determined by RNA-Seq analysis, where the red-color symbols signify higher expression in FL chickens and green-color gene symbols indicate higher expression in LL chickens. The false discovery rate (FDR) and fold-difference cutoff were not used in this functional analysis of the highest-expressed (HE) genes, which is simply based on a high number of reads mapped to known transcripts. Each gene was assigned a shape and function by IPA as shown in the "Network Shapes" legend inset. The direct gene interaction network in the top panel (A) was functionally annotated by IPA as "Lipid Metabolism, Molecular Transport, and Small Molecular Biochemistry", which emphasizes transcriptional regulation of adipogenesis and lipogenesis. The direct gene interaction network in panel (B) was functionally annotated by IPA as related to "Cellular Development, Cellular Growth and Proliferation, and Cellular Movement". Fisher's Exact Test as described by Ingenuity. The percent overlap and ratio were calculated for the number of observed genes compared to the number of known genes in the Ingenuity Knowledge Base for that category. Adipose genes with positive expression ratios are expressed at higher levels in FL cockerels, whereas genes with negative expression ratios are more abundant in abdominal fat of the LL. The top canonical pathway populated by 23 DE genes found in abdominal fat of the FL and LL chickens at 7 wk was the "Adipogenesis Pathway" (Table 3). Among the 8 DE genes up-regulated in the FL were frizzled class receptor 6 (FZD6), fibroblast growth factor receptor 3 (FGFR3), lipin 1 (LPIN1), CCAAT/enhancer binding protein (C/EBP) alpha (CEBPA), and perilipin 1 (PLIN1) (S7 Table), A total of 15 DE genes were expressed higher in the LL chickens; and among these, 12 genes encode transcription factors including Kruppel-like factor 5 (KLF5), nuclear receptor subfamily 2, group F, member 2 (NR2F2 or COUPTFII), early B-cell factor 1 (EBF1), SMAD family member 5 (SMAD5) and hypoxia inducible factor 1, alpha (HIF1A), The "Glucocorticoid Receptor (GR) Signaling" pathway contained 36 DE genes (10 genes up-regulated in the FL; 26 genes up-regulated in the LL). The highly-expressed DE genes in the FL included annexin A1 (ANXA1), plasminogen activator, urokinase (PLAU), CEBPA, Harvey rat sarcoma viral oncogene (HRAS), two RNA polymerases (TAF13 and POLR2H), and two heat-shock proteins (HSPA2 and HSPA8). A total of 14 transcription regulators related to GR signaling were more abundant in visceral fat of the LL, including progesterone receptor (PGR), androgen receptor (AR), nuclear receptor coactivator 2 and 3 (NCOA2; NCOA3), nuclear factor of activated T cells, cytoplasmic, calcineurin dependent 2 and 3 (NFATC2; NFATC3), and E1A-binding protein p300 (EP300). Of the 49 DE genes assigned to the "Axonal Guidance Signaling" pathway, 16 genes were up-regulated in the FL and 33 genes were more abundant in the LL, including 6 peptidases and 13 kinases. The "Hepatic Fibrosis" pathway also shows over-representation of 22 DE genes in visceral fat of the LL that encode 8 types of collagen, myosin, heavy chain (MYH10), 8 growth factor receptors [fibroblast growth factor receptor 2 (FGFR2), platelet-derived growth factor receptor, alpha and beta polypeptides (PDGFRA; PDGFRB), interleukin 1 receptor, type I (IL1R1), insulin-like growth factor 1 receptor (IGF1R), interferon gamma receptor 1 (IFNGR1), transforming growth factor, beta receptor II (TGFBR2), and angiotensin II receptor, type 1 (AGTR1)], and 4 growth factors (ligands) [(transforming growth factor, beta 2 (TGFB2), PDGFRB, tumor necrosis factor (ligand) superfamily, member 10 (TNFSF10) and c-fos induced growth factor (FIGF)]. In contrast, the FL chickens showed up-regulation of only five genes associated with fibrosis, namely three receptors (tumor necrosis factor receptor superfamily, member 1B (TNFRSF1B), endothelin receptor type B (EDNRB) and interleukin 1 receptor accessory protein-like 2 (IL1RAPL2), SMAD7 and platelet derived growth factor C (PDGFC)]. The "RXR Activation" pathway was similarly overrepresented with highly-expressed DE genes from the LL (23 genes, including retinoic acid receptor beta (RARB), EP300, and Janus kinase 2 (JAK2), whereas only retinoid X receptor gamma (RXRG), aldehyde dehydrogenase family 1, subfamily A2 (ALDH1A2) and SMAD7 were expressed higher in the FL birds.
RNA-Seq analysis revealed 29 GCPRs that were differentially expressed in abdominal fat of the FL and LL chickens at 7 wk of age (Fig 4B). Eight GCPR genes (EDNRB, FZD6, G proteincoupled receptor 132 (GPR132), latrophilin 3 (LPHN3), MC5R, neuropeptide Y receptor 2 (NPYR2), SSTR2 and UTS2R) were expressed at higher levels in the FL birds; whereas, 21 GCPRs were over-expressed in abdominal fat of the LL. Many of these differentially expressed  Fig 5B) and AR (Fig 5C) in the LL were identified by the Ingenuity Upstream Regulatory Analysis, which predicts that these two upstream regulators are inhibited (blue-colored symbols and edges), since their direct target genes are also down-regulated in the FL (i.e., expressed higher in the LL). This network of up-regulated genes, which interact with six transcription factors and numerous up-regulated direct targets of FOXA2 and AR in abdominal fat of the LL (Fig 5) clearly implicates intense transcriptional regulation of the lean phenotype. In particular, FOXA2 appears to play a critical role in restricting visceral adipose accretion in the LL, since it directly up-regulates expression of fibrinogen beta (FGB, a key blood-clotting protein), glucagon (GCG, a lipolytic pancreatic hormone), PCK1 (the key enzymatic regulator of gluconeogenesis), HMGCS2 (which catalyzes the first reaction in ketogenesis), two major transport proteins (ALB and TTR), and alpha-2-macroglobulin (A2M, a protease inhibitor and clinical biomarker of type 2 diabetes).
Another interesting gene interaction network involved in growth factor signaling was also highly represented by DE genes found in the LL adipose tissue (Fig 6A).  Table 3, under the IPA A comparison between the DE and HE gene sets has identified numerous candidate genes that control lipid metabolism and expansion of visceral fat (Table 4) Verification of RNA-seq analysis by quantitative RT-PCR Based on biological function, several candidate genes were selected from the RNA-Seq analysis for qRT-PCR verification (Table 5). Of the genes shown, 41 were differentially expressed in the FL/LL by RNA-Seq analysis (FDR-adjusted P-value 0.05). All 41 genes were significantly different (P 0.05) between FL and LL birds by qRT-PCR analysis. Neuropeptide Y (NPY) was a candidate gene not significantly different between the FL and LL chickens by either RNA-Seq or qRT-PCR analysis. The glucagon receptor (GCGR) was not differentially expressed by RNA-Seq; however, the abundance of GCGR in abdominal fat was 1.6 fold higher in the LL chickens by qRT-PCR analysis (P0.05). The short isoform of chicken growth hormone (scGH) was 1.5-fold higher (P0.05) in FL chickens by qRT-PCR analysis, while differential expression of scGH was not indicated by RNA-Seq analysis. Albumin (ALB; not shown) was highly up-regulated in abdominal fat of LL chickens by both RNA-Seq (38-fold) and qRT-PCR (73-fold) analyses. Furthermore, the fold change values for this gene set, across both RNA-Seq and qRT-PCR analyses, were highly correlated (r = 0.90) and significant (P0.001) as determined by Pearson correlation analysis. A complete list of the FL and LL expression means and their standard errors (SEM) for all genes in the qRT-PCR analysis is provided in S9 Table. The expression of six candidate genes [glucagon (GCG), glucagon receptor (GCGR), glucagon-like peptide 1 receptor (GLP1R), lysophosphatidic acid receptor 1 (LPAR1), attractin-like 1 (ATRNL1) and tubby (TUB)] was initially examined by qRT-PCR analysis at three different ages (3, 7 and 9 wk of age) to verify highest expression at 7 wk (Fig 7). A two-factor analysis of variance [genotype (FL and LL) x age (3, 7 and 11 wk)] shows significant (P0.05) main effects of genotype and age for all six genes. The abundance of all six genes reached a peak at 7 wk of age and was higher in abdominal fat of the LL than the FL chickens at 7 and 9 wk.

Discussion
The FL and LL chickens were originally developed by Leclercq et al [25] as experimental genetic models to identify the mechanisms controlling abdominal fatness, a complex polygenic trait likely governed by interactions among multiple genes controlling different endocrine and metabolic pathways. Between 7 and 9 wk of age, FL and LL chickens exhibit the greatest difference (2.5-to 2.8-fold) in abdominal fatness while maintaining similar body weights [24,41] and feed intakes [54]. This large divergence in abdominal fatness between FL and LL chickens at 7 wk also corresponds to the highest levels of gene expression in abdominal fat tissue [24]. To maintain a similar body weight despite large differences in abdominal fatness, FL chickens appear to favor partitioning of nutrients (particularly dietary amino acids) into abdominal fat; they also have a higher protein turnover rate [54]. Recent metabolomics studies have clearly demonstrated that the FL and LL chickens are able to maintain their fat and lean phenotypes independent of dietary energy sources, albeit the fatty acid composition of abdominal fat was affected by genotype [32,55]. On the other hand, LL chickens have a greater lean muscle mass and conversely a lower abdominal fat content than do the FL chickens. Several metabolic (liver, adipose tissue, and skeletal muscle) and regulatory (hypothalamus) tissues of the FL and LL chickens have been analyzed using numerous transcriptional methods (i.e., differential mRNA display, qRT-PCR, low-and high-density microarrays, etc.) to unravel the underlying mechanisms of excessive fatness [2,18,24,38,41,42,56,57]. To our knowledge, the present study represents the first deep RNA-Seq analysis of abdominal fat in these unique FL and LL chickens. The Ingenuity Upstream Regulator Analysis predicts inhibition of forkhead box A2 (FOXA2; blue-colored symbol with a z-score of -2.52), since 9 of its 17 direct targets have down-regulated expression ratios (FL/LL), which is consistent with inhibition of FOXA2 (i.e., downregulated in FL or up-regulated in LL). (C) Likewise, the predicted inhibition of the androgen receptor (AR; blue-colored symbol with a z-score of -3.04) would lead to inhibition (blue arrows) of its direct targets, since 15 of 27 genes are inhibited in the FL (i.e., a reduced FL/LL ratio or up-regulated in the LL). According to RNA-Seq analysis, the expression of both FOXA2 (fold-change of -10.7) and AR (fold-change of -1.4) was higher in the LL. Understanding the molecular mechanism of excessive abdominal fatness in the chicken requires examination of differential gene expression using different approaches. Only a few microarray studies have described the abdominal fat transcriptome of chickens [22,58,59]. The most relevant to the present study was the transcriptional analysis of abdominal fat in commercial broiler chickens after a short (5 h) period of acute fasting or immunoneutrilization of insulin [22]. Eighty DE genes identified by RNA-Seq analysis of abdominal fat in FL and LL chickens in the present study were also differentially expressed in visceral fat of broiler chickens subjected to 5 h of fasting [22]. Additionally, five DE genes (AGTR1, AIMP2, EEPD1, GCG and ICA1) that were identified in abdominal fat after acute insulin immunoneutralization of chickens [22] were also DE genes found presently in FL and LL chickens. Among these genes was GCG, the major contra-insulin regulator of glycemia in chickens [60][61][62], which was expressed higher in abdominal fat of the LL chickens and similar to that observed in chickens after acute insulin immunoneutralization [22]. The endocrine pancreas is the major site of GCG synthesis and secretion in birds [63], while adipose tissue is a major target of GCG, which stimulates lipolysis and release of non-esterified fatty acids into circulation. The GCG receptor, GCGR, showed a similar higher expression level in abdominal fat of LL chickens. In agreement with up-regulation of glucagon production and its lipolytic action in abdominal fat of LL chickens, the transcript for the somatostatin receptor (SSTR2), which when activated by somatostatin potently inhibits the production of GCG [64,65], was highly expressed and up-regulated in FL chickens. Furthermore, FOXA2, a transcriptional regulator of differentiation of glucagon producing cells [66], was also up-regulated in visceral fat of the LL chickens. We found that the glucagon-like peptide 1 receptor (GLP1R) was nearly 4-fold greater in LL chickens, whereas the glucagon-like peptide 2 receptor (GLP2R) was undetectable in abdominal fat of both LL and FL chickens. The absence of GLP2R transcripts in our RNA-Seq analysis could be of importance since its expression in adipose tissue of chickens was reported earlier [67]. The presence of GLP1R (and absence of GLP2R) in FL and LL chickens suggest that differential regulation occurs through the preproglucagon class A transcript (which contains GCG and GLP1) in abdominal fat of these birds rather than the preproglucagon class B transcript (which contains GCG, GLP1 and GLP2) [67]. Furthermore, only the preproglucagon class B transcript was upregulated in abdominal fat of chickens after insulin immunoneutralization [22]. Our observation of enhanced GCG and GLP1 signaling within abdominal fat of LL chickens supports the idea that visceral fat serves as an important endocrine organ, which displays"ectopic' expression of several endocrine factors and/or their respective receptors (i.e., GCG/GCGR, GLP1R, NPY/NPY2R, SSTR2, scGH, TUB, ATRNL1, CNR1 and MC5R).

Altered lipid metabolism in abdominal fat of FL and LL chickens
A recent review provides a comprehensive description of growth and energy metabolism in our divergently-selected FL and LL chickens, amassed from more than three decades of intensive investigation [28]. The "glucose-insulin imbalance" originally described by Simon et al [30,31] appears to be the major physiological difference noted between the FL and LL chickens, where the FL exhibit hypoglycemia and only slight hyperinsulinemia, which resembles the 'pre-obese' condition of juvenile mammals without true insulin resistance [26]. In fact, the FL chicken appears to be more sensitive to insulin [31] as indicated by their higher rate of glucose  utilization and enhanced hepatic lipogenesis [26,35], higher secretion of very-low-density lipoproteins (VLDL) [27], and their greater storage capacity of lipids in abdominal fat [68]. Furthermore, the divergence in abdominal fatness in the FL and LL chickens does not depend on the type of dietary energy, which indicates genetic control [32,55]. Since the liver is considered as the primary site of lipogenesis in birds [10][11][12], most previous studies of the FL and LL were directed at understanding the genetic difference in hepatic lipogenesis, rather than the lipogenic capacity of abdominal fat. A major finding of our time course microarray analysis of abdominal fat in FL and LL chickens [24] was the discovery of numerous DE genes involved in lipid synthesis. We have confirmed the up-regulation of many DE lipogenic genes in abdominal fat of FL chickens by both RNA-Seq and qRT-PCR analyses. Further, the abundant expression of lipogenic genes in abdominal fat in the FL was validated by the report that acute fasting [22] depressed expression of several lipogenic genes in the chicken. For example, the lipogenic enzyme ME1, which is critical for the synthesis of lipids via generation of NADPH required for the conversion of malonyl-CoA to palmitic acid (by FASN), was down-regulated in LL chickens and similar to its lower expression in adipose tissue of fasted chickens. Both the soluble (MDH1) and mitochondrial (MDH2) forms of malate dehydrogenase, which are expressed higher in adipose tissue and liver [57] of FL chickens, were among the top six percent of the highest expressed genes found by RNA-Seq analysis of abdominal fat in FL and LL chickens at 7 wk. Similar to the FL chickens [68], the expression of MDH1, MDH2 and ME1 show a strong positive correlation with increased adipocyte volume in genetically fat pigs [69]. Another gene differentially expressed in response to either genetic or nutritional perturbation is lipin 1 (LPIN1), which was identified as a reciprocal regulator of triglyceride synthesis and hydrolysis in adipocytes of LPIN1 knockout mice [70]. In chickens, adipose LPIN1 expression is elevated in response to short-term (5 hr) food deprivation [22], a 30% energy restriction [71] or in "leaner" and slowgrowing Leghorn and Fayoumi breeds when compared to "fatter" broiler chickens [23]. However, LPIN1 was expressed higher in abdominal fat of our genetically fat (FL) chickens as indicated by both RNA-Seq and qRT-PCR analyses (see Fig 6B). The present study clearly shows that divergent genetic selection of broiler chickens for a 2.5-2.8 fold difference in abdominal fatness results in higher expression of LPIN1 in visceral fat of the FL, but not the LL chickens. Of our cross-validated candidate genes, LPIN1 is of particular interest for its potential to control several important metabolic processes including adipogenesis, lipogenesis and eicosanoid signaling. Involvement of LPIN1 in the regulation of adipogenesis occurs via control of phosphatidic acids levels, which influence PPARG expression [72] and LPIN1 is required for PPARG-driven adipogenesis both in vitro and in vivo [73,74]. Further, LPIN1 is also a co-regulator of nuclear receptor signaling by interaction with the PPARs (PPARA, PPARG and PPARGC1A) and CEBPA [75,76]. We found that atrial natriuretic peptide receptor 3 (NPR3) and AMP-activated protein kinase (AMPK) were both up-regulated in abdominal fat of the FL chickens (see Table 4). Atrial natriuretic peptide (ANP) provokes a strong lipolytic action in human adipocytes by activating AMPK and increasing mitochondrial biogenesis [77].
In the present study, we found that the expression of KLF5 was 2.5-fold greater in abdominal fat of the LL than in the FL at 7 wk (Table 3; S8 Table). The transcription factor KLF5 has been implicated in adipogenesis, since knockout of KLF5 retards expansion of white (visceral) adipose tissue [78]. Furthermore, this study shows that over-expression of KLF5 stimulates adipocyte differentiation via interactions with CEBPB and PPARG, which are key regulators of adipocyte differentiation. A recent study of KLF2 expression in abdominal fat during juvenile development (1-12 wk) of a different population of divergently-selected FL and LL broiler chickens shows that the abundance of KLF2 was greater in the LL at 1 wk, but higher in fat of the FL between 3-8 wk [79]. Functional in vitro analyses indicate that over-expression of KLF2 inhibits preadipocyte differentiation and expression of CEBPA and PPARG. These observations support their conclusion that KLF2 plays a role in divergent genetic selection of this distinct population of FL and LL chickens. Clearly, the role of multiple transcription factors in the divergence of abdominal fat accretion in FL and LL chickens needs further investigation.
Hormone signaling appears to play a major role in controlling lipid metabolism in abdominal fat of FL and LL chickens. For example, growth hormone (GH) signaling is highly active in abdominal fat of the LL birds (see Fig 2B). We found that the peripheral receptor for melanocyte-stimulating hormone and adrenocorticotropic hormone, MC5R, was expressed higher in abdominal fat of FL chickens. MC5R regulates fatty acid oxidation in skeletal muscle of mice [80] and in hepatocytes of sea bass, where a MC5R agonist stimulates lipolysis and the release of free fatty acids from cultured sea bass adipocytes [81]. Furthermore, MC5R has been associated with obese phenotypes in humans [82]; and the MC5R gene locus is highly conserved in the genomes of human, mouse and chicken [83]. Additionally, the membrane-bound progesterone receptor (PGRMC1) was highly expressed in abdominal fat of FL chickens (see Table 4). The PGRMC1 directly regulates cholesterol synthesis and metabolism of steroid hormones in mammals [84]. Progesterone signaling seems to be enhanced in abdominal fat of LL chickens as indicated by the 2-fold up-regulation of the progesterone receptor (PGR), which regulates Comparisons are shown for 41 DE genes from the RNA-Seq analysis and three additional genes that were not differentially expressed by RNA-Seq analysis, but of interest [NPY, GCGR, and short chicken growth hormone (scGH)]. Fold changes were derived from FL/LL expression ratios of 4 birds/ genotype at 7 wk. Positive values correspond to up-regulation in FL chickens, while negative ratios indicate higher expression in abdominal fat of LL cockerels. Pearson correlation analysis showed a high correlation (r = 0.90) and significance (P0.001) between the relative expression ratio obtained by the two independent analytical methods (qRT-PCR and RNA-Seq). No significant difference (N.S.) was found for abundance of NPY transcripts between FL and LL chickens by either qRT-PCR or RNA-Seq analyses. The FL and LL expression means and their standard errors are provided in S9 Table for genes verified by qRT-PCR analysis. doi:10.1371/journal.pone.0139549.t005 the expression of 25 DE target genes in visceral fat. Progesterone administration in rats increases body and inguinal fat mass, a response only observed in females [85], whereas the significance of up-regulation of PGR in abdominal fat of LL cockerels in the present study remains unknown.
Androgen signaling also appears to be up-regulated in abdominal fat of LL chickens, where the AR was expressed 40 percent higher compared to their fatter counterparts. Androgen signaling in adipose tissue of murine models protects against obesity and regulates insulin action and glucose homeostasis, and AR knockout mice are more susceptible to high-fat diet-induced visceral obesity [86]. The up-regulation of AR and the 46 DE genes that are direct targets of AR found presently in abdominal fat of the LL suggests a similar mechanism for AR signaling in adipose tissue of chickens. Unfortunately, a recent study of gene expression in castrated (capon) chickens did not include a transcriptional analysis of abdominal fat in the capons, which become fatter, nor in the leaner intact or testosterone-implanted groups [87]. This report clearly shows the importance of androgen in limiting excessive fattening in the chicken, albeit more attention was focused on over-expression of PCK1 in the liver of capons rather than in abdominal fat as a direct target of androgen signaling.
Another intriguing finding in abdominal fat of FL and LL chickens was the abundance of genes involved in eicosanoid signaling. This includes production of eicosanoid precursors as well as downstream signaling events. In human breast cancer MCF7 cells, loss of FADS2 function blocks normal polyunsaturated fatty acid biosynthesis resulting in the FADS1 generation of polyunsaturated fatty acids which lack the 8-9 double bond of eicosanoid signaling precursors (i.e., arachidonic acid and eicosapentaenoic acid) [88]. Both FADS2 (higher in FL chickens) and FADS1 (higher in LL chickens) are differentially expressed in the present study. Another enzyme involved in eicosanoid signaling, DAGLA (higher in LL chickens), catalyzes the hydrolysis of diacylglycerol (DAG). The product of this reaction is arachidonic acid precursor 2-arachidonoyl-glycerol (2-AG), a major peripheral endocannabinoid signaling molecule [89]. In adipose tissue of rodents, 2-AG activates the CNR1 (a gene also up-regulated in LL chickens), an event that is up-regulated in fat treatment groups [89][90][91]. Furthermore, FAAH, the enzyme that degrades the active endocannabinoid 2-AG [92] was up-regulated in abdominal fat of FL chickens. Taken together, these findings provide evidence for enhanced endocannabinoid signaling in adipose tissue of the LL chickens.

Hemostatic gene expression in abdominal fat of chickens
The predominant expression of numerous hemostatic genes in abdominal fat of FL and LL chickens is clearly illustrated by an overlay on the coagulation cascade from the Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/kegg/) website (Fig 8). Seven hemostatic genes, up-regulated in the LL, were identified earlier in our time-course (1-11 wk) microarray analysis of abdominal fat in the FL and LL chickens [24]. The present RNA-Seq analysis of abdominal fat in the FL and LL chickens at 7 wk of age verifies and extends the remarkable involvement of the hemostatic system in divergent genetic selection of abdominal fatness in meat-type chickens. From our RNA-Seq analysis, an additional 11 genes [A2M, END1, END2, FGB, FGG, F2R, protein C receptor, endothelial (PROCR), SERPIND1, SER-PINB5, THSD1, and THSD4] were up-regulated in the LL chickens, whereas another 4 genes (PROS1, PLAU, SERPINF2 and EDNRB) were expressed at higher levels in visceral fat of the FL at 7 wk. The extraordinary involvement of the coagulation system in limiting expansion of fat mass in the LL chickens is further demonstrated by the high abundance of additional six hemostasis-related genes (VFW, TFPI2, F13A1, SERPINH1, THBS1 and THBSR) across both genotypes. Furthermore, our time-course microarray analysis of liver samples from the same FL and LL chickens used for transcriptional analyses of their abdominal fat failed to show any effects of the FL or LL genotype on hepatic expression of these coagulation factors [Unpublished work; see NCBI GEO accession GSE8812]. Interestingly, the coagulation factor XIII, A1 polypeptide (F13A1) was previously identified as a novel obesity candidate gene in humans [93]. Collectively, our previous time-course microarray analysis [24] and the present RNA-Seq analysis of abdominal fat in the FL and LL chickens at 7 wk support the idea that coagulation and angiogenic factors are critical components in divergent genetic selection of visceral fat mass, which seem to function locally (paracrine or autocrine) and independently of the hemostatic system found in systemic circulation, since the liver is the main source of coagulation factors found in the bloodstream. The recent RNA-Seq analysis of breast muscle from broilertype chickens afflicted with "Wooden Breast disease" has implicated the chicken's coagulation system in progression of this novel muscular disorder [48]. Birds suffering from this myopathy exhibit down-regulation of eight hemostatic genes (A2M, FGA, FGB, FGG, KNG1, SERPINC1 and VWF) and up-regulation of four other coagulation genes (F5, F10, PROS1 and F2R).
The serine protease thrombin (F2), which we identified as a DE gene in the time course microarray analysis of abdominal fat of the FL and LL chickens [24], is an extremely potent platelet agonist. In the present RNA-Seq study at 7 wk, F2 was not significantly different (1.3-fold higher in LL chickens; P 0.07; data not shown) between the genotypes. Interestingly, we found that the receptor for F2 [F2R or proteinase-activated receptor-1 (PAR1)] was expressed higher in adipose tissue of LL chickens. PAR1 and PAR4 are expressed in human adipose tissue, where treatment with F2 induces expression and secretion of several adipokines (IL-1B, IL-6, MCP-1, and VEGF) [94]; however, this is mediated through the PAR4 rather than PAR1. In another study, Kajimoto et al. [95] found that the interaction of F2 and PAR1 in 3T3-L1 adipocytes stimulates FABP4 which regulates the expression of interleukin 6 (IL-6) and vascular endothelial growth factor (VEGF). Although there have been no studies describing F2-PAR1 interactions in adipose tissue of chickens, an earlier report indicates that interaction between thrombin and PAR1 is required for "endothelial mesenchymal transdifferentation" in the aorta of chick embryos [96]. Our present study shows that genes involved in thrombin signaling are differentially expressed in visceral fat of FL and LL chickens. The function of these coagulation factors in regulating expansion of visceral fat mass in the chicken certainly needs further definitive investigation.
The up-regulation of numerous genes involved in hemostasis (pro-coagulation and anticoagulation) in visceral fat of LL chickens was originally reported in our time-course microarray analysis [24]. Although prothrombotic genes are associated with increased fatness in humans [97][98][99][100], our studies suggest a novel role of several hemostatic factors in limiting abdominal fat accretion in the LL chickens [24]. The present RNA-Seq analysis (at 7 wk) and our previous time-course (1-11 wk) microarray study [24] clearly show differential expression of numerous pro-and anti-coagulation genes in abdominal fat of FL and LL chickens (see Fig  8). Most remarkably, abdominal fat of LL chickens appears to be in a pro-coagulation state, since five genes driving thrombosis are up-regulated (F9, F2, FGA, FGB, and FGG). RNA-Seq analysis identified seven additional genes (F2R, FGB, FGG, PROCR, PROS1, PLAU, and SERPINF2) that were not found by our previous microarray analysis. Four related genes [serpin peptidase inhibitor, clade H, member 1, (SERPINH1), serpin peptidase inhibitor, clade B, member 5 (SERPINB5), thrombospondin, type I, domain containing 1 (THSD1) and thrombospondin, type I, domain containing 4 (THSD4)] were also identified by RNA-Seq analysis of abdominal fat. Likewise, five additional genes [endothelin 1 (EDN1), endothelin 2 (EDN2), endothelin receptor type B (EDNRB), thrombospondin 1 (THBS1) and thrombospondin receptor (CD36)] were found in the present study; these genes are usually involved in vasoregulation or activation of platelets and endothelial cells. We propose that these proteases and protease inhibitors found in visceral fat of the chicken could control the proteolytic activation or deactivation of adipokines and other endocrine factors. The down-regulation of PROS1, an inhibitor of blood coagulation, in the LL supports a pro-coagulation state in abdominal fat of these chickens. Further, the higher expression of PROC [24] and the PROC receptor (PROCR) in visceral fat of the LL chickens suggests increased activation of PROC and PAR1 (or F2R) by the multifunctional PROCR [101]. Fibrinogen (FGA, FGB and FGG) is mainly produced in the liver and secreted into systemic circulation. Interestingly, fibrinogen expression in FL chickens is nearly undetectable whereas all three fibrinogen subunits are expressed 10-fold higher in visceral fat of the LL. Similarly, the level of FGG in plasma was inversely related to adiposity in subcutaneous and visceral fat in humans [102] and in rats fed a high-fat diet for 8 wk [103]. Aside from its primary function in blood coagulation, fibrinogen has been identified as a binding surface for several proteins involved in vascular homeostasis [104]. The over-expression of fibrinogen and albumin (ALB) in abdominal fat of LL chicken could contribute to increased binding, uptake and/or transport of adipose-derived proteins. The differential expression of the large number of hemostatic factors observed in abdominal fat, but not in liver of our FL and LL chickens was not observed during nutritional or hormonal perturbation in the chicken [22], which suggests that these coagulation genes could have a novel role in divergent selection and the 2.5-fold difference in abdominal fatness observed in our unique experimental lines.
These hemostatic factors could contribute to local (paracrine/autocrine) and/or endocrine control of fat accretion and/or vasomotor tone, which would have important consequences on adipogenesis [105], angiogenesis and regulation of blood flow [106]. This apparent vascular regulation in visceral fat seems highly complex and under control of several mechanisms. One such mechanism is the renin-angiotensin system (RAS), which usually serves as a major systemic regulator of blood pressure; now, it is known that many components of RAS are also expressed in adipose tissue [107]. Support for this idea is provided by the present study and our observation that 11 of the 13 DE genes found in the renin-angiotensin signaling pathway, are expressed higher in abdominal fat of the LL chickens. As an active RAS ligand, angiotensin II is involved in lipid metabolism via regulation of insulin signaling in adipose tissue [105,108]. RAS plays an important role in the regulation of local and systemic blood pressure via vasoconstriction of peripheral arterioles. The G-coupled protein receptor for angiotensin II, AGTR1, was up-regulated in chickens after insulin immunoneutralization or acute fasting [22]. In the present study, we found that both AGTR1 and AGTR2 were down-regulated in adipose tissue of FL chickens (see Table 4). The down-regulation of angiotensin II receptors in adipose tissue of anabolic chicken models (FL, fed or normal insulin activity) suggests increased vasodilation in abdominal fat, a process that accompanies angiogenesis [109], allowing for increased blood flow and expansion of adipose mass. Furthermore, an overactive RAS could inhibit adipocyte differentiation [105] resulting in down regulation of adipogenesis in LL chickens.
Vascular tone is also mediated through the release of vasoconstrictors or vasodilators independent of the RAS, including catecholamines, endothelins and nitric oxide. Reduced expression of potent vasoconstrictors (EDN1 and EDN2) in adipose tissue of FL chickens strengthens the importance of increased vasodilation and therefore enhanced blood flow to support higher accretion of abdominal fat in FL chickens. While the endothelin receptor B (EDNRB) transcript was slightly higher in FL chickens (1.5-fold), the abundance of EDNRA was more than 10-fold higher, which suggests that EDNRA is the more active isoform in adipose tissue of FL chickens. In mammals, the major effects of vasodilation are mediated through either down regulation of vasoconstrictors or increased generation of nitric oxide (NO). The production of NO by nitric oxide synthase (NOS) is usually driven by hypoxia [110], where hypoxic conditions in adipose tissue usually lead to overproduction of and ultimately resistance to NO [111]. Since HIF1A was down-regulated in visceral fat of FL chickens and the abundance of NOS similar between genotypes, hypoxia and over-production of NO do not seem to be an issue. Rather, regulation of NO signaling in these chickens could be controlled by phosphodiesterases, which are expressed at much lower levels in FL chickens (PDE1C, PDE3A, PDE5A, PDE9A and PDE10A). Phosphodiesterases catalyze the degradation of cGMP [112], ultimately inhibiting the vasodilatory action of NO, and have been implicated in the treatment of cardiovascular diseases [113], respiratory diseases [114] and metabolic syndrome [115]. For example, PDE1C, a non-selective phosphodiesterase, is amongst the highest expressed phosphodiesterases in rat β-islet cells and knockdown of this gene significantly increases insulin secretion [116]. Higher expression of another nonselective phosphodiesterase, PDE3A, in adipose tissue in humans was correlated with weight loss after gastric bypass surgery [117]. Furthermore, inhibition of cGMP-specific phosphodiesterases (i.e., PDE5A and PDE9A) is effective in cardio-protection [113]. Taken together, altered RAS and cGMP signaling within adipose tissue of the FL chickens could ultimately result in vasodilation and increased blood flow to support their higher fat deposition.
The up-regulation of growth factors (PDGFC and FGFR3; see Fig 6A) and GREM1 in abdominal fat of FL chickens could support enhanced angiogenesis. The higher expression of FGFR3 in FL chickens is similar to that observed in adipose tissue of ApoE knockout mice, which exhibit increased adiposity when fed a high fat diet [118]. The expression of GREM1, a bone morphogenetic protein antagonist, correlates with increased angiogenesis in humans with pancreatic neuroendocrine tumors [119]; whereas, the knock-down of GREM1 in human HK-2 cells increases BMP7 signaling activity [120]. Correspondingly, BMP5 and BMP7 were down-regulated in abdominal fat of the FL chickens, whereas BMP15 expression was strongly up-regulated in these birds. The expression of BMP5 is down-regulated under highly vascularized conditions, like those simulated by tumor cell lines (i.e., adrenocortical carcinoma and adrenocortical tumor cell lines [121]). Another angiogenesis-related gene endoglin (ENG), a co-receptor for TGFβ, was expressed higher in abdominal fat of the LL chickens (see Fig 3B  and Table 4). Knockout of this gene (Eng -/-) is lethal to mid-gestation mouse embryos due to severely impaired angiogenesis and cardiac development [122,123]. Further, plasma insulin and hepatic triglyceride levels are reduced in heterozygous Eng +/mice fed a high-fat diet [124]. In addition, NPY-stimulated angiogenesis and adipogenesis in the mouse are mediated through the NPY2R receptor, a pathway augmented by glucocorticoids [125]. This study using mice provides strong evidence that up-regulation of NPY secretion and NPY2R expression within visceral fat contributes to obesity and metabolic syndrome in mammals. Our qRT-PCR analysis clearly shows overexpression of both NPY and NPY2R in visceral fat of FL chickens, which could contribute to their enhanced adipogenesis and angiogenesis. This idea is supported by a recent report of higher expression of NPY and its receptors (NPY1R and NPY5R) in abdominal fat of fatter high-growth chickens than that of the slow-growing leaner chickens [126]. Furthermore, NPY treatment of cultured chicken adipocytes boosts NPY2R-mediated proliferation, adipogenesis and lipid accumulation, which implicate a functional role for adipose-derived NPY and the NPY2R in expansion of visceral fat in the chicken [127]. The combined over-expression of angiogenic genes and under-expression of vasoactive factors identified in the present study could promote increased blood flow to support enhanced growth of visceral adipose tissue in FL chickens.
Collectively, our previous time-course microarray study [24] and the present RNA-Seq analysis of abdominal fat in FL and LL chickens have identified several novel features of the transcriptional response to divergent genetic selection for extremes in visceral fatness. The assumption that avian and mammalian species share similar mechanisms that regulate energy metabolism and adipogenesis must be tempered with a heightened awareness of novel avian characteristics. Namely, the loss of multiple adipokine genes (LEP [7], PAI-1, TNFA, resistin and omentin [8]) from the chicken genome has a major impact on our understanding of the avian regulatory systems that control energy intake, metabolism and adiposity. Most discussions of insulin resistance, type-2 diabetes, inflammation, and obesity in mammals are centered on involvement of these major adipokines. Many of the unique aspects of metabolism, including fasting hyperglycemia and insulin insensitivity, are likely reflected by absence of these and perhaps other "mammalian" regulators in the chicken-the premier avian model [1][2][3][4].
In summary, our RNA-Seq analysis of visceral adipose tissue in meat-type chickens divergently selected for a large difference in abdominal fatness reveals higher expression of numerous hemostatic and lipolytic genes in genetically lean chickens, whereas abundant expression of lipogenic, angiogenic and adipogenic genes were found in fat chickens. In both genotypes, abdominal fat shows "ectopic" expression of numerous endocrine factors and/or receptors that could contribute to their divergence in visceral fatness. We propose that over-expression of multiple hemostatic genes encoding serine proteases and serine protease inhibitors in visceral fat of lean chickens could be involved in processing of adipokines and other endocrine factors. The present RNA-Seq analysis of abdominal fat in genetically fat and lean chickens provides strong support for the idea that adipose tissue is an important endocrine organ in chickens. In addition, the abundance and differential expression of several key lipogenic and adipogenic genes in genetically fat chickens suggests that abdominal fat could make a more important contribution to lipid synthesis in birds than previously recognized. This idea is supported by an earlier paper [128], which indicated that fatty acids synthetized in liver of chickens are subject to oxidation by liver and muscle, whereas fatty acids synthetized by abdominal fat are retained there for storage. The importance of functional candidate genes identified in visceral fat from the present transcriptional profiling and bioinformatics analyses will require further definitive studies to verify their importance in lipogenesis and adipogenesis, and to gain mechanistic insight to genetic control of adiposity in the chicken.
Supporting Information S1 Fig. Power analysis of the FL and LL abdominal fat RNA-Seq dataset. (A) A power analysis was conducted to demonstrate adequate biological samples size using the web-based software program called "Scotty" (http://euler.bc.edu/marthlab/scotty/scotty.php). Using the average of 38.5 million reads per sample, the power was calculated at 1.5, 2, or 3-fold change (-X change) differences between FL and LL chickens at a significance of P0.01. We achieved the power to detect 80% genes with 1.5-fold differences as indicated by the red broken line. (B) The "Scotty" program also performed a hierarchical cluster analysis using the Spearman correlation as the distance metric to demonstrate relatedness among the eight individual (4 FL and 4 LL) birds used in the RNA-Seq analysis. (TIF) S1 Table. Information on design of primers used for qRT-PCR analysis. For each primer used for qRT-PCR analysis, gene symbol, gene name, forward and reverse primer sequences, GenBank accession, and amplicon size (bp) are provided. Primers were designed from indicated NCBI GenBank accessions using Primer Express 2.1 software (Applied Biosystems, Inc.). (XLSX) S2 Table. Detailed summary of RNA-Seq analysis. This table provides a summary for all samples used in the RNA-Seq analysis. Information for each biological sample includes: multiplexing number, genotype, NCBI GEO sample ID, total input reads, single reads after trimming, paired-end reads after trimming, total reads mapped, single reads mapped, paired-end reads mapped, total reads unmapped, single reads unmapped, paired-end reads unmapped, genes with reads per kilobase of transcript per million mapped reads (RPKM) > 0, and transcripts with RPKM > 0. (XLSX) S3 Table. Highest expressed (HE) genes in abdominal fat of FL and LL chickens at 7 wk. The table contains a list of the top 5% of highly-expressed genes found in adipose tissue. Information provided for each gene includes feature ID, gene symbol, Entrez Gene name, fold change and average reads (i.e., number reads averaged across FL and LL chickens for each gene). No fold-difference or statistical cutoff was applied to the Highest Expressed (HE) Gene list (Worksheet 1). Two additional worksheets provide a list of HE genes with a < ±1. 2 This table provides the feature ID, gene symbol, gene description, FL Table. DE genes in abdominal fat of FL and LL Cockerels associated with lipid metabolism. This table contains 607 genes that are known to be involved in lipogenesis or lipolysis according to their annotation in the Ingenuity Knowledge Base. The "Average reads" column contains the number of reads from the RNA-Seq analysis of abdominal fat for that gene averaged across FL and LL genotypes. (XLSX) S6 Table. Top canonical pathways highly populated by DE genes from RNA-Seq analysis of abdominal fat in FL and LL cockerels at 7 wk. This Excel file contains five worksheets that list DE genes that populate the five top canonical pathways identified by IPA analysis (see Table 3). The top five canonical pathways from IPA analysis are: "Adipogenesis Pathway, Glucorticoid Receptor Signaling, Axonal Guidance Signaling, Hepatic Fibrosis, and RAR Activation". Each worksheet lists the following information: gene symbol, Entrez gene name, foldchange (FL/LL), cellular location, and type of molecule. (XLSX) S7 Table. Top toxicology (tox) functions identified by IPA analysis of DE genes in abdominal fat of FL and LL cockerels at 7 wk. This Excel file contains five worksheets of DE genes that populate the five toxicology functions identified by IPA analysis (see Table 2). The top five toxicology functions from IPA analysis are: "Cardiac Hypertrophy, RAR Activation, Liver Proliferation, Oxidative Stress, and TGF-β Signaling". Each worksheet lists the following information: gene symbol, Entrez gene name, fold-change (FL/LL), cellular location and type of molecule. (XLSX) S8 Table. Sixteen DE genes identified by Ingenuity Upstream Regulator Analysis are known transcription regulators. This Excel file contains a worksheet that list 16 DE genes identified by IPA as "Upstream regulators". The worksheet list the following information: gene symbol, fold-change (FL/LL), function of regulator, IPA activation z-score, P-value of overlap, the number of target genes, and list of all direct target genes of that transcription regulator in the DE dataset from RNA-Seq analysis. (XLSX) S9 Table. Verification of gene expression from RNA-Seq of abdominal fat in FL and LL cockerels by qRT-PCR analysis. This Excel file contains a worksheet that lists the genes used for qRT-PCR verification of differential gene expression, the average expression level ±SEM, and the significance level between four FL and 4 LL cockerels (see Table 5 for FL and LL means ±SEM). (XLSX)