If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Variants in the 3′ untranslated region of the ovine acetyl-coenzyme A acyltransferase 2 gene are associated with dairy traits and exhibit differential allelic expression
The Roslin Institute and Royal (Dick) School of Veterinary Studies, University of Edinburgh, EH25 9RG Midlothian, United KingdomAnimal and Veterinary Sciences, Scotland's Rural College, EH25 9RG, Midlothian, United Kingdom
The acetyl-CoA acyltransferase 2 (ACAA2) gene encodes an enzyme of the thiolase family that is involved in mitochondrial fatty acid elongation and degradation by catalyzing the last step of the respective β-oxidation pathway. The increased energy needs for gluconeogenesis and triglyceride synthesis during lactation are met primarily by increased fatty acid oxidation. Therefore, the ACAA2 enzyme plays an important role in the supply of energy and carbon substrates for lactation and may thus affect milk production traits. This study investigated the association of the ACAA2 gene with important sheep traits and the putative functional involvement of this gene in dairy traits. A single nucleotide substitution, a T to C transition located in the 3′ untranslated region of the ACAA2 gene, was used in mixed model association analysis with milk yield, milk protein yield and percentage, milk fat yield and percentage, and litter size at birth. The single nucleotide polymorphism was significantly associated with total lactation production and milk protein percentage, with respective additive effects of 6.81 ± 2.95 kg and −0.05 ± 0.02%. Additionally, a significant dominance effect of 0.46 ± 0.21 kg was detected for milk fat yield. Homozygous TT and heterozygous CT animals exhibited higher milk yield compared with homozygous CC animals, whereas the latter exhibited increased milk protein percentage. Expression analysis from age-, lactation-, and parity-matched female sheep showed that mRNA expression of the ACAA2 gene from TT animals was 2.8- and 11.8-fold higher in liver and mammary gland, respectively. In addition, by developing an allelic expression imbalance assay, it was estimated that the T allele was expressed at an average of 18% more compared with the C allele in the udder of randomly selected ewes. We demonstrated for the first time that the variants in the 3′ untranslated region of the ovine ACAA2 gene are differentially expressed in homozygous ewes of each allele and exhibit allelic expression imbalance within heterozygotes in a tissue-specific manner, supporting the existence of cis-regulatory DNA variation in the ovine ACAA2 gene. This is the first study reporting differential allelic imbalance expression of a candidate gene associated with milk production traits in dairy sheep.
); therefore, the improvement of milk production is the most important breeding objective. In Mediterranean countries, however, most of the ovine milk produced is used for the production of cheese commercialized as products of protected designation of origin (PDO) and other quality labels (
). Although traditional breeding programs have achieved appreciable genetic gains mainly for milk yield, application of selection schemes assisted by molecular information could expedite improvement (
). Moreover, marker-assisted selection could be of special interest for dairy sheep due to the great regional diversity of breeds, funding limitations, organizational difficulties (
To date, few reports exist of genome-wide association studies and genome scans based on linkage mapping that detect QTL and quantitative trait mutations for ovine dairy traits (reviewed by
Identification of quantitative trait loci underlying milk traits in spanish dairy sheep using linkage plus combined linkage disequilibrium and linkage analysis approaches.
detected suggestive QTL for milk, fat, and protein yields, mapped in a region of the ovine chromosome 23 harboring the acetyl-CoA acyltransferase 2 (ACAA2) gene. The ACAA2 gene encodes an enzyme of the thiolase family, also known as 3-oxoacyl-CoA thiolase or mitochondrial 3-ketoacyl-CoA thiolase. The ACAA2 enzyme catalyzes the last step in mitochondrial fatty acid β-oxidation, thus playing a central role in the supply of energy for the animal (
and its functional role in lipid metabolism, it was regarded as a putative functional and positional candidate gene that may affect milk yield and composition.
Genes encoding enzymes of the thiolase family have been correlated with production traits in other livestock species. Single nucleotide polymorphisms detected in the porcine ACAA2 gene are reported to be associated with daily weight gain and loin muscle area (
I: Segregation, Cloning, SNP Detection and Its Association with Economic Traits of 5 Bovine Genes. II: Quantitative Trait Loci Analysis of Meat Quality Traits in Pigs. PhD thesis. Chinese Academy of Agricultural Sciences,
Beijing, China2008
). An important paralog of the ACAA2 gene, the acetyl-CoA acetyltransferase 2 (ACAT2) gene, has been associated with production and fertility traits (milk protein content, productive life, and conception and pregnancy rates) in Holstein cattle (
), whereas SNP within the swine ACAT2 gene were suggested to influence the metabolic functions of the corresponding enzyme and thus may affect growth performance (
Our previous study showed that the entire mRNA (coding and untranslated regions) of the ovine ACAA2 gene is monomorphic in Chios sheep, one of the most productive and extensively used breeds in Greece and Cyprus (
). The SNP was significantly associated with milk yield at first lactation and across first to third lactations in Chios sheep. Animals from a closed nucleus research flock at the Agricultural Research Institute of Cyprus (Nicosia) carrying the g.2982TT or g.2982CT genotype had significantly higher milk yield than those with the g.2982CC genotype, and the g.2982T>C SNP explained 10% of the additive genetic variance for milk yield when data up to third lactation from a single flock were analyzed (
) have shown that alleles of nonimprinted genes are not expressed equally at the mRNA level in heterozygous animals, a phenomenon called allelic expression imbalance (AEI). Allelic expression imbalance is the outcome of the presence of at least 1 cis regulatory element in the regulatory sequences of a gene (
); therefore, AEI is one of the possible mechanisms underlying the effect of causative genetic variations that are not located on the translated region of a gene.
The objective of the current study was to provide novel insights into the association of the ovine ACAA2 gene with important sheep traits. First we performed an association analysis of the previously identified g.2982T>C SNP with milk yield, fat and protein contents, fat and protein yields, and litter size in a population of Chios sheep from all available farms keeping production records in Cyprus. Upon confirmation of the association of the SNP with total milk production and additional detection of association with milk protein percentage, we tested the hypothesis that the ACAA2 is a putative functional gene affecting dairy traits by comparing the expression of the gene in different genotypes and by developing an assay to test the presence of AEI.
MATERIALS AND METHODS
Animals and Phenotypic Data
Data were collected from 742 purebred Chios ewes from 5 farms. The farms included 1 governmental farm (Orites, Cyprus) and the only 4 commercial farms in Cyprus keeping phenotypic records of purebred Chios sheep according to the regulations of the International Committee for Animal Recording (
Section 2 ICAR rules, standards and guidelines for dairy production recording. Approved by the General Assembly held in Berlin, Germany, May 2014. ICAR (International Committee for Animal Recording),
Rome, Italy2014
For all animals, individual records included month of lambing, year of lambing, lactation number, and age of lambing. Phenotypic data were obtained for lambing years between 2009 and 2016 and included 1,514 individual records of 742 milking ewes; this data included total lactation milk yield, 1,242 observations for litter size at birth, 1,203 measurements of milk fat content and yield, and 615 measurements of milk protein content and yield, respectively. Total lactation yield was calculated for each animal with the Fleischmann method with monthly tests on actual yields (sum of a.m. and p.m. records,
Section 2 ICAR rules, standards and guidelines for dairy production recording. Approved by the General Assembly held in Berlin, Germany, May 2014. ICAR (International Committee for Animal Recording),
Rome, Italy2014
). Milk samples were obtained for fat and protein analysis by combined thermo-optical procedures (LactoStar 3510, Funke Gerber, Berlin, Germany), previously calibrated for protein with the Lowry protein assay and for fat with the method 989.05 (
Whole-blood samples were obtained from all 742 ewes for DNA extraction and genotyping of the g.2982T>C SNP of the ACAA2 gene. Genomic DNA was isolated from all samples using the Genomic DNA Blood kit (Macherey-Nagel, Düren, Germany), according to the manufacturer's instructions, and DNA quality and quantity were estimated by UV absorption at 260 and 280 nm. Genotyping of the Chios ewes was performed by a cost effective direct DNA sequencing protocol (
Each trait was analyzed separately with the mixed linear model
[1]
where Y = lactation milk yield, fat yield/percentage, protein yield/percentage, or litter size record n of animal m; µ = overall population mean for the trait; F = fixed effect of flock i (1–5); YS = fixed effect of year (2009–2016) by season (1–2) of lambing interaction j; L = fixed effect of number of lactation k (1–4); b1 = linear regression on age at lambing (age); b2 = linear regression on lactation duration (dur), yield traits only; G = fixed effect of genotype l in the g.2982T>C locus (1–3; CC, CT, TT); A = random effect of animal m; and e = random residual effect. All lactation records were analyzed for all traits. In a separate analysis, only first-lactation records were considered for milk yield. In the latter case, the fixed effect of number of lactation and the random animal effect were removed from the model.
In all cases, predicted trait values for each SNP genotype and respective standard errors were derived; these values were reflective of the marginal genotypic effect on each trait adjusted for all other effects fitted in the model. The predicted trait values were used to estimate additive and dominance SNP effects on traits and the proportion of phenotypic variance for each trait accounted for by the SNP locus. The equations used were
where a is the additive effect; d is the dominance effect; VP is the percentage of phenotypic variance due to SNP; TT, CC, and CT are the predicted trait values for each genotype class; and p and q are the allele frequencies at the SNP locus. Variance components were estimated with model [1] after excluding the genotype effect. All statistical analyses were conducted with the ASReml3 software (
As multiple traits were analyzed, we assessed the number of independent tests to adjust for multiple testing. For this reason, we conducted multivariate analyses based on model [1] after excluding the genotype effect to assess the correlation among traits. Very high correlations were found among the 3 production traits (milk, fat, and protein yield; 0.89–0.99), and between the 2 milk concentration traits (fat and protein percentage; 0.75). Correlations between the 2 groups of traits were near zero. Correlations of these 2 groups with litter size ranged from 0.05 to 0.19. Consequently, we regarded 3 separate trait groups (production, concentration, and litter size) that corresponded to 3 distinct, independent hypothesis tests. This was confirmed with a principal component analysis of the studied traits (Supplementary Figure S1; https://doi.org/10.3168/jds.2016-12326). Subsequently, a Bonferroni adjustment for multiple tests was implemented based on the Holm-Bonferroni method (
). This method works sequentially by testing first the lowest nominal P-value against the threshold value of 0.05/n, where n is the number of independent hypotheses tests (3 in the present study). If this hypothesis is rejected, the next lowest nominal P-value is compared with 0.05/(n − 1) and so on, until a hypothesis is not rejected (
Animals and Tissue Sampling for Expression Analysis
Blood samples were collected from 161 first-parity ewes from the biggest commercial farm in Cyprus, allowing the selection of age-, parity-, and lactation stage-matched animals. Based on the methodology described in our previous work (
), we identified animals with the genotypes g.2982TT, g.2982CT, and g.2982CC, and chose 3 ewes from each genotype for subsequent biopsies and RNA extraction. All 9 selected ewes were first parity, 15 mo old, and had given birth to 2 lambs each.
Mammary and liver biopsies were obtained under anesthesia by a professional veterinarian at 42 ± 2 d after lambing, a week after weaning. Liver parenchyma was sampled via puncture biopsy using ultrasonography for the selection of the appropriate biopsy site and to avoid large vessels and other organs. Udder biopsies were taken from either the left or the right rear gland. The biopsy site was carefully selected to avoid large subcutaneous blood vessels. Preparation of the site involved shaving and washing with dilute betadine solution followed by sanitizing with ethanol (70%). Ewes were given intravenous xylazine before anesthetizing the biopsy site by subcutaneous injection of lidocaine hydrochloride. An incision was made (∼0.5–1.0 cm) on the outside of the quarter using a scalpel blade (size 22). A Bard Magnum core biopsy instrument (Bard Peripheral Vascular Inc., Tempe, AZ) with a Bard Magnum core tissue biopsy needle (MN1210, 12 gauge × 10 cm) was used.
In addition, liver and udder ex vivo tissue samples were obtained from 16 randomly selected Chios ewes from an abattoir directly after slaughtering under sterile conditions. Although tissue samples from the abattoir were obtained at different and unknown parities and stages of lactation, they were taken to increase the number of heterozygous samples for allelic expression imbalance analysis without increasing the veterinary costs required for biopsies. All samples were snap-frozen in liquid N2 immediately after the tissue was obtained under sterile conditions and stored at −80°C until RNA extraction.
All animal sampling and handling procedures in our study were carried out in strict accordance with the national legislation (
) and no animals were euthanized for the purposes of this study.
DNA, Total RNA Extraction, and Reverse Transcription
Frozen tissue biopsies and ex vivo samples were homogenized using a mortar and pestle constantly covered under liquid N2 and aliquoted in tubes containing ∼20 mg of tissue biopsy each. The DNA was extracted using the Genomic DNA Nucleospin Tissue kit (Macherey-Nagel), whereas ∼20 mg of frozen homogenized tissue was subjected to RNA extraction using the RNA isolation Nucleospin RNA kit (Macherey-Nagel) according to manufacturer instructions (http://www.mn-net.com/Portals/8/attachments/Redakteure_Bio/Protocols/RNA%20and%20mRNA/UM_TotalRNA.pdf). Before reverse transcription, extracted RNA was incubated with 2 U/μL of rDNAase (Macherey-Nagel) for 10 min at 37°C to eliminate putative contaminating genomic DNA and then purified using the RNA clean up Nucleospin RNA Clean up XS kit (Macherey-Nagel). The concentration and quality of DNA and RNA for all samples was measured using a Nanodrop 1000 UV/VIS spectrophotometer (Thermo Fisher Scientific, Waltham, MA). The RNA integrity was assessed by electrophoretic analysis of the 28S and 18S rRNA subunits.
Complementary DNA was synthesized from 0.5 µg of RNA, using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Thermo Fisher Scientific), using random hexamer primers following the manufacturer's recommendations in a final volume of 20 μL. The lack of genomic DNA (gDNA) contamination was verified by a PCR amplification using primers amplifying intron 2 of the prolactin gene (PRL), as described by
, and by including a control without reverse transcriptase during cDNA synthesis of each sample. These controls were used for conventional PCR amplification of the ACAA2 exon 10, as we have previously described (
), and for subsequent real-time PCR and verified lack of genomic contamination.
ACAA2 mRNA Expression
The ACAA2 relative expression in liver and udder was assessed by quantitative reverse transcriptase PCR (qRT-PCR). The ACAA2 qRT-PCR was performed using primers ACAA2_SNP forward: AACCAGGTGACCTTCAGAGC and ACAA2_SNP reverse: AATGGGTCTTTGCTTCACCA, designed and synthesized by VBC-Biotech Service (Vienna, Austria), and expression levels were normalized by the expression of reference genes, such as ovine β-actin (ACTB) for liver samples and ovine β2 microglobulin (ovB2M) for udder samples. We utilized primers for ovACTB forward (GCAAAGACCTCTACGCCAAC) and ovACTB reverse (TGATCTTGATCTTCATCGTGCT) (
) and for ovB2M forward (CTGTCGCTGTCTGGACTGG) and ovB2M reverse (TTTCCATCTTCTGGCGGGTG; designed using the NCBI/Primer-BLAST tool; https://www.ncbi.nlm.nih.gov/tools/primer-blast/). For each gene and cDNA sample qRT-PCR reactions were performed in triplicates. Reactions of no-template negative control and 5-point 1:5 serial dilutions of reference cDNA, to obtain the reaction standard curves and calculate amplification efficiencies, were also performed in triplicates. Reactions for ACAA2 expression were performed in a final volume of 20 μL containing 1× KAPA SYBR fast qPCR Master Mix (Kapa Biosystems, Wilmington, MA), 100 nM of each of ACAA2_SNP forward and reverse primers, 1× ROX (Kapa Biosystems, Wilmington, MA) low, and cDNA. Reactions for ACTB and B2M were performed in a final volume of 10 μL containing 1× KAPA SYBR fast qRT-PCR Master Mix (Kapa Biosystems), 200 nM of each of the ovACTB or ovB2M forward and reverse primers, 1× ROX low, and cDNA. Cycling conditions included a 3-min initial denaturation step at 95°C followed by 40 cycles of (1) denaturation for 3 s at 95°C and primer annealing and template extension for 30 s at 60°C for ACAA2, (2) denaturation for 3 s at 95°C, primer annealing at 62°C for 20 s and template extension at 72°C for 30 s for ACTB, and (3) denaturation for 3 s at 95°C, primer annealing at 64°C for 15 s and template extension at 72°C for 30 s for B2M. For all 3 genes, a melting curve analysis was performed following the amplification cycles consisting of 15 s at 95°C, 30 s at the respective annealing temperature for each gene, and continuous heating and data collection up to 95°C at a rate of 1% temperature increase per 30 s to evaluate specificity of the amplification. Reactions were optimized for amplification efficiency to be between 85 and 95% and linear standard curve fit (R2) to be greater than 0.990 for all genes. Raw data were analyzed with the 7500 software v2.3 (Applied Biosystems), and mean Ct values and PCR reaction efficiencies were exported. The normalized expression (nQ) of ACAA2 was calculated according to the delta-Ct method (
); that is, nQ = EACAA2(minACAA2_Ct − sampleACAA2_Ct)/Eref(minref_Ct − sampleref_Ct), where E is the PCR efficiency for a given gene (values are greater than 1 and 2 corresponds to ideal doubling of templates per cycle, 100% = 2) estimated based on standard cDNA dilution series reactions, mingene_Ct is the minimum Ct value among all samples for a given gene, and samplegene_Ct is the Ct value for a given gene and sample. To express ACAA2 nQ as fold change relative to nQ in CC genotype, nQ for each sample was divided by mean nQ in CC samples. Data were then analyzed using the mixed models procedure of SAS software (
) with genotype as the fixed effect. Mean comparisons were performed with Tukey's adjustment with significance level set at 0.05.
AEI Analysis
Sequencing of the ACAA2 gene from gDNA and cDNA of each sample was initially used as a semiquantitative approach to detect differential expression of the ACAA2 alleles. The PCR products amplified from heterozygous individuals were sequenced and data were analyzed using Applied Biosystems Sequencing Analysis. The SNP was identified by visually inspecting each base in sequencing traces. Allelic variation was estimated by comparing the proportions of the peak heights of the 2 alternative alleles of the SNP.
Subsequently, the detection of AEI was based on quantitative analysis of mRNA transcripts using a TaqMan probe qRT-PCR assay (
) to detect deviations from the null hypothesis expecting equimolar ratio between the 2 alleles in heterozygous samples. TaqMan qPCR reactions were performed using the ACAA2 primers described above and probes SNP_2982_C (5'FAM/3′BHQ1 labeled) CACGGAGAGTGACCAGGTTTGGGTAAGG and SNP_2982_T (5'-JOE/3′BHQ1 labeled) CACGGAGAGTGACTAGGTTTGGGTAAGG, designed and synthesized by VBC-Biotech Service, in 10-μL reaction volumes containing 1× Type-it Fast SNP Probe PCR Master Mix (Qiagen, Hilden, Germany), 400 nM of each PCR primer and 200 nM of each TaqMan probe, and liver or udder cDNA or gDNA from ACAA2 heterozygous ewes. A 7500 Real Time PCR system (Applied Biosystems) was used with a cycling profile that included a 5-min denaturation step at 95°C and 40 cycles of 15 s at 95°C and 30 s at 62°C.
Pooled gDNA from 3 ACAA2 heterozygous animals (genotype confirmed by several replications of genotyping reactions) was used as standard gDNA, and 5 points of 1:4 serial dilutions were used to construct the standard curve for relative ACAA2 T and C allele quantification (Figure 1A-C). The T:C allele copy ratio was expected to equal 1 in all dilutions.
Figure 1Allelic expression imbalance TaqMan assay for the quantification of the ACAA2 allele T:C transcription ratio. (A) Amplification plots of TaqMan probe PCR from ACAA2 TT (top), TC (middle), and CC (bottom) individuals. (B) Amplification plot of serial dilution of standard genomic DNA pooled from 3 ACAA2 TC ewes. (C) Log10 quantity and threshold cycle (Ct) relationship graph for T and C alleles. Red dots represent standard genomic DNA dilution series reactions. Linear standard curves for T and C alleles (targets), standard curve slope, Y-intersection, R2, and efficiency values are shown.
The quantitation relative to standard curve method was employed using the homonymous option of the ABI 7500v2.3 software (Applied Biosystems, Foster City, CA) to extract the value for a normalized quantity of T allele at each qPCR reaction, having set the C allele as the endogenous control. The normalized quantity of the T allele was computed as the quantity of T allele divided by the quantity of C allele, and thus corresponds to the ratio of T:C ACAA2 allele quantity. The T and C allele quantities in each qRT-PCR reaction were computed from the respective threshold cycle values and the linear prediction function for each of the ACAA2 alleles, computed by the software based on standard gDNA dilution series. The quantity of T and C alleles in the standard gDNA was set at the same arbitrary value; thus, the T:C ratio equaled 1 at all points of the dilution series of the standard gDNA. The T:C ratio for each sample (i.e., the mean normalized quantity of the 3 replicates performed with each cDNA or gDNA) was divided by the average normalized quantity of the T allele in all the gDNA samples to obtain the corrected T:C ratio for each cDNA sample. Average corrected T:C ratio in the gDNA samples equaled 1 and ranged between 0.90 and 1.09. Standard gDNA dilution series reactions were performed in each qRT-PCR run. In all runs, efficiency values for both T and C allele amplification were very similar (102.01% for T and 101.23% for C in run 1, and 102.66% for T and 102.74% for C in run 2). Linear standard curve r2 values were 0.997 for T and 0.996 for C in run 1, and 0.998 for T and 0.998 for C in run 2.
Statistical analysis of the T:C ratios in cDNA or gDNA was performed using the mixed models procedure of SAS software (
) with tissue (liver or udder) as the fixed effect. Means are presented as ordinary means and standard error. Least squares means pairwise comparisons were performed with Tukey's adjustment with significance level set at 0.05.
RESULTS
Association of ACAA2 with Sheep Traits
Allelic frequencies in the g.2982C/T SNP locus were 0.54 for the T allele and 0.46 for the C allele; genotypic frequencies were 0.27, 0.54 and 0.19 for TT, CT and CC, respectively. Genotypic frequencies were found to deviate from the Hardy-Weinberg equilibrium (P = 0.019).
Marginal predicted means for the 3 genotype classes for all production traits are presented in Table 1. After applying Holm-Bonferroni correction, the SNP was significantly associated with milk yield at first lactation (P < 0.01) and all lactations (P < 0.025), with respective additive effects of 10.61 ± 3.56 and 6.81 ± 2.95 kg and respective positive dominance effects of 13.02 ± 4.26 and 8.67 ± 3.53 kg (Table 2). Significant differences were found between the CC and CT genotype pairs (P < 0.01) and between the CC and TT classes (P < 0.025) at the g.2982T>C SNP locus for milk yield (Table 1). These results overall suggest a complete dominance effect at the locus, as heterozygous CT animals exhibit similar predicted mean values compared with homozygous TT animals. Based on the estimated allelic effects and the allele frequencies observed in the sample, it was estimated that the g.2982T>C SNP explained 2.25 and 0.62% of the total phenotypic variance for first lactation and for all lactations, respectively (Table 2). Overall, these results suggest a stronger association in first-lactation milk yield, whereas the effect was reduced in subsequent lactations.
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
a–d Means within a column with different superscripts differ; a,bP < 0.025; a,cP < 0.010; a,dP < 0.001 after Holm-Bonferroni adjustment. The corresponding thresholds for each test are provided in Supplementary Table S1 (https://doi.org/10.3168/jds.2016-12326).
1 Marginal genotype means (±SE) predicted from the mixed model analyses.
After applying Holm-Bonferroni correction, the SNP was also associated with protein percentage (P < 0.01), with a significant additive effect of −0.05 ± 0.02 (Table 1). Pairwise contrasts between the predicted protein percentage values showed significant differences between the genotype classes CC and CT (P < 0.01) and between CC and TT (P < 0.01), with CC animals exhibiting higher protein percentage, whereas differences between CT and TT were not significant (Table 1). It was estimated that the SNP explained 1.68% of the phenotypic variance for protein percentage (Table 1). Additionally, we noted a dominance effect of 0.46 ± 0.21 kg (P = 0.027) and 0.41 ± 0.21 kg (P = 0.048) of the SNP on fat and protein yields, respectively. Homozygous TT and heterozygous CT animals exhibited significantly higher milk fat yield compared with homozygous CC ewes (P < 0.025 and P < 0.001, respectively). No significant associations of the SNP genotype with fat percentage or litter size at birth were found.
mRNA Expression Analysis
The expression of ACAA2 was found to be significantly increased in TT compared with CC udder and liver biopsies from age-matched ewes at the same lactation stage and parity that gave birth to 2 lambs each. In particular, homozygous TT animals showed a 2.8-fold increase in mRNA expression levels compared with homozygous CC animals in liver (P < 0.05), whereas TT ewes exhibited an 11.8-fold increase in mRNA of the ACAA2 gene in the udder compared with CC ewes (P < 0.05; Figure 2A). The ACAA2 mRNA expression of heterozygous animals showed a tendency for reduction compared with homozygous TT animals only in liver (P = 0.07), whereas it did not significantly differ from homozygous CC animals. Heterozygous mRNA expression in the udder exhibited a very high standard deviation (Figure 2A) and did not differ significantly from homozygote expression.
Figure 2Expression of the ACAA2 gene in liver and udder. (A) ACAA2 expression in liver and udder samples (taken 42 ± 2 d after lambing) of CC, CT, and TT first-parity, 15-mo-old ewes that gave birth to 2 lambs each relative to the mean normalized expression in CC ewes. Expression was normalized by β-actin and β-microglobulin expression in liver and udder, respectively. n = 3 for all means except udder CT (n = 2). (B) Mean normalized expression of the ACAA2 gene in ex vivo liver and udder samples from CC, CT, and TT ewes, normalized by corresponding β-actin transcription and relative to mean normalized transcription in the liver of TT individuals. CC: n = 1, CT: n = 10 (udder) n = 9 (liver), TT: n = 4 (udder) n = 5 (liver). Bars represent mean relative normalized transcription and SEM is shown; different letters (a,b) represent significant differences (P < 0.05).
The expression of ACAA2 was also assessed in liver and udder ex vivo samples from 1 CC, 10 CT, and 5 TT ewes (Figure 2B). These ewes were selected at random and were not age-matched nor at the same milking season or lactation stage. Expression of ACAA2 was found to be significantly increased in the liver compared with the udder (P < 0.0001), but no significant differences were observed between CT and TT ewes in either the liver or the udder (P > 0.05). As only 1 ewe from the slaughterhouse was genotyped as CC, it was not possible to make comparisons between homozygotes.
AEI Analysis
To test whether the differences observed between homozygous TT and CC animals was due to AEI between the T and C alleles, the transcription of T and C alleles was quantified in liver and udder of heterozygous individuals. Sequencing of gDNA and cDNA from 12 heterozygous samples (3 from biopsies and 9 from the slaughterhouse) was initially performed to evaluate semiquantitatively whether AEI was present. The allele-specific expression levels were evaluated by comparing the peak heights of C and T in cDNA samples from heterozygous individuals compared with the peak heights of the respective gDNA. Although in all gDNA samples from heterozygotes peak heights of both nucleotides at the SNP position were tangential, suggesting no preferential amplification, peak heights from heterozygous cDNA samples were different, with the T allele exhibiting a higher peak compared with the C allele (Figure 3).
Figure 3Sequencing of heterozygous samples from genomic DNA (gDNA) and cDNA from liver (A) and udder (B). The T/C SNP is shown in light blue shade. The peak heights of the T and C alleles are the same in gDNA, whereas the peak height of the T allele is higher compared with the C allele in cDNA of the same samples. Color version available online.
To further confirm the initial results from the semiquantitative method consistent with AEI in ACAA2 gene, we developed an AEI assay using TaqMan probes to estimate the T:C transcription ratio in cDNA samples corrected by the mean T:C ratio in gDNA samples and ultimately test for deviations from the expected ratio of 1 in the absence of AEI. Results of the assay for the control gDNA reactions showed that observed T:C ratios were, as expected, very close to 1 (average noncorrected ratio was 1.052). The ACAA2 AEI was observed in individual ewes in both the liver and the udder (Figure 4). The average corrected T:C ratio in the liver cDNA was 1.026, not significantly different from the average corrected gDNA ratio (mean = 1, P = 0.67). On the contrary, corrected T:C ratio of transcription in the udder of heterozygous animals was 1.18, significantly increased compared with both gDNA and liver cDNA mean corrected T:C ratios (P = 0.0054 and 0.0134, respectively). Therefore, we observed a moderate and tissue-specific imbalance in the allelic expression of the 2 ACAA2 alleles in favor of allele T in the udder, but not in the liver of heterozygous animals.
Figure 4Allelic expression imbalance of the ACAA2 gene. Graph shows TaqMan assay corrected T:C ratios in the liver and udder of heterozygous ewes. Horizontal lines represent mean corrected T:C allele in the cDNA from liver (mean = 1.03, n = 12) and udder (mean = 1.18, n = 11) and in genomic DNA (gDNA; mean = 1, n = 11) of ACAA2 heterozygous ewes. Open circles represent corrected T:C ratios for individual ewes. *P < 0.05, **P < 0.01.
Our study confirmed the previously observed association of the HM537015:g.2982T>C ACAA2 SNP with milk yield in an extended population of Chios sheep from multiple flocks, and showed that the g.2982T>C SNP of the ACAA2 gene is also associated with milk protein percentage and milk fat yield. Investigation of the expression of the ACAA2 gene from the 3 different genotypes and AEI analysis of heterozygous samples supported the hypothesis that ACAA2 is a functional gene affecting dairy traits. To the best of our knowledge, this is the first study showing differential allelic imbalance expression of a candidate gene associated with milk production traits in dairy sheep.
), the g.2982T>C SNP, was significantly associated with milk yield, with the T allele exhibiting positive additive and dominance effects, mainly attributed to first-lactation production data. Confirmation of a previously detected association using a bigger and different data set strengthens the evidence for the observed association (
, possibly due to the increased variation introduced from the use of multiple flocks managed in different ways, as a substantial fraction of the environmental variance for production traits is attributed to farm (
). In the current analysis, we adjusted for the systematic effect of the flock and, therefore, the estimates are more representative and indicative of the true effect of the gene in the population. The reasons why the effects are more important in the first lactation than in the rest of lactations are not clear at this stage; therefore, further research is needed.
In the present study, additional evidence about the correlation of the ACAA2 gene with important sheep traits is provided, as the g.2982T>C SNP was found to be significantly associated with milk protein percentage. Homozygous CC animals exhibited superior values for protein percentage compared with both homozygous TT and heterozygous CT animals, in contrast to their inferior values estimated for total milk yield (Table 1). This is consistent with the negative genetic correlation between those 2 traits (
). As genotypic classes for protein yield did not significantly differ, the decrease in protein content could be attributed to a dilution effect due to the increase of milk yield (
). However, although the marginal predicted mean for the fat content from CC ewes was higher compared with that of the other genotypes in the present study, the differences were not significant and the investigated SNP was not significantly associated with milk fat percentage (Table 1, Table 2), consistent with the results of
), in the current study standard errors for fat content were twice as much as those for protein content. Therefore, low precision of data for fat content may be a reason for not detecting the association of the studied SNP with fat percentage; however, we found a significant association of the SNP with fat yield (Table 2). In addition, homozygous CC animals exhibited significantly lower fat yields compared with heterozygous CT and homozygous TT animals (Table 1), which could be attributed to decreased milk yield with similar fat content.
As fat and protein content are crucial for cheesemaking, the potential use of the studied SNP in a selection scheme requires applied research to identify whether each additional unit of milk production compensates for the lower cheese yield due to lower protein percentage, based on farm prices and cheese yield for a certain type of cheese. Using current prices in Cyprus for sheep milk and the values of each genotype in Table 1, the producer is going to get an extra €1.012 per additional liter, whereas the price of the milk is reduced only by €0.012 (from €1.024 to €1.012) per liter due to the lower protein and fat content. Therefore, selecting for the T allele is expected to be beneficial; however, before incorporation of the SNP in any breeding program, negative effects of selection on other important traits (e.g., fertility) need to be excluded. To facilitate further research, the SNP has been incorporated into the sheep HD SNP chip developed by the International Sheep Genomics Consortium (ISGC) for functional studies (James Kijas, Csiro, Brisbane, Australia, personal communication).
Allele frequencies of the g.2982T>C polymorphism genotyped were similar to those previously observed (T:0.56; C:0.44;
). However, genotypic frequencies were found to deviate from the Hardy-Weinberg equilibrium in the studied population, in contrast to the previously observed frequencies from a single experimental flock (
), due to higher frequencies of TT and CT ewes carrying the favorable for milk yield T allele. This could possibly be explained by the fact that farmers select animals mainly based on milk yield, due to low meat prices during the last decade, whereas the selection indices used at the experimental flock (
Selection Index and Expected Genetic Progress in Chios Sheep. Agricultural Research Institute, Ministry of Agriculture and Natural Resources,
Nicosia, Cyprus1991
The candidate gene approach followed in the current study is an alternative to QTL and genome-wide association studies that can be very powerful in the identification of loci even with small effects on the trait if the candidate gene represents a true causative gene (
). Therefore, the polymorphism in the 3′ UTR mRNA of the ACAA2 gene that was found to be significantly associated with dairy traits can be either the functional polymorphism to which the genetic variance can be attributed or can be in linkage disequilibrium with the functional polymorphism or both. As the ACAA2 g.2982T>C SNP is in the untranslated region of the gene, and thus does not affect the protein produced, it may be a functional polymorphism playing a regulatory role by causing alterations in the expression levels of the gene. Identification of regulatory variants requires studying 2 alleles of a gene under identical circumstances and comparing the expression associated with each allele (
). For that reason, the expression of the ACAA2 variants was investigated in first-parity and age-, lactation stage-, and litter size-matched ewes that were managed, fed, and sampled uniformly to eliminate confounding sampling and environmental effects. The practical difficulties of obtaining those uniform conditions resulted in the selection of 3 animals per genotype. Three biological replicates are accepted as a minimum number of animals for mRNA expression analysis, as shown by recent published studies (
Liver X receptor a promotes the synthesis of monounsaturated fatty acids in goat mammary epithelial cells via the control of stearoyl-coenzyme A desaturase 1 in an SREBP-1-dependent manner.
). The expression of the gene in both the udder and the liver was significantly influenced by the animals' genotype, with ewes homozygous for the T allele showing significantly increased expression, by several orders of magnitude, compared with ewes homozygous for the C allele. The reasons why heterozygous mRNA expression in the udder exhibited a very high standard deviation needs to be elucidated by further research.
To assess whether the observed genotype effects in the expression of the ACAA2 gene could be attributed to AEI, an AEI assay using TaqMan probes was designed. The advantage of measuring AEI compared with total transcript levels is the reduction of the confounding effect of trans-acting factors, because the alleles are compared within and not across individuals (
). The results obtained are in support of differential expression of the 2 alleles of the ACAA2 gene due to a cis-acting mechanism. Differential allelic expression has been found in a large proportion of human (
) genes, and growing evidence suggests that polymorphisms in regulatory DNA and the resulting variability in gene expression can explain a significant proportion of disease susceptibility and quantitative trait phenotypic variance (
), the differential allelic expression of the ACAA2 gene varied among individuals in postmortem random samples, with an overall observed increase of T relative to C allele expression in the udder of heterozygous ewes close to 20%. Previous studies in cattle indicate that even small gene expression imbalances may result in large phenotypic variance of complex traits (
). Interestingly, mean AEI of the ACAA2 gene was only significant in the udder but not in the liver. Such tissue-specific AEI has also been reported for cattle genes (
). This finding, together with the observed difference between mammary gland and liver in the average mRNA gene expression in TT relative to CC ewes (11.8- and 2.8-fold increase, respectively), suggest an organ-specific differential expression. Preferential or higher expression in organs related to a quantitative trait is an important criterion for the selection of functional candidate genes (
As no other polymorphisms have been detected in the mRNA sequence of the ACAA2 gene, the differential mRNA expression of homozygotes and the observed AEI in the mammary gland could be attributed to the g.2982T>C SNP acting in cis to modulate gene transcription or mRNA survival (
). Analysis of the region harboring the SNP, using the MicroInspector program (Institute of Molecular Biology and Biotechnology, Foundation for Research and Technology-Hellas, Heraklion, Crete, Greece) found differences in potential miRNA binding sites between the 2 alleles (
). However, as this in silico analysis requires functional experimental support, the possibility of another cis-acting regulatory polymorphism in linkage disequilibrium with the g.2982C > T SNP (e.g., in the promoter region of the gene) cannot be excluded. Moreover, it is also likely that other trans-acting factors or mutations in other genes in linkage disequilibrium with ACAA2, found here to be associated with milk yield and protein content, may explain part of the QTL for milk, fat, and protein yields observed in the region harboring the ACAA2 gene (
The ACAA2 gene is involved in mitochondrial fatty acid elongation and degradation (KEGG database: http://www.kegg.jp/), by catalyzing the last step of the respective β-oxidation pathway. In agreement with our data, although ACAA2 is predominantly expressed in liver, it has also been reported to be expressed in the mammary gland (http://www.genecards.org/cgi-bin/carddisp.pl?gene=ACAA2;
). In mammals, excess energy is stored primarily as triglycerides, which are mobilized when energy demands arise. During periods of underfeeding or in early lactation, ruminants cover their increased energy demands by mobilizing fat from the adipose tissue (
). Fatty acids coming from triglycerides are taken up by the liver, where they are either used as energy source or converted to ketone bodies that may be released into the blood and used as energy or substrates for de novo fatty acid synthesis in the mammary gland. In addition, it is estimated that gluconeogenesis is usually increased 2- to 3-fold during early lactation to meet the demands of the mammary gland for lactose and triglyceride synthesis (
). The increased energy needs for gluconeogenesis and triglyceride synthesis are met primarily by increased fatty acid oxidation. Additionally, a recent transcriptomics analysis of the ovine mammary gland has shown increased expression of genes of the β-oxidation pathway during late pregnancy (
). This contributes to the development of the udder that is crucial for subsequent lactation, as the number of active secretory cells primarily determines milk yield (
). Therefore, although the mechanism by which the ACAA2 gene could be linked to increased milk yield and decreased protein percentage needs to be elucidated, increased ACAA2 expression at the mRNA level is likely to result in higher levels of the enzyme and, thus, elevated amount of energy and carbon substrates for mammary development and lactation.
CONCLUSIONS
In the current study, we demonstrated that the variants in the 3′UTR of the ACAA2 gene, which are associated with milk yield, protein percentage, and fat yield, are differentially expressed in homozygous ewes of each allele and exhibit AEI within heterozygotes in a tissue-specific manner, suggesting the existence of a cis-acting regulatory DNA mechanism. These findings support the hypothesis that the ACAA2 gene is a functional candidate affecting dairy traits.
ACKNOWLEDGMENTS
This work was funded by the Cyprus Research Promotion Foundation (Nicosia, Cyprus), the European Structural Funds, and the Cyprus University of Technology. G. Banos and A. Psifidi were supported by the Biotechnology and Biological Sciences Research Council (BBSRC, Swindon, Wiltshire, UK) Institute Strategic Programme Grant (BB/J004235/1; UK) and G. Banos also by the Rural & Environment Science & Analytical Services Division of the Scottish Government (Edinburgh, Scotland, UK). We thank the professional veterinarian John Pierroua (JnP Vet clinic, Avdimou, Cyprus) for obtaining biopsies. We also thank the head of the sheep and goat farming sector of the Ministry of Agriculture, Paniko Christoforou, the head (Giorgo Aspromalli) and technical staff of the governmental farm (Orites, Department of Agriculture, Paphos, Cyprus, for both Panico Christoforou and Giorgo Aspromali) and the owners of the commercial farms for giving us access to their farms and help in sampling and collection of data.
Identification of quantitative trait loci underlying milk traits in spanish dairy sheep using linkage plus combined linkage disequilibrium and linkage analysis approaches.
Section 2 ICAR rules, standards and guidelines for dairy production recording. Approved by the General Assembly held in Berlin, Germany, May 2014. ICAR (International Committee for Animal Recording),
Rome, Italy2014
I: Segregation, Cloning, SNP Detection and Its Association with Economic Traits of 5 Bovine Genes. II: Quantitative Trait Loci Analysis of Meat Quality Traits in Pigs. PhD thesis. Chinese Academy of Agricultural Sciences,
Beijing, China2008
Selection Index and Expected Genetic Progress in Chios Sheep. Agricultural Research Institute, Ministry of Agriculture and Natural Resources,
Nicosia, Cyprus1991
Liver X receptor a promotes the synthesis of monounsaturated fatty acids in goat mammary epithelial cells via the control of stearoyl-coenzyme A desaturase 1 in an SREBP-1-dependent manner.
30588-X&id=gr1.jpg){kind=link}
30588-X&id=gr2.jpg){kind=link}
30588-X&id=gr3.jpg){kind=link}
30588-X&id=gr4.jpg){kind=link}