Comparative Response of the Hepatic Transcriptomes of Domesticated and Wild Turkey to Aflatoxin B1

The food-borne mycotoxin aflatoxin B1 (AFB1) poses a significant risk to poultry, which are highly susceptible to its hepatotoxic effects. Domesticated turkeys (Meleagris gallopavo) are especially sensitive, whereas wild turkeys (M. g. silvestris) are more resistant. AFB1 toxicity entails bioactivation by hepatic cytochrome P450s to the electrophilic exo-AFB1-8,9-epoxide (AFBO). Domesticated turkeys lack functional hepatic GST-mediated detoxification of AFBO, and this is largely responsible for the differences in resistance between turkey types. This study was designed to characterize transcriptional changes induced in turkey livers by AFB1, and to contrast the response of domesticated (susceptible) and wild (more resistant) birds. Gene expression responses to AFB1 were examined using RNA-sequencing. Statistically significant differences in gene expression were observed among treatment groups and between turkey types. Expression analysis identified 4621 genes with significant differential expression (DE) in AFB1-treated birds compared to controls. Characterization of DE transcripts revealed genes dis-regulated in response to toxic insult with significant association of Phase I and Phase II genes and others important in cellular regulation, modulation of apoptosis, and inflammatory responses. Constitutive expression of GSTA3 was significantly higher in wild birds and was significantly higher in AFB1-treated birds when compared to controls for both genetic groups. This pattern was also observed by qRT-PCR in other wild and domesticated turkey strains. Results of this study emphasize the differential response of these genetically distinct birds, and identify genes and pathways that are differentially altered in aflatoxicosis.


Introduction
Aflatoxin B 1 (AFB 1 ) is a ubiquitous hepatotoxic, hepatocarcinogenic, and immunosuppressive mycotoxin. Poultry and other livestock are exposed to AFB 1 by consuming contaminated feed. Many agricultural feed commodities (corn, cottonseed, peanuts, and sorghum) and other foods (figs, tree nuts, and spices) are at especially high risk of being contaminated [1]. AFB 1 is practically unavoidable types, but other processes and pathways (i.e., apoptosis, cellular regulation, immune responses) are also important. Whereas, understanding the effects of AFB 1 on developing embryos is important in poultry production, the manifestation of AFB 1 toxicity is likely to be different in more mature birds with fully developed gastrointestinal systems. The purpose of this study was to compare the hepatic transcriptome response to dietary AFB 1 in juvenile (three weeks of age), susceptible (domesticated), and more resistant (wild) turkeys. We hypothesized that transcriptome responses in juvenile birds would reflect the more mature status of the gastrointestinal and antioxidant systems than those of embryos.

Results
Liver measurements were collected at the end of the exposure trial to characterize phenotypic effects of AFB 1 toxicity. Livers of domesticated (DT) birds (average = 20.54 g) were nearly three times the mass of those from Eastern wild (EW) birds (8.3 g) primarily due to differences in body size (average = 1147.5 g and 396.1 g, respectively). Liver weights of AFB 1 turkeys were smaller than those of the control (CNTL) groups. In DT, average liver mass at the conclusion of the trial ranged from 14.76 g to 23.98 g in the AFB 1 group (mean = 20.02 ± 2.44 g) and from 17.39 g to 25.11 g in controls (21.00 ± 2.12 g). Although this difference in liver mass was not significant (t-test p = 0.1962), when corrected for body weight (% BW) the livers of birds from the AFB 1 -treated group were significantly smaller (p = 0.0098). Livers of AFB 1 -treated EW birds (7.19 ± 1.06 g) were similarly smaller than those of control birds (9.43 ± 1.11 g). This difference was significant for absolute mass (p = 0.0005) and nearly so for % BW (p = 0.0531).
Sequencing of RNA libraries produced over 195 M reads. The number of reads per library ranged from 10.8 M to 14.6 M (average 12.2 M, Table 1). After trimming and filtering, median Q scores were consistently high and ranged from 36.4 to 37.4 among the forward and reverse reads. The number of reads per treatment group was balanced and ranged from 11.5 to 12.7 M, with an average of 12.21 ± 0.5 M. Approximately 91% of the quality-trimmed reads mapped to the annotated turkey gene set (NCBI Annotation 101, Table 1). This percentage was consistent across treatment groups and the percentage of aligned read pairs exceeded 89.7% (average = 90.8%); the majority of reads (average = 85.3%) mapped concordantly (Table 1). Based on mapping, the estimated mean library insert was 191.7 bp.
Principal component analysis (PCA) of normalized read counts visualized variation among the treatment groups ( Figure 1). Groups clustered distinctly according to treatment (AFB 1 versus control) within the first two principal components, accounting for approximately 95% of the observed variation. Hierarchical clustering of groups by Euclidean distance reiterated the relationships shown by PCA ( Figure S1). After segregating by AFB 1 treatment, groups secondarily cluster by type (domesticated versus wild), with the exception of samples N3L (domesticated) and EW1L (wild) that clustered with the opposite bird type. Significant differences in overall gene counts among groups are shown in the heat map of co-expressed genes. For each library the total number of concatenated reads, median read qualities (R1 and R2), estimated mean insert length (bp), number of and percentage of aligned reads, percentage of concordant reads, and the number and percentage of observed genes (mapped reads > 1) and expressed genes (mean group normalized read count > 3.0) are given.

Gene Expression
Evidence of expression (mean mapped reads ≥ 1.0 in at least one treatment group) was detected for 19,764 genes (tRNAs excluded), with an average of 17,137.5 genes being detected per group (81.56% of the turkey gene set) (Table S1). Mean read depth was 394.8 reads per gene. When limited to an average number of mapped reads ≥ 3.0, the number of expressed genes ranged from 14,399 to 17,230 among treatment groups (average 15,979.25, Table 1). Distribution of unique and shared expressed genes is illustrated in Figure S2. A total of 14,373 genes (81.2%) was co-expressed among all groups, and the number of co-expressed genes within the EW and DT lines was 14,908 and 14,669, respectively. Each treatment group had distinct sets of uniquely expressed genes, but that number was considerably greater in the AFB 1 -treatment groups (2169, 12.3%) when compared to controls (300, 1.6%, Table 1 and Table S1). had distinct sets of uniquely expressed genes, but that number was considerably greater in the AFB1treatment groups (2169, 12.3%) when compared to controls (300, 1.6%, Tables 1 and S1).

Differential Transcriptomic Expression: AFB 1 Effects
Table S2 provides the full list of genes showing significant differential expression (DE). DE was observed for 9620 genes (FDR p-value < 0.05, log 2 FC = −9.670 to 9.358) in wild turkeys exposed to AFB 1 when compared to control birds with 4176 genes having |log 2 FC| > 2.0 (Table 2). Similarly, 11,325 DE genes were observed for the AFB 1 -treated DT turkeys (FDR p-value < 0.05, log 2 FC = −14.133 to 12.676) with 4621 genes having |log 2 FC| > 2.0. The majority of DE genes (3380) were shared in both bird types, with 796 being unique to wild and 1241 unique to the domesticated birds ( Figure 2). The majority of DE genes was up regulated by AFB 1 (Table S3). In mammals, ANGPTL3 is in part involved in regulation of lipid and glucose metabolism by inhibiting the lipolysis of triglyceride-rich lipoproteins [29,30]. This transcript was highly down regulated in the AFB 1 -treated turkeys (log 2 FC = −7.94 and −14.13 in EW and DT, respectively). Similarly down regulated was GC an important protein in vitamin D transport and storage, actin-scavenging, and enhancement of complement component 5a activity for neutrophils in inflammation and during macrophage activation [31]. Comparison analysis in IPA found the most significant canonical pathways to include "Axonal Guidance Signaling" and "Hepatic Fibrosis/Stellate Cell activation". Hepatic stellate cells are closely linked to the progression of hepatic fibrosis [32]. For each comparison, the treatment groups, number of genes with significant FDR p-value, and the numbers of significant genes that also had |log 2 fold change| > 1.0 and > 2.0 are given.  Table S2 provides the full list of genes showing significant differential expression (DE). DE was observed for 9620 genes (FDR p-value < 0.05, log2FC = −9.670 to 9.358) in wild turkeys exposed to AFB1 when compared to control birds with 4176 genes having |log2FC| > 2.0 (Table 2). Similarly, 11,325 DE genes were observed for the AFB1-treated DT turkeys (FDR p-value < 0.05, log2FC = −14.133 to 12.676) with 4621 genes having |log2FC| > 2.0. The majority of DE genes (3380) were shared in both bird types, with 796 being unique to wild and 1241 unique to the domesticated birds ( Figure 2). The majority of DE genes was up regulated by AFB1 treatment in both the wild (2717, 65%) and domesticated birds (2914, 70%  (Table S3). In mammals, ANGPTL3 is in part involved in regulation of lipid and glucose metabolism by inhibiting the lipolysis of triglyceride-rich lipoproteins [29,30]. This transcript was highly down regulated in the AFB1-treated turkeys (log2FC = −7.94 and −14.13 in EW and DT, respectively). Similarly down regulated was GC an important protein in vitamin D transport and storage, actinscavenging, and enhancement of complement component 5a activity for neutrophils in inflammation and during macrophage activation [31]. Comparison analysis in IPA found the most significant canonical pathways to include "Axonal Guidance Signaling" and "Hepatic Fibrosis/Stellate Cell activation". Hepatic stellate cells are closely linked to the progression of hepatic fibrosis [32]. For each comparison, the treatment groups, number of genes with significant FDR p-value, and the numbers of significant genes that also had |log2 fold change| > 1.0 and > 2.0 are given.

Shared DE Genes
Changes in expression of the 3380 shared significant DE genes were highly correlated (r 2 = 0.909, F = 0.010, Figure S1) and essentially linear, except for genes with the greatest down regulation where log 2 FC in the domestic birds tended to be of greater magnitude than observed for the wild birds. This is consistent with a common physiological response to AFB 1 exposure. Of the 3380 shared DE genes, 1609 IDs mapped to the G. gallus gene REFLIST and statistical overrepresentation tests (PANTHER) of these shared DE genes found the greatest enrichment in the Biological Process category for "amino acid processes" and "negative regulators of hemostasis and wound healing" (Table S4). Comparison analysis in IPA found the most significant toxicology functions consistent with cellular damage. Categories with the greatest number of included genes included "liver hyperplasia/hyperproliferation" (p = 5.10 × 10 −41 ), "cardiac hypertrophy" (p = 2.91 × 10 −12 ), "renal necrosis/cell death" (p = 7.78 × 10 −12 ), and "liver steatosis" (p = 1.75 × 10 −11 ). Highest activation (Z) scores were obtained for "Integrin Signaling", "Rho Family GTPase signaling", and "NFAT regulation of immune response" pathways.
Only two loci among the 3380 shared significant DE genes showed opposite directional expression changes between wild and domesticated birds. CD96, a T cell-specific receptor, was significantly up regulated in response to AFB 1 in the EW birds (log 2 FC = 3.83) but down regulated in DT (log 2 FC = −2.15). CD96 may play a role in the adhesive interactions of activated T and NK cells when actively engaging diseased cells within areas of inflammation [33]. A second locus (LOC104911020, serum amyloid A protein-like) was significantly down regulated in response to AFB 1 in EW birds (log 2 FC = −3.99), but up regulated in DT (log 2 FC = 2.22). Serum amyloid A (SAA) proteins are a family of apolipoproteins produced primarily by the liver and are associated with high-density lipoprotein in plasma [34].

Unique Responses
Although the greatest number of DEGs was shared between the domesticated and wild turkey comparisons, 796 DEGs were uniquely affected in the wild birds exposed to AFB 1 in comparison to their controls ( Figure 2). Up-regulated genes with the greatest fold change in the EW birds (Table S3) (Table 3). Effected Cellular Component groups included "elements of the sarcolemma" (GO:0042383) and "proteinaceous extracellular matrix" (GO:0005578) enriched by 7.65-and 3.92-fold, respectively. A greater number of DEGs (1241) were uniquely affected in the domesticated birds exposed to AFB 1 , in comparison to their controls. Of the 50 DEGs with the greatest fold change in the domesticated birds, only two were up regulated; SMIM24 (small integral membrane protein 24) and LOC100546964 (cis-aconitate decarboxylase-like) ( Table S3). The function of SMIM24 is currently unknown. In humans, cis-aconitate decarboxylase (ACOD1 = IRG1) is highly expressed in mammalian macrophages during inflammation where it catalyzes itaconic acid production [35]. Among the most down-regulated loci in DT were ANGPTL3 (angiopoietin-like 3) which in humans, is expressed predominantly in the liver and functions in angiogenesis, and Fibrinogen (FGA, FGB, and FGG) and other coagulation components like coagulation factor IX (F9). This set of down-regulated DEGs also included LOC100547030 (a cytochrome P450 2W1-like gene). In humans, CYP2W1 is able to metabolically activate several pro-carcinogens, including AFB 1 , into cytotoxic products [36]. Due to its selective expression, CYP2W1 is suggested as a potential prognostic biomarker in hepatocellular and other carcinomas [37].
GO analysis of the DEGs unique to the domestic turkey liver indicate a number of distinctive responses to AFB 1 . Of the 1241 DEGs, 636 mapped to IDs in G. gallus REFLIST and overrepresentation tests in PANTHER found the greatest enrichment for biological process categories "organelle fission" (GO:0048285), "oxidation-reduction process" (GO:0055114), and "regulation of immune system process" (GO:0002682) ( Table 4). Cellular component categories were enriched for mitochondrial and membrane components.

Control
Comparison of the transcriptomes of control EW and DT birds found 774 DEGs (FDR p-value < 0.05, log 2 FC = −8.826 to 8.213), with 184 having log 2 FC > 2.0 ( Figure 3, Table S5). Of the 184 genes, seven were shared in common in the EW vs. DT AFB 1 comparisons (Figure 3). The shared loci included 5 genes up regulated in EW birds (ANGPTL3, CAMK4, LOC100538933 [probable ATP-dependent RNA helicase DDX60], LOC100545362 [ncRNA], LOC104912934 [ncRNA]), and two that were down regulated (LOC104910139 [ncRNA], LOC104915640 [KIAA1755 homolog]) when compared to DT. These shared genes are primarily metabolic and transcriptional regulators. In mammals, ANGPTL3 [angiopoietin-like 3] is a hepatokine involved in regulation of lipid and glucose metabolism and in the regulation of angiogenesis [37], CAMK4 [calcium/calmodulin-dependent protein kinase IV] is implicated in transcriptional regulation in immune and inflammatory responses [39], and DDX60 positively regulates DDX58/RIG-I-and IFIH1/MDA5-dependent type I interferon and interferon inducible gene expression [40].  [37], CAMK4 [calcium/calmodulin-dependent protein kinase IV] is implicated in transcriptional regulation in immune and inflammatory responses [39], and DDX60 positively regulates DDX58/RIG-I-and IFIH1/MDA5-dependent type I interferon and interferon inducible gene expression [40]. Of the 177 DEGs unique to the control group birds the majority (69%) were up regulated in the wild birds as compared to the domesticated birds. GO analysis found only a single biological process category ("single-organism process" GO:0044699) enriched among these genes (fold enrichment = 1.66, p-value = 4.15 × 10 −2 ). Interestingly, the up-regulated DEGs observed in the control-diet comparison also included glutathione S-transferase A3 (GSTA3). This gene was expressed at a 5.3-fold higher level (log2FC = 2.313) in the EW birds when compared to the DT birds ( Figure 4A), suggesting a higher constitutive expression in the former. Of the 177 DEGs unique to the control group birds the majority (69%) were up regulated in the wild birds as compared to the domesticated birds. GO analysis found only a single biological process category ("single-organism process" GO:0044699) enriched among these genes (fold enrichment = 1.66, p-value = 4.15 × 10 −2 ). Interestingly, the up-regulated DEGs observed in the control-diet comparison also included glutathione S-transferase A3 (GSTA3). This gene was expressed at a 5.3-fold higher level (log 2 FC = 2.313) in the EW birds when compared to the DT birds ( Figure 4A), suggesting a higher constitutive expression in the former. The genes showing the greatest positive fold change was the Interferon-inducible iron-sulfur clusterbinding antiviral protein RSAD2, radical S-adenosyl methionine domain containing 2, with log2FC = 8.2 and the interferon alpha-inducible protein 27-like 2B (IFI27L2B, LOC100550948). In mammals, RSAD2 can inhibit a wide range of DNA and RNA viruses, and has also been shown to play a role in CD4+ T-cells activation and differentiation [41], whereas IFI27L2B mediates virus-induced apoptosis [42]. These suggest heightened immune system activity in the wild birds.
Two genes involved in cell cycle control and chromosomal replication were also expressed at a significantly higher level in the wild bird controls. CD1, chromatin licensing and DNA replication factor 1 and MCM3, minichromosome maintenance complex component 3 had log2FC = 2.74 and 2.4, respectively (Table S5). Interestingly several other genes involved "cell cycle control" and "chromosomal replication" The genes showing the greatest positive fold change was the Interferon-inducible iron-sulfur cluster-binding antiviral protein RSAD2, radical S-adenosyl methionine domain containing 2, with log 2 FC = 8.2 and the interferon alpha-inducible protein 27-like 2B (IFI27L2B, LOC100550948). In mammals, RSAD2 can inhibit a wide range of DNA and RNA viruses, and has also been shown to play a role in CD4+ T-cells activation and differentiation [41], whereas IFI27L2B mediates virus-induced apoptosis [42]. These suggest heightened immune system activity in the wild birds.
Two genes involved in cell cycle control and chromosomal replication were also expressed at a significantly higher level in the wild bird controls. CD1, chromatin licensing and DNA replication factor 1 and MCM3, minichromosome maintenance complex component 3 had log 2 FC = 2.74 and 2.4, respectively (Table S5). Interestingly several other genes involved "cell cycle control" and "chromosomal replication" were also significantly up regulated in the wild birds, although with log 2 FC < 2.0 (Table S2). These included CDC45 (cell division cycle 45), DNA2 (DNA replication helicase/nuclease 2), the minichromosome maintenance complex components MCM2, MCM4, MCM5, MCM6, and MCM8, and POLE (DNA polymerase epsilon, catalytic subunit). These together with CD1 and MCM2 were components of the most significant Canonical Pathway identified by IPA (p = 1.1 × 10 −8 , ratio 0.25). Genes showing the greatest negative fold change were GYG2 (glycogenin 2), IQCD (IQ Motif Containing Protein D), and LOC100540418 (BPI fold-containing family C protein-like) (Table S5). GYG2 is involved in the initiation reactions of glycogen biosynthesis [43], and IQCD has been shown to interact with the RXR nuclear hormone receptor, and is thought to function as a transcriptional coactivator [44].

AFB 1 Treatment
Comparison of the transcriptomes of EW and DT birds fed the AFB 1 diet revealed 903 DEGs (FDR p-value < 0.05, log 2 FC = −7.987 to 9.262), of which, 143 had log 2 FC > 2.0 ( Figure 3, Table S5). Of these 143 DEGs, 136 were unique to the AFB 1 -treated birds. The majority (76%) of the DEGs were up regulated in EW relative to DT in a similar fashion to that seen for the control-group comparison. Unlike the control-fed groups, GO analysis found the DE genes from the AFB 1 between-line comparison highly enriched in several biological process categories (Table 5) including processes related to coagulation, inflammatory response and apoptosis. For example, up regulated in EW relative to DT were three components of the blood clotting factor fibrinogen (FGA, FGB, and FGG) and several other coagulation-related genes (vitamin K-dependent coagulation factor IX (F9), serpin peptidase inhibitor, member 10 (SERPINA10), histidine-rich glycoprotein (HRG), serpin peptidase inhibitor, clade C (antithrombin), member 1 (SERPINC1)). Although down-regulated genes did not show significant enrichment for a particular bioprocess, EW birds may produce less of the pro-inflammatory cytokine IL-17A (LOC100546746).
Expression of many enzymes responsible for xenobiotic metabolism, with an emphasis on those with specificity towards AFB 1 , was significantly altered by the AFB 1 treatment (Table 6). With a few exceptions, AFB 1 largely caused down regulation of these genes in both types of turkeys when compared to controls; within the control groups, the number of DEGs was considerably smaller and the magnitude of expression changes was also smaller, but directionally similar. Hepatic expression of CYP and GST genes was typically greater in EW vs. DT. Five cytochrome P450 loci; CYP1A5 (cytochrome P450, family 1, subfamily A, polypeptide 5, log 2 FC = 2.404), LOC100547030 (cytochrome P450 2W1-like, log 2 FC = 5.87), LOC100542486 (cytochrome P450 1A4, log 2 FC = 3.150), LOC100548433 (cytochrome P450 2K1-like, log 2 FC = 3.025), and LOC104915479 (cytochrome P450 2H1-like, log 2 FC = 2.135) were among the unique DEGs in the EW vs. DT comparison.
The turkey possesses six α-GST cluster genes, all of which possess detectable enzymatic activities toward prototype substrates in a recombinant expression system [19], unlike hepatic forms. Expression of the GSTAs was significantly altered by AFB 1 treatment. With the exception of GSTA3, GSTA1.1, 1.3, 2 and 4 were down regulated, while GSTA3 increased with dietary AFB 1 in DT but not EW (log 2 FC = 1.4667, Table 6). It is noteworthy that GSTA3 expression was significantly higher in EW birds when compared to DT birds for both control (log 2 FC = 2.3130) and AFB 1 (log 2 FC = 0.9211) group comparisons (Table 6, Figure 4A).
Expression differences in GSTA3 observed in RNA-seq read counts were confirmed by qRT-PCR. GSTA3 expression varied widely among treatment groups with experiment-wise threshold values (∆Ct) ranging from 17.27 to 26.06. Expression of GSTA3 transcripts was significantly higher in AFB 1 -treated birds than controls for both genetic groups (EW, p = 0.0061 and DT, p = 0.0036). Relative GSTA3 expression was also similarly variable in the other commercial (BB) and wild-type birds (RGW) ( Figure 4B) where GSTA3 expression was higher in AFB 1 -treated birds. This difference, however, was only significant in the BB comparison (p = 0.0015). Relative expression in BB birds was slightly higher than observed in the DT (Nicholas strain) birds and also higher in RGW birds (Rio Grande wild) when compared to EW. This result demonstrates that the differences observed in GSTA3 expression between the EW and DT birds is not unique to these genetic lines but is a broader, wild versus domesticated-bird phenomenon.   Genes included were significant (FDR p-value < 0.05) in at least one comparison. Comparisons highlighted in green are down regulated and those in red up regulated. Cytochrome P450 (CYP) and glutathione S-transferase (GST) family members shown to have in vitro activity towards AFB 1 and its metabolites in turkey are indicated in bold. * Due to similarity, these likely include transcripts assignable to GSTA1.1, A1.2, and A1.3.

Discussion
When compared to their domestic relatives, wild turkeys are relatively resistant to aflatoxicosis. This difference is largely due to functional hepatic GSTA-mediated detoxification activity of the bioactive electrophilic AFBO intermediate that is completely lacking in domesticated birds [20]. The present data indicates other pathways may also account for difference in AFB 1 susceptibility, such as cellular regulation, modulation of apoptosis, inflammatory responses, and other pathways relevant to AFB 1 pathogenesis. The liver is the principal organ of AFB 1 bioactivation and detoxification [6,21,22,24,45]. In turkeys, AFB 1 causes reduced feed intake, weight gain, and immunological function in a dose-dependent fashion [46,47]. Dietary exposure in poultry causes lipid accumulation, resulting in hepatomegaly and increases in liver:body weight ratios [48][49][50]. During the 14 day exposure, decreased relative liver mass initially occurred in both EW and DT consistent with that observed in chickens [49] and wild turkeys [9].
Numerous significant DEGs occurring in the livers of AFB 1 -treated birds have potential roles in lipid metabolism or accumulation. AFB 1 is known to alter lipid metabolism and increase lipid content resulting in pale or yellowed pigmentation [46]. Dietary AFB 1 primarily down regulated several hepatic apolipoprotein genes (cofactors in lipid binding and transport) in the turkey, and dis-regulation of genes, such as ANGTPL3, would have direct effects on lipids. Significant up regulation of ANGTPL3 was observed for both EW and DT birds treated with AFB 1 . This would likely stimulate synthesis of plasma triglycerides (TG) via the inhibition of lipoprotein lipase (LPL) activity. In both AFB 1 -treated groups, LPL was significantly down regulated (log 2 FC = −2.905 and −6.032 in EW and DT birds, respectively). LPL functions in the hydrolysis of triglycerides in lipoproteins and is essential to lipid metabolism and storage. Significant down regulation of LPL was also observed in our previous analyses of AFB 1 -treated domesticated Orlopp turkeys [26] and decreased expression of LPL occurs in AFB 1 -treated chickens [50].
As expected, the significant hepatic DEGs included the Phase I and II detoxifying enzymes that we have shown are relevant to AFB 1 exposure in turkeys (Table 6). Previous studies have demonstrated efficient epoxidation by hepatic turkey cytochromes CYP1A5 and CYP3A37 [16]. At environmentally-relevant hepatic concentrations (<50 uM) CYP1A5 bioactivates the majority (~98%) of AFB 1 [17,21], whereas CYP3A37 predominates at much higher substrate concentrations unlikely to be achieved in the livers of exposed animals [16]. Based on RNA-seq, it is clear that dietary AFB 1 significantly down regulated CYP1A5 in both EW and DT birds, but more significantly so in DT. This result is at odds with our earlier findings in another strain of DT (Orlopp) where almost no expression change was observed for CYP1A5 and CYP3A37, and where none of the transcripts associated with CYP genes had significant DE as a result of AFB 1 treatment [26]. Significant down regulation of CYP1A5 in response to AFB 1 was also observed in ducks, another avian species with high AFB 1 susceptibility [51]. Several other P450 genes in addition to 1A5 and 3A37 had significant DE in the present study (Table 6), including both CYP2W1 and CYP2K1. Interestingly, these genes have been shown in other species to activate AFB 1 into cytotoxic products [52,53]. We have found CYP2W1-like transcripts to have significant DE in DT embryos challenged with AFB 1 [28]. Down regulation of CYP1A5 in both EW and DT birds could affect their overall ability to bioactivate AFB 1 . However, as this expression change was seen in both bird types, it does not account for the differences seen in AFB 1 susceptibility [16].
Expression of GSTs with affinity toward AFBO is a known predictor of relative AFB 1 resistance [20]. Constitutive expression of GSTA3, the ortholog to the putative AFB 1 -protective GSTA3 isoform in mice [18] was significantly higher in EW than in DT birds. Dietary AFB 1 caused significant down-regulation of hepatic α-class GSTs, with the exception of GSTA3, where increased expression of this isoform was observed in the AFB 1 -treated DT group. This pattern was also observed in the qRT experiments of other wild (RGW) and domesticated (BB) turkeys. A similar pattern of GSTA3 expression in response to AFB 1 was also observed in turkey embryos early after exposure, where small increases were observed in DT [28] and in ducks [51].
Expression of GSTA3 mRNA in turkeys is not correlated with AFB 1 sensitivity in that domesticated birds lack hepatic GST-mediated AFBO conjugating activity [19], despite expression of GSTA3. Hepatic cytosols isolated from wild turkeys possess functional AFBO-trapping GSTs [20]. While hepatic GSTs in DT lack detoxification activity, with or without AFB 1 treatment, increased GSTA3 expression in DT in response to AFB 1 may reflect a greater inflammatory response or perhaps an indicator of hepatocyte injury. Although GSTs are toxicologically important for their role in "trapping" electrophilic intermediates by conjugating with the nucleophilic GSH, they may also play a role in cell signaling through binding of non-substrate ligands to mediate cell proliferation and cell death [54]. Up regulation of GSTAs may also reflect antioxidant functions as AFB 1 , exposure in poultry can lead to oxidative stress and lipid peroxidation [55,56]. When combined, these results support the hypothesis that the greater ability of wild turkeys to detoxify AFB 1 is related to higher constitutive expression of GSTA3, coupled with an inherited (genetic) difference in functional expression in domesticated birds. Expression in these CYP and GSTAs suggests that the physiological response to AFB 1 is mediated through genes not experimentally linked in the turkey to AFB 1 metabolism.
Up regulation of transcription factors and metabolic inhibitors characterized the shared response to AFB 1 . Taken together, these are genes that comprise the molecular mechanisms underlying aflatoxicosis. Recurrent themes amongst the many DEGs of AFB 1 -treated birds are linked by functional analysis to inflammation, apoptosis, the cell cycle (cancer), or lipid regulation, suggesting common underlying regulation. For example, recent studies of AFB 1 -induced hepatocellular carcinoma have examined regulatory ncRNAs (miRNA and lncRNA) [57,58]. Studies in the rat, another AFB 1 -susceptible species, have found coincident DE of transcripts that are related to these same functions and specific lncRNAs in hepatocellular carcinomas [59,60]. Our study of miRNA expression in the same turkey liver tissues used in the present study is currently underway (Coulombe, unpublished).
Transcriptome analysis not only includes genes responding to the presence of AFB 1 , but also reveals genes dis-regulated as a response to toxic insult. Significant up regulation was seen for several vasoactive peptides, including, neuropeptide Y (NPY), somatostatin (SST), substance P (tachykinin, TAC1), and vasoactive intestinal peptide (VIP), suggesting altered sinusoidal blood flow with AFB 1 treatment. Also, affected were extracellular matrix proteins including glycoproteins (e.g., HAPLN1 and HAPLN3), protein receptors (KERA, LAMB3, LUM, LRRN2, and LRRN3), proteinase inhibitors (TIMP4), signaling molecules (SFRP1, Wnt6 and 7a), and structural proteins (COL10A1, FRAS1). Expression of the majority of the ADAM metallopeptidases was altered in AFB 1 treatment (Table S2). Some of these proteases are thought to be involved in regulating matrix degradation [61]. Unique response in the EW birds was seen in genes that negatively regulate cellular processes, components of the extracellular matrix and accumulation of coagulation factors. DT birds showed greater up regulation of genes responding to inflammation, which was likely due to the reduced ability to detoxify AFBO. Dis-regulation of extracellular matrix proteins is a resulting effect of chronic liver injury [32]. Aflatoxin inhibits cell-mediated immunity in domestic poults [47,62] with the suppression of lymophoblastogenesis [9], T-helper, or cytotoxic T-cell activity [63].
Multiple genes involved in pathways of coagulation (FGA, FGB, FGG, F9, HRG, SERPINA10, and SERPINC1) were expressed at higher levels in EW as compared to DT, where they were among the genes with the highest negative fold change. Lower expression of coagulation factors was also seen in livers of domesticated turkey embryos after just 5 days of exposure to AFB 1 [28]. AFB 1 has been shown to increase blood clotting times in poultry [64,65] and activities of coagulation factors, such as F9, were reduced by dietary AFB 1 in chickens [50,65]. Effects on hemostasis are more dramatic in turkeys than chickens [66]. In comparison, only small non-significant increases in prothrombin times were seen in wild turkeys exposed to AFB 1 [9], which is consistent with the gene expression patterns observed in the liver transcriptomes.
In a previous comparison of EW and DT after in ovo exposure [28], we used RNA-seq to examine gene expression responses to AFB 1 in the embryonic hepatic transcriptomes and identified gene expression effects dependent on exposure time and turkey type. Most notable in turkey embryos was the more rapid response of the DT, which was likely due to their lack of GST activity towards the AFBO-epoxide. The present study was designed to contrast gene expression responses in the hepatic transcriptome of growing domesticated and wild turkeys during AFB 1 exposure. In conclusion, our findings emphasize the differential response of these genetically distinct birds, demonstrating significant differences in expression of Phase I and Phase II genes and in genes important in cellular regulation, modulation of apoptosis, and inflammatory responses. The molecular basis for the differences in AFB 1 detoxification observed between EW and DT birds, and the mechanism of GSTA silencing in DT remain under investigation.

Materials and Methods
This study used two turkey subspecies previously demonstrated to vary in AFB 1 -detoxifying GST activity. Eggs from domesticated (DT = Nicholas) and a wild subspecies (Eastern wild = EW, Meleagris gallopavo silvestris) were obtained from Privett Hatchery (Portales, NM, USA) and hatched at Utah State University. Birds were sexed by PCR [67]. Male turkey poults were maintained on an ad libitum standard grow-up soy-based diet and acclimated to the facility for two weeks. At the end of this period, males from each line (n = 8 for EW and n = 10 for DT) were equally assigned to one of two treatment groups and subjected to a short-term AFB 1 -treatment protocol [21,68]. For the AFB 1 treatment, the diet of challenge birds was amended beginning on day 15 with 320 ppb AFB 1 (Sigma-Aldrich, Inc., St. Louis, MO, USA) that continued for 14 days. Control birds continued on the standard diet with AFB 1 levels below detection limits (<10 ppb), based on testing of 50 g of feed extracted and cleaned using Mycosep 112 AflaZON cleanup columns (Romer Labs., Union, MO, USA), and examined by HPLC. Birds were weighed three times per week and feed and water availability checked daily. At the conclusion of the 14 day challenge period, birds were sacrificed by CO 2 asphyxiation and blood collected by cardiac puncture for DNA and serological analysis. Livers were removed, examined, weighed, sampled, and fixed in neutral buffered formalin for histological examination. Portions of the liver tissues infused with RNAlater (ThermoFisher Scientific, Waltham, MA, USA) for RNA-Seq analysis. All of the procedures were under the authority and institutional approval of Utah State University's Animal Use and Care Committee. Ethical approval code: 2670, Date of approval: 26 September 2016.

RNA Isolation and Sequencing
Total RNA was isolated from each sample by TRIzol extraction (Ambion, Inc., Foster City, CA), DNase-treated (Turbo DNA-freeTM Kit, Ambion, Inc.), and stored at −80 • C. Initial RNA concentration and quality was assessed by spectrophotometry (Nanodrop 8000). RNA samples were submitted for library preparation and sequencing at the University of Minnesota Genomics Center (UMGC). Replicate samples were sequenced from each treatment group (n = 4). Each sample was quantified by RiboGreen assay (Invitrogen Corp., Carlsbad, CA, USA) and RNA integrity confirmed on a 2100 Bioanalyzer (Aligent Technologies, Santa Clara, CA, USA). Each sample had clear 18S and 28S peak separation on the electropherograms and an average RNA Integrity Number (RIN) of 6.3. Indexed libraries (n = 16) were constructed with 1 µg of total RNA/sample with the TruSeq RNA Sample Preparation Kit version 2 (Illumina, Inc., San Diego, CA, USA) and size selected for approximately 200 bp inserts. Libraries were multiplexed, pooled, and sequenced over two lanes on the HiSeq 2000 using v3 chemistry (Illumina, Inc.) to produce 101-bp paired-end reads. Data are deposited in the NCBI Gene Expression Omnibus (GEO) repository as part of SRA BioProject 346253.

RNA-seq Data Analyses
Sequence reads were groomed (Trimmomatic, [69]) and quality checked (FastQC, [70]) prior to read mapping (Bowtie v2.2.4.0) on the turkey genome (UMD 5.0, NCBI Annotation 101). Read counts were normalized in CLC Genomics Workbench (CLCGWB v. 8.0.2, CLC Bio, Aarhus, Denmark) by dividing the total read counts by the group sample sum with the results being expressed as reads per 12.2 M. Hierarchical clustering of samples was performed (based on Euclidean sample distances with single linkage) in CLCGWB. Principal component analysis (PCA), Volcano plots, and Venn diagrams were used to visualize expression data and the results of significance testing. Empirical analysis of differential gene expression and ANOVA were performed in CLCGWB on EdgeR-normalized read counts. Pair-wise comparisons between treatment groups were made in CLCGWB following the standard workflow Wald test with multi-comparison p-values < 0.05 being considered as significant (Bonferroni and FDR corrected). In each pair-wise comparison, significant DE genes were used to investigate affected gene pathways using Ingenuity Pathway Analysis (IPA, Ingenuity Systems, Redwood City, CA, USA). Gene enrichment tests were performed using the PANTHER Overrepresentation Test (GO Consortium release 20150430, [38]).

Quantitative Real-Time PCR
To more broadly examine the expression profile response of GSTA3 to dietary AFB 1 , quantitative real-time PCR (qRT-PCR) was performed on both domesticated and wild turkey liver samples. Samples included AFB 1 -treated and control animals (six per group) from the domesticated Nicholas turkey (DT) and Eastern Wild (EW) experiment, plus AFB 1 -treated and control animals (six per group) of a parallel experiment that included domesticated Broad Breasted White (BB), and birds of the Rio Grande subspecies (RGW, M. g. intermedia) of wild turkey. Four of the six samples for the DT and EW groups were in common with the RNA-seq study.
Briefly, cDNA was synthesized from DNase-treated liver mRNA (TRIzol extracted) using Invitrogen Super Script IV First-strand synthesis kit (Invitrogen, Carlsbad, CA, USA). Expression analysis of gene-specific amplicons was performed with the iTaq Universal SYBR Green Supermix (BioRad, Hercules, CA, SA) with the CFX96 touch real time detection system (BioRad, Hercules, CA, USA). Primers were designed within Primer3 software (http://www.ncbi.nlm.nih.gov/tools/primerblast) from accessioned genomic DNA sequence (NM_001303157.1) to span an exon/exon junction and at least one intron in the amplicon. RefFinder software was utilized to determine the most stable reference gene. Several normalizing genes were tested for uniformity between treatments and lines and RNA polymerase II subunit D (POLR2D, XM_003208947) was found to have the highest stability value (0.848). Target gene reactions were conducted in triplicate, and POLR2D, no template and gDNA controls were run in duplicate. All of the reactions included a disassociation curve to confirm a single product and to preclude the possibility of dimers amplifying. Expression in each RNA sample was normalized first to the control gene POLR2D. Results were interpreted using the Double Delta Ct Analysis (∆∆Ct, [71]) and a comparative Ct approach. Expression analysis was performed within the Biorad CFX Maestro software package following the standard ∆∆Ct workflow.
Supplementary Materials: The following are available online at www.mdpi.com/2072-6651/10/1/42/s1. Figure S1: Hierarchical clustering of samples based on Euclidean distance reiterated relationships shown by PCA. Figure S2: Distribution of genes expressed in turkey liver by treatment group. Figure S3: Differential fold change in DEGs shared by Eastern wild (EW) and domesticated (DT) birds exposed to AFB 1 in comparison to controls. Table S1. Mean quality-trimmed RNA-seq read counts for genes with mapped reads in the livers of wild and domesticated turkeys. Table S2. Summary of pairwise differential gene expression analysis of liver transcriptomes. Table S3. Fifty genes showing the greatest differential expression in each pairwise comparison of treatment groups. Table S4. Summary of PANTHER Overrepresentation test of the 3380 differentially expressed genes shared in the AFB 1 -treated birds as compared to controls. Table S5. Significant differentially expressed genes (FDR p-values < 0.05 and |log 2 FC| > 2.0) identified in each pairwise comparison of genetic groups (Eastern wild and domesticated turkey).