Figures
Abstract
Adaptation of global food systems to climate change is essential to feed the world. Tropical cattle production, a mainstay of profitability for farmers in the developing world, is dominated by heat, lack of water, poor quality feedstuffs, parasites, and tropical diseases. In these systems European cattle suffer significant stock loss, and the cross breeding of taurine x indicine cattle is unpredictable due to the dilution of adaptation to heat and tropical diseases. We explored the genetic architecture of ten traits of tropical cattle production using genome wide association studies of 4,662 animals varying from 0% to 100% indicine. We show that nine of the ten have genetic architectures that include genes of major effect, and in one case, a single location that accounted for more than 71% of the genetic variation. One genetic region in particular had effects on parasite resistance, yearling weight, body condition score, coat colour and penile sheath score. This region, extending 20 Mb on BTA5, appeared to be under genetic selection possibly through maintenance of haplotypes by breeders. We found that the amount of genetic variation and the genetic correlations between traits did not depend upon the degree of indicine content in the animals. Climate change is expected to expand some conditions of the tropics to more temperate environments, which may impact negatively on global livestock health and production. Our results point to several important genes that have large effects on adaptation that could be introduced into more temperate cattle without detrimental effects on productivity.
Citation: Porto-Neto LR, Reverter A, Prayaga KC, Chan EKF, Johnston DJ, Hawken RJ, et al. (2014) The Genetic Architecture of Climatic Adaptation of Tropical Cattle. PLoS ONE 9(11): e113284. https://doi.org/10.1371/journal.pone.0113284
Editor: Olav Rueppell, University of North Carolina, Greensboro, United States of America
Received: April 15, 2014; Accepted: October 21, 2014; Published: November 24, 2014
This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.
Data Availability: The authors confirm that all data underlying the findings are fully available without restriction. We have submitted data to the Animal QTL database where they will be publicly available. We have received the repository identifier www.animalgenome.org/repository/pub/USDA2013.1112/associated with this data deposition. The gene expression data used in this paper were previously published and are available from the GEO database under accession numbers GSE44030 and GSE2554. A copy of the data will be available at CSIRO Data Access Portal https://wiki.csiro.au/display/dmsdoc/Home.
Funding: The data collection was funded by the Cooperative Research Centre for Beef Genetic Technologies, the Cooperative Research Centre for Cattle and Beef Quality and their core partners. The Cooperative Research Centres have closed and funding for the analysis and write-up has come from non-specific funding from CSIRO through its Food Futures Flagship. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: William Barendse is a current academic editor of PLOS ONE. This does not alter our adherence to PLOS ONE Editorial policies and criteria. Rachel Hawken is an employee of Cobb-Vantress and Kishore Prayaga is an employee of Zoetis Inc. We have read the journal's policy, and declare the following: This does not alter our adherence to PLOS ONE policies on sharing data and materials. Furthermore, there is no company confidential information or intellectual property such as patents or trade secrets or commercial know-how retained by these authors because the work was performed while they were employed by CSIRO, before they left for their current places of employment.
Introduction
Cattle are a key component of tropical rural sustainability, wealth, and food production for small holder farms, which are vital for global food security [1], [2]. Tropically adapted cattle breeds, usually of indicine ancestry, are less productive than European taurine cattle under favourable conditions [3], [4]. However, they are easier to maintain under conditions of low inputs, and remain productive under the onslaught of parasites, tropical diseases and high heat loads [5], [6], [7]. Nevertheless, lower innate productivity is a serious limitation to sustainable intensification of cattle productivity in the tropics.
The barriers to uptake of more productive cattle in tropical environments are lack of quality nutrition during the dry season, high temperature and humidity during the wet season, high water requirements, and devastating tropical parasites and infections, amongst others [8], [9], [10], [11]. Some breeders, particularly in developed countries, have responded by generating taurine-indicine composite cattle and then selecting the best animals to form new, composite breeds for particular environments or production systems.
Indicine and European taurine cattle differ substantially in many respects, not just in their productive potential in tropical environments, but also in traits such as temperament and appearance. First generation taurine x indicine crossbred cattle show heterosis including resistance to parasites and increased productivity, allowing them to outperform either taurine or indicine purebred cattle in tropical production systems [12], [13]. However, the causes of the decline of adaptation and productivity in subsequent generations may depend on the nature of the heterosis. This may be caused in part either 1) by genes of large, dominant effect that are near to being fixed for opposite alleles in the two subspecies, or 2) due to epistatic or pleiotropic interactions between genes within one subspecies that are then upset by new genetic combinations generated in the F2 or subsequent generations, or 3) due to recessive deleterious mutations in each group that are diluted in the first generation but through segregation emerge in the second generation, or 4) due to overdominance at a large number of loci, in which individuals in the F1 generation would on average have more loci with overdominant effects than individuals in the F2 generation [14], [15].
There has been anecdotal observation that some F2 crossbred cattle either perform substantially less favourably than F1 cattle, or seem to inherit the weak points of both parental breeds, which might suggest genes of large effect. However, evidence for genes of large effect in genome wide analyses has so far been lacking. Previous studies of crossbred animals have focussed on meat quality and efficiency of production in feedlots, and so far, only meat tenderness and reproductive performance have shown any evidence of genes of moderate effect [16], [17].
We have used a large sample of crossbred cattle of two types in genome wide association studies (GWAS) to analyse 10 traits that are important in extensive tropical environments. Such traits include overall productivity and yearling weight, heat and parasite tolerance, and temperament. So far, GWAS of parasite resistance traits have been performed with small sample sizes and have not identified any genes of large effect [18]. Here we show that contrary to the general expectation of a large number of genes each of infinitesimally small effect, that 9 of the 10 traits were affected by genes of moderate to large effect size. This suggests that for tropical productivity, the genetic differences between indicine and European taurine cattle includes the segregation of genes of large effect.
Results and Discussion
To determine whether there are genes of large effect segregating in taurine x indicine composite cattle we performed a genetic analysis of ten traits that are relevant for improved performance in tropical production systems (Table 1, Table S1). Here we consider the SNP to have a large size of effect if it accounted for more than 5% of the genetic variance, which translates to approximately 0.3 S.D. in allele substitution effect, due to the overall rarity of such genetic effects when measured in large sample sizes [16], [19], [20], [21]. We performed a genome wide association study (GWAS) using 4,662 pedigreed cattle and high density bovine single nucleotide polymorphism (SNP) arrays. There were 2,112 animals of the Brahman breed, which is largely of indicine origin [22] and 2,550 cattle of a set of composite cattle with varying amounts of European and African taurine and indicine ancestry, all raised in tropical environments (Figure 1A). We also made use of genotypic datasets of 81 Angus (taurine) and 91 Nelore (indicine) cattle from previous experiments including the Bovine HapMap [23], [24], [25] to quantify the population structure and breed composition of our samples (Figure 1B, C). The Brahman sample showed small differences to the Nelore sample in a multi-dimensional scaling analysis of the genotypes. In comparison, the Tropical Composites showed substantial variation not only between the taurine and indicine extreme captured by the first component, but also substantial variation captured by the second component, which may be due to African ancestry of the animals descended from the Belmont Red [26], Africander, and Senepol breeds.
(b) Distribution of indicine% across the two samples showing a mode at ∼30% for Tropical Composite (red) and ∼95% for Brahman (blue) cattle. (c) Multi-dimensional scaling plot showing a clustering of the Brahman (blue) and Tropical Composite (red) cattle relative to the reference samples of Angus (black) and Nelore (green). (d) Genetic relationship matrices based on genotype similarity among individuals for Tropical Composite (left) and Brahman (right) samples where the diagonal blocks show paternal half-sib families.
Trait measurements included 1) measures of coat colour, score (length), and rectal temperatures (Table S1), which are relevant to heat tolerance, 2) temperament and sheath/navel score, which are relevant to the ease of care and management needed in tropical agriculture, 3) response to parasites (ticks, worms, flies), which relate to tropical diseases and parasites, and 4) yearling weight and body condition score, which measure the overall response of the animal to the conditions. The estimated heritabilities and genetic correlations were found to be similar in both samples (Table 2). Seven out of the 10 traits (Table 1, Table 2 and Table S1) showed moderate to high heritability (i.e., between 0.30 and 0.60), and the genetic correlation between traits in males and females was almost always close to one (Table S2). Two of the 10 traits were lowly heritable in both samples, the size of lesions due to fly infestation and rectal temperature. The first depends on Stephanofilaria infection from the fly salivary glands [27] and the second is a trait under negative feedback control. Only tick burden showed a moderate heritability in the Tropical Composite sample and a low heritability in the Brahman sample. This trait had the smallest number of observations, so the difference could be due to sampling. Previous estimates of heritabilities and genetic correlations for these traits have been estimated on smaller subsets of these data [28]. Overall, when comparing the 45 estimates of genetic correlations across the two samples, we did not find a significant difference (2 tailed t-test, P = 0.874) between them. This suggests that although the animals differed substantially in the amount of taurine or indicine ancestry, the amount of genetic variation and the additive genetic relationships between traits was similar in the two samples.
To determine the role of indicine origin in the genetic architecture of these traits we quantified the extent of admixture percentage of each animal. Using genotypes for 228,268 SNP for the two reference samples of Angus and Nelore, we estimated the genomic proportion that can be attributed to indicine (indicine%) ancestry in each animal. When Hereford, Shorthorn and Gir were used as the reference samples they provided highly correlated estimates (r>0.91), demonstrating that the specific taurine or indicine breed was not important for quantifying indicine% in these samples (Figure S1). The mean indicine% of the Brahman sample was 96.7% while that of the Tropical Composite sample was 27.2% (Figure 1B). The estimated effect of indicine% was statistically significant (P<0.005) for 7 out of the 10 observed traits in the Tropical Composite sample, and for 4 traits in the Brahman sample (Table 1 and Table 3). In both samples, indicine% lightened coat colour, short coat score and caused more pendulous penile sheaths and increased weight and body condition score. In addition, indicine% reduced worm infestation and worsened temperament in the Tropical Composite sample. These changes were consistent with the general scientific and anecdotal description of indicine and European taurine cattle. Nevertheless, although these effects were statistically significant, the effect sizes were small. The observation of the empirical density having a clear mode at ∼30% of indicine% (Figure 1b) encouraged the exploration of the effect of indicine content separate for each group. The switching of size and/or sign in the estimated effect of indicine% in the two groups indicates a possible non-linearity of indicine gene effect in the transition from mostly taurine to mostly indicine. This was clearly the case for EPG.
In the GWAS, nine of the ten traits showed SNP of genome-wide significance (P<5×10E-8) (Table 4, Figure 2, Table S3, Table S4 and Figure S2 and Figure S3). Genetic regions identified by these SNP explained between 3.6% and 71.3% of the genetic variance (Table 4), with the proviso that these values were not shrunk through whole genome simultaneous estimation. Evidence of genetic control including genetic variants of large effect, as opposed to only due to genes of small effect, was observed for coat colour, coat score, penile sheath, yearling weight, and body condition (Figure 2A). While coat colour in cattle is known to be controlled by a small number of genes of large effect [29], the finding for other traits including penile sheath and yearling weight has not been demonstrated before. Some of the alleles affecting parasite resistance exceeded 5% of the genetic variance, larger than previous estimates based on smaller sample sizes [18].
(b) Manhattan plots of the significance (-logP on the y-axis) of each SNP in genome order (x-axis) for Sheath score (left panels) and Coat Colour (right panels) for Tropical Composite (upper panels) and Brahman (lower panel) cattle. Note the similarity on BTA5 for Sheath score in both samples but different genes for colour, (c) Average frequency of the forward allele for sliding windows of 100 consecutive SNP along BTA5 at 1 SNP pace for Brahman (blue) and Tropical Composite (red) cattle relative to the reference samples of Angus (black) and Nelore (green). The insert shows regions of divergent allele frequencies between taurine and indicine cattle on BTA5 from 20 to 60 Mb. (d) Heat map of LD (r2) on BTA5 in the Tropical Composite sample, with red dots corresponding to r2>0.1, showing LD blocks spanning several Mb. (e) The effect of BTA5 on Sheath score in Tropical Composite and Brahman cattle was not due to indicine% but to a major gene. 50 SNP were selected from 5 regions of BTA5 with divergent alleles in Angus and Nelore cattle. The average additive association (based on –log(P), y-axis) of these 50 SNP across the ten traits was calculated (green bars) and compared with the association of 50 SNP (orange bars) also with divergent alleles in Angus and Nelore cattle but randomly located to regions other than BTA5 or BTAX. A total of 10 random selections were chosen and the average plotted. (f) Genes close (< = 3 Kb) to SNP significantly associated (P<0.0001) with Sheath in both breeds.
Although the size of the genetic variance in each sample was similar for all traits, as measured by the estimated heritability, the regions of the genome with the largest effects were generally different between the two samples, except for sheath score. For example, SNP related to the genes ASIP and PLAG1, which had important effects on coat colour or body condition score in both samples, accounted for different amounts of genetic variance in each sample (Table 4). Animals of the two groups show significant differences in coat colour and structure, response to parasites, size and body shape (Figure 1A and Table 1), and so major genetic differences between them may be expected. Traits with low heritability showed few associations that were genome wide significant and little to no overlap between the two samples, suggesting that even larger sample sizes will be needed to analyse these complex traits in detail.
The most notable results were for penile sheath and coat colour in the Tropical Composite sample, in which almost the entire length of bovine chromosome (BTA) 5 was associated to both traits (Figure 2B). Even restricting the analysis to a SNP association value of –logP>20 resulted in more than 50 Mb of the chromosome associated to either trait, with a massive overlap of these BTA5 segments between traits. We took several approaches to determine if this was an artefact of analytical inputs. We found evidence of significant differences in allele frequency between Angus and Nelore in several regions of BTA5. However, other regions of the genome with such differences in allele frequency did not show effects on coat colour or penile sheath, so the allele frequency differences do not explain the observations (Figure 2C and E). Although there was evidence of long range linkage disequilibrium (LD) across most of BTA5 (Figure 2D) in this sample, and such LD might explain correlated –logP values over extended regions of a chromosome, BTA5 is not special compared to other chromosomes in this data set and indeed shows lower levels of LD between SNPs at different distances than the average of other autosomes (Dataset S1 taken from Porto-Neto and coworkers [30]). Regions of the genome containing PLAG1 [31] and the Slick locus [32] gene also show broad regions of significance in this sample, with PLAG1 on BTA14 showing –logP>8 over 4.4 Mb for yearling weight. The region containing the Slick gene showed –logP>20 over 28 Mb. In the Tropical Composite breeds studied here, the gene is likely to be derived from crossing with bulls of part Senepol ancestry, a known source of this trait. Since the colour and sheath score genes are unlikely to have a similar defined origin, this suggests that more than LD is involved. To determine whether these effects were due to recent crossbreeding, we excluded all animals of recent crossbred origin, restricting the analysis to 921 animals that have been closed to introgression for many generations and found the same broad extent of association. We examined the haplotype structure of BTA5 in this restricted sub-population of animals from both breed types and found that it contained significant haplotypes associated to both coat colour and penile sheath score extending over at least 10 Mb in size (Dataset S2). Moreover, in addition to these two traits, this region is also associated to yearling weight, body condition score, and parasite resistance (Table 5, Table S5) and female reproductive success in previous studies [17]. Given this combination of traits in one genomic region, we speculate that breeders have maintained haplotypes of genes for coat colour and penile sheath conformation because they are linked to parasite resistance and productivity in tropical environments, possibly because these traits are signatures of Bos indicus ancestry. However, we found that genetic variants leading to more pendulous penile sheaths, derived from Brahman breed ancestry, were associated with reduced resistance to parasites, lower yearling weight and better body condition scores in these samples. One would have expected improved resistance to parasites from the Brahman alleles.
To determine whether any of the linked effects may instead be pleiotropic in origin, in particular effects that might be antagonistic, we evaluated the effect of each SNP on all traits using the summation of standardized effects (Table 5 and Table S5). Genes on BTA5, 14 and 20 in particular, showed antagonistic effects for body condition and yearling weight, that is, the favourable allele for body condition was the unfavourable allele for yearling weight rate, suggesting an underlying genetic antagonism consistent with the slightly negative genetic correlation between these traits in both samples. Furthermore, some of these same favourable alleles for body condition or yearling weight have unfavourable effects on several of the parasite resistance traits, again consistent with the slightly negative genetic correlations between these traits. This suggests that selection for these genes of large effect on weight or body condition will come with a trade-off of reduced parasite resistance, which will need to be accommodated by selection for parasite resistance using the rest of the genome. Two of the genes that have antagonistic effects on more than one trait have been identified before, due to their large effects on the visible phenotype, viz. PLAG1 on BTA14 associated with hip height and female reproductive performance [31], [33], and the Slick gene on BTA20 associated with coat score, ability to sweat, rectal temperatures [32], and hence ability to be active and feed during the heat of the day in the tropics. In this study we also identified an extended genetic region centred around the MSRB3 gene on BTA5 as having large effects on several traits.
To identify the genes that underpin effects on yearling weight and body condition, we combined information on location of significant SNP effects and gene expression profiles of Tropical Composite cattle under tropical dry season conditions, where typical reductions in weight and body condition occur. SNPs that were significantly associated with yearling weight or body condition score in both samples identified a list of genes that had been examined for gene-expression differences using muscle biopsies before and after a period of low food availability for the Belmont Red breed [34], [35]. Of the 175 genes associated with significant SNP for condition score, 63 showed significant differences in gene expression before and after the period of low food availability, while 52 of 213 genes for yearling weight showed differences in gene expression (Table 6). Of great relevance for periods of undernutrition, some of these genes were linked to energy metabolism (LDHA, IGFBP4, ATP7A, STEAP2 and MON2), muscle metabolism (MYH2, ASPH and APOBEC2), fertility (INHBC) and the immune system (LYN, CHD7, ADORA2B, IRAK3, CORO2A and TOX). A large number of these genes were associated with the region containing MSRB3 on BTA5 and PLAG1 on BTA14 as noted above. This suggests not only that some of the genetic variation affects the phenotype through differences in gene expression, that clusters of genes affecting these traits are co-ordinately expressed, but also that genetic variation in physiologically relevant gene networks are deployed by the animal to adapt to tropical conditions.
Climate change is expected to expand conditions of the tropics to more temperate environments, which may impact negatively on global cattle production unless steps are taken to improve cattle adaptation and productivity. Not only will temperate regions become warmer but climatic variability will increase and the range of devastating tropical and subtropical parasites will extend into temperate regions [36]. Due to the limits on arable land [37] and increased requirement by other industries for grain and other feedstuffs currently fed to livestock, cattle production systems must become more efficient [1]. Cattle will be expected to improve productivity and increase adaptation to heat, diseases, and poor quality feed on marginal land, to add to the total food available. Our results point to several important genomic regions that have large effects on adaptation without compromising productivity. They also show that the genetic architecture for adaptation to the tropics includes alternative alleles of large effect rather than co-adaptation of the genome that breaks down during crossbreeding between the subspecies. This suggests that a focus on these genes would rapidly improve the predictability of crossbreeding of taurine and indicine cattle. Finally, sustainably intensifying livestock systems in tropical regions may need genetic safeguards to ensure that productivity is raised while also adapting to climate change. The results of this study are a necessary step towards achieving that goal.
Materials and Methods
Animals
Animal Care and Use Committee approval was not obtained for this study because no new animals were handled in this experiment. The experiment was performed on trait records and DNA samples that had been collected previously. The resource population used in this study had been established by the Cooperative Research Centre for Beef Genetic Technologies (Beef CRC) to understand the genetic links between adaptability and components of herd profitability in northern Australia [38]. DNA had been extracted from venous blood samples from each animal using Qiagen kits as previous described [39].
Tropical adaptation-related traits
Tropical cattle are exposed to several environmental stressors such as hot and humid conditions, restricted water supply, periodic challenges from exo- and endo-parasites and the diseases they transmit, and seasonal under-nutrition due to low protein availability and feed digestibility. In this study, we examined 10 traits (Table S1): flight time (FT), rectal temperature (TEMP), endoparasite eggs per gram measured in faeces (EPG), penile sheath score (SHEATH) expressed as the correlated trait navel score in females, coat colour (COLOUR), buffalo fly lesions (FLY), tick infestation (TICK), coat score (COAT), condition score (COND) and yearling weight (YWT) [28]. While some of the traits are obviously related to tropical conditions, such as direct measures of heat or parasite resistance, all traits have relevance to tropical agriculture. Flight time is a measure of temperament, how excitable or nervous the animal is when exposed to humans, which is relevant because in much Australian tropical agriculture the animals are not regularly exposed to humans and hence can react and cause injury to themselves or their handlers. Coat colour and length are directly related to heat absorption and radiation and affect heat tolerance and water usage. Sheath or navel score are related not only to the performance of the bull, as animals with pendulous sheaths are more likely to be injured, but it is thought by many cattle breeders that pendulous sheaths, navels, dewlaps and long ears help to radiate heat, although this has not been proven. To improve consistency between subjective scores and across different sites, these scores were taken by only four trained individuals. Repeated measures on these tropical adaptive traits were collected at various post-weaning ages with specific measures defined based on the biological significance of the age at which the measurement was taken and to maximise the number of records available for analyses across all cohorts.
Trait heritability and correlations
The estimates of genetic and phenotypic variances, heritability of each trait and phenotypic and genetic correlations between traits were derived from bivariate animal model analyses of Brahman and Tropical Composite performed independently using VCE 6.0.2 (ftp://ftp.tzv.fal.de/pub/vce6/). Genetic parameter estimates were calculated based on the genotyped animals and their pedigree consisting of 5 generations of ancestors. Analytical models included the fixed effects of contemporary group (combination of sex, year and location), age of dam and estimated indicine percent of the individual as covariate; as well as animal as random additive effect and the random residual component. Age of dam was found to be statistically significant as a covariate. Animals that are older, in a tropical area, will have had a longer time to adjust to heat, parasites, diseases, poor fodder, than younger animals, who would also produce less milk than older mothers. These differences will include the development of antibodies, some of which will be passed on to their offspring. Therefore, dams that are older will likely be more successful in raising offspring, and indeed, we see strong significant effects for this covariate. For YWT the additional covariate of AGE12 was also fitted into the model. AGE12 is the average age of the animal at days in which body weight was recorded (from all the weight measurements taken between 300 and 420 days of age) and so is a proxy for yearling age. The genetic variance and correlation for each trait was also calculated between males and females. Large differences in these statistics might be due to differential responses of the animals to the environment or because the traits may have different genetic architectures in males and females. To compare the correlation coefficients found for these traits between Brahman and Tropical Composite animals, the differences between correlation coefficients for each pair of traits between Brahman and Tropical Composite were compared using a paired t-test. This generated a single test for the 45 comparisons, rather than comparison the differences of the correlation coefficients to the standard errors and then applying a False Discovery Rate model for the number of tests performed.
Genotypes
For the present study, we used 2,112 Brahman and 2,550 Tropical Composite cattle from the resource population genotyped using either the BovineSNP50 [40] or the BovineHD BeadChip (Illumina Inc., San Diego, CA) that includes more than 770,000 SNP.
All SNP had been mapped to the UMD build assembled by the University of Maryland [41] using the updated version 3.1 of the genome (available from Genbank accession DAAA00000000.2 and at http://www.cbcb.umd.edu/research/bos_taurus_assembly.shtml). Animals genotyped using the lower density array had their genotypes imputed to higher-density based on the genotypes of relatives, consisting of 589 Tropical Composite and 304 Brahman animals including all available sires, that had been genotyped using the BovineHD BeadChip. The imputation was performed within breed using as reference 519 Brahman and 351 Tropical Composites genotypes using the BovineHD (Illumina) and 30 iterations of BEAGLE [42], which resulted in a final dataset of 729,068 SNP genotypes per individual as reported in Bolormaa et al. [24]. For the estimation of indicine content, the full genotype dataset was filtered by linkage disequilibrium (LD) to reduce redundant information and optimise computation utilization. The LD filter was applied using PLINK v1.07 [43] in a sliding window consisting of 50 adjacent SNP, and if r2>0.5 was detected between a pair of SNP one of the SNP was removed, and then LD for the window was re-calculated. Once no more pairs of SNP had r2>0.5 the window moved 10 SNP along the chromosome. This procedure yielded a dataset containing 227,085 SNP. For SNP association analyses, the SNP were consistently encoded for all traits as AA, AB and BB using the TOP/BOT encoding of Illumina (http://http://res.illumina.com/documents/products/technotes/technote_topbot.pdf checked 30 July 2014) and then converted into numerical values of 0, 1, and 2 B alleles.
Identification of population sub-structure and indicine content
We used PLINK v1.07 to calculate multi-dimensional scaling (MDS) and genetic relationship matrices based on the genotypes to quantify the sample substructure for both the full dataset as well as the LD filtered dataset. The Angus and Nelore data were used as representatives of pure European taurine and pure indicine animals respectively. We note that indicine percent lines up with the first principal component of a principal component analysis, whereas African ancestry lines up with the second principal component, where the full diversity of cattle is included in the same analysis [44].
The Tropical Composite sample consisted of beef industry lineages formed using indicine, Sanga and taurine cattle [38], and included animals with African taurine and no reported indicine ancestry, such as some sectors of the Belmont Red breed. The Brahman breed in Australia started in the 19th century from various Indian cattle including animals from the Melbourne Zoo, upgraded by American Brahman and the Indu-Brazilian breeds so also includes a small proportion of taurine ancestors [22], [45], [46]. It therefore includes breeds such as the Kankrej (Guzerat), the Ongole (Nelore), and the Gir (Gyr). To estimate the amount of taurine and indicine content in the Tropical Composite and Brahman animals, HD genotypes for 81 Angus (Beef CRC) and 91 Nelore cattle were used as a reference. Estimates of indicine ancestry were also obtained using 55 Hereford, 54 Shorthorn, 44 Angus and 50 Gir. Some Angus, and the Hereford and Shorthorn samples were obtained from the Beef CRC database and some Angus, and the Nelore and Gir samples were obtained from the Bovine Hapmap [23], with some Nelore and Gir animals sampled from Brazil [25]. Ancestry was estimated using Admixture software [47] on the LD filtered dataset (∼228 K SNP) under supervised mode. Most of these Tropical Composite cattle were descended from the Hereford or Shorthorn breeds as taurine ancestors. Thereafter to evaluate the potential impact of breed on the indicine estimates, these were also obtained using different combinations of Hereford, Shorthorn and Gir. The correlation between different estimates of indicine content were>0.91 for these comparisons.
To estimate the effect of indicine percent on adaptability-related traits a linear model was fitted using SAS (SAS Inst., Cary, NC). The covariates used for the estimation of genetic parameters were used as fixed effects. Given the number of comparisons, the type I significance threshold was Bonferroni adjusted to α = 0.005.
Genome-wide association studies
Genome-wide association studies (GWAS) were performed separately within each breed and for each of the ten traits using the final dataset with 729,068 SNP. The GWAS were performed one SNP at a time using the same univariate linear mixed models stated above, which included the fixed effects of contemporary group (combination of year and location), age of dam and estimated indicine percent of the individual as covariate; as well as animal as random additive effect and the random residual component and the SNP genotype (recoded as 0, 1 or 2) as an additional linear covariate. Solutions for the SNP effects and associated P-values were obtained using Qxpak5 [48]. As the coefficient for each SNP is provided as a signed number, the signs of significant coefficients of the same SNP for different traits can be compared to determine whether the SNP shows effects in the same or opposite directions between two traits. To determine whether fitting indicine percent substantially changed the outputted results, we compared GWAS output with and without indicine percent in the Brahman sample and found that the allele effects across traits had an average correlation of 0.97 (Table S6).
False discovery rate
Following Bolormaa et al. [24], the false discovery rate was calculated as
where P is the P-value tested (e.g., 0.0001), S is the number of SNP significant at the P-value tested and T is the total number of SNP tested (i.e. 729,068).
Percentage of genetic variance explained by each SNP
The percentage of the genetic variance explained by the i-th SNP was estimated for each sample separately according to the following formula:
where pi and qi are the allele frequencies for the i-th SNP calculated for each breed, is the estimated additive effect of the i-th SNP on the trait under analysis,
is the REML estimate of the (poly-)genetic variance for the trait, and
is the estimated additive variance of the i-th SNP in the absence of dominance [49].
Test for pleiotropy
In both the Brahman and the Tropical Composite sample, the effects of 729,068 SNPs estimated from the GWAS were divided by their associated standard errors to obtain a t-value corresponding to the studentized SNP effects [50]. A multi-trait test of the i-th SNP was performed by storing its studentized effects across the 10 traits in the 10×1 vector ti. Then, the quadratic form ti'V-1ti, where V is the correlation matrix among the SNP effects, is distributed approximately as a chi-squared with 10 df under the null hypothesis that the SNP does not affect any of the traits. The correlation matrix V was calculated using the estimated SNP effects across the 729,068 SNPs. Custom code in FORTRAN 95 was developed to perform these operations. SNP were taken as having a significant pleiotropic effect when the quadratic form t' V-1 t exceeded 29.588 (i.e., P<0.001 from a Chi-squared distribution with 10 degrees of freedom). Here we are looking for pleiotropic effects for SNP that had already been demonstrated to be significantly associated to the traits of interest after correction for multiple testing, so no further correction for multiple testing was required.
Gene-expression of positional candidate genes in muscle biopsies under environment stress
Positional candidate genes were defined as those genes within 3,000 bp of a SNP significantly associated (P<0.01) to either body condition score or yearling weight, and in which the same allele was favourably associated with the trait in both samples. To further explore those associated genes and link their responses to a potential environmental stressor, we evaluated their expression using data from a previously described cattle undernutrition experiment [34], [35]. The artificial feed reduction was to mimic natural drought conditions when food is of low protein composition and poor digestibility and cattle lose weight. Muscle biopsies from 12 Tropical Composite animals were collected before and after an undernutrition period of 114 days. Gene-expression was assessed using microarrays (ViaLactia Bioscience in collaboration with Agilent) that included more than 21,000 probes representing around 19,500 protein-coding genes. Genes that had significant (P<0.05) expression change in the compared period were selected and matched to those candidate genes derived from the GWAS.
Estimated effect of the number of “taurine alleles” on adaptation-related traits
There are many SNP across the cattle genome that have one allele fixed or near fixation in indicine cattle while the alternative allele is fixed or near fixation in taurine animals. If a SNP had an allele with frequency>0.95 in Angus (taurine) and <0.05 in Nelore (indicine) animals, that allele was called the “taurine” allele and the alternative allele the “indicine” allele and labelled as a highly divergent SNP.
We noted a QTL of large effect on BTA5. To determine whether highly divergent SNP per se were associated to trait values or were contributing to the apparent size of the QTL effect, we collected genotypes for 50 highly divergent SNP distributed across 5 regions of BTA5. As a control, we collected 10 random selections each of genotypes for 50 highly divergent SNP for all sectors of the genome except BTA5 and BTAX. Then, the “indicine” and “taurine” alleles were identified and the number of “taurine” alleles was summed for each animal. To estimate the effect of the sum of “taurine” alleles on adaptability-related traits a linear model was fitted using SAS (SAS Inst., Cary, NC). The model contained the same fixed effects as the model for estimated indicine percent.
Supporting Information
Figure S1.
Correlation of breed composition (indicine%) estimated using different sets of animals as reference populations.
https://doi.org/10.1371/journal.pone.0113284.s001
(DOCX)
Figure S2.
Manhattan plots of genome wide association for FT, TEMP, EPG, SHEATH, and COLOUR.
https://doi.org/10.1371/journal.pone.0113284.s002
(DOCX)
Figure S3.
Manhattan plots of genome wide association for FLY, TICK, COAT, COND, and YWT.
https://doi.org/10.1371/journal.pone.0113284.s003
(DOCX)
Table S1.
Trait definitions and abbreviations.
https://doi.org/10.1371/journal.pone.0113284.s004
(DOCX)
Table S2.
Bi-variate analysis for same trait in different sexes: Additive genetic variance (Va), heritability (h2) and genetic correlation (rg).
https://doi.org/10.1371/journal.pone.0113284.s005
(DOCX)
Table S3.
Number of significant SNP from GWAS at different levels of stringency for each trait and false discovery rate between brackets.
https://doi.org/10.1371/journal.pone.0113284.s006
(DOCX)
Table S4.
Summary table of SNP association: overlap regions between Brahman and Tropical Composite (P<0.0001 both breeds).
https://doi.org/10.1371/journal.pone.0113284.s007
(DOCX)
Table S5.
Selected Pleiotropy results: Studentized SNP effects across the 10 phenotypes for candidate SNP and its closest gene. Highlighted in yellow are effects> 2.0 in absolute value.
https://doi.org/10.1371/journal.pone.0113284.s008
(DOCX)
Table S6.
Correlation between SNP additive effects (n = 729,068), and statistical model R2 estimated in the Brahman population with and without fitting the covariate “% indicine” in the model.
https://doi.org/10.1371/journal.pone.0113284.s009
(DOCX)
Dataset S1.
Linkage Disequilibrium decay on BTA5 compared to the other autosomes.
https://doi.org/10.1371/journal.pone.0113284.s010
(XLSX)
Dataset S2.
Haplotype distribution in Tropical Composite animals for BTA5 affecting NAVEL and COLOUR.
https://doi.org/10.1371/journal.pone.0113284.s011
(XLSX)
Acknowledgments
This work was performed within the Cooperative Research Centre for Cattle and Beef Quality, the Cooperative Research Centre for Beef Genetic Technologies and their core partners including Meat and Livestock Australia. The resource population was constructed with assistance from AgForce Queensland, the Australian Agricultural Company, C & R Briggs, Consolidated Pastoral Company, J & S. M Halberstater, S. Kidman & Co, GE & J McCamley, MDH Pty Ltd, North Australian Pastoral Company and Stanbroke Pastoral Company. We acknowledge the contributions of staff of the Queensland Department of Agriculture, Fisheries and Forestry (Dick Holroyd, Brian Burns, Tim Grant, Mick Sullivan, Tracy Longhurst, Neil Cooper, Peggy Olsson, Ken Olsson, Alena Ladd, Rebecca Hall, Trudy Obst, Kate Pratt, Steve O'Connor, John Hughes, Debra Corbet, Russ Tyler, Cindy McCartney, Lisa Hutchinson, Ruth Coates, Jared Heath, Megan Debney, Joanne Campbell, plus stock handlers), CSIRO (Warren Sim, Rob Young, Graeme Halford (dec), John Quilty, Lindsay Cranston, Scott Prewett, Peter Corbett, Ernie Young, Adele Weaver, Keith Harder, Wayne Flintham, Phil Orchard, Paul Williams, Nick Corbet), the Animal Genetics and Breeding Unit in Armidale (Steve Barwick, Andrew McCann, Jim Cook) and the University of Queensland (Dianne Vankan, Emily Piper). We thank Mario Herrero and James Kijas for discussion.
Author Contributions
Conceived and designed the experiments: AR WB HMB KCP DJJ. Performed the experiments: GF TSS HMB. Analyzed the data: LRPN AR JMH EKFC DJJ. Contributed reagents/materials/analysis tools: TSS JFG SB MEG RJH GF SAL WB. Wrote the paper: WB AR LRPN SAL.
References
- 1. Foley JA, Ramankutty N, Brauman KA, Cassidy ES, Gerber JS, et al. (2011) Solutions for a cultivated planet. Nature 478:337–342.
- 2. Herrero M, Thornton PK, Notenbaert AM, Wood S, Msangi S, et al. (2010) Smart investments in sustainable food production: Revisiting mixed crop-livestock systems. Science 327:822–825.
- 3. Frisch JE, Vercoe JE (1977) Food intake, eating rate, weight gains, metabolic-rate and efficiency of feed utilization in Bos taurus and Bos indicus crossbred cattle. Animal Production 25:343–358.
- 4. Frisch JE, Vercoe JE (1984) An analysis of growth of different cattle genotypes reared in different environments Journal of Agricultural Science. 103:137–153.
- 5. Davis GP (1993) Genetic parameters for tropical beef cattle in northern Australia: a review. Australian Journal of Agricultural Research 44:179–198.
- 6. Menjo D, Bebe B, Okeyo A, Ojango J (2009) Analysis of early survival of Holstein-Friesian heifers of diverse sire origins on commercial dairy farms in Kenya. Tropical Animal Health and Production 41:171–181.
- 7. Guis H, Caminade C, Calvete C, Morse AP, Tran A, et al. (2012) Modelling the effects of past and future climate on the risk of bluetongue emergence in Europe. Journal of the Royal Society Interface 9:339–350.
- 8. Blackshaw JK, Blackshaw AW (1994) Heat stress in cattle and the effect of shade on production and behaviour: a review. Australian Journal of Experimental Agriculture 34:285–295.
- 9. McLean RW, McCown RL, Little DA, Winter WH, Dance RA (1983) An analysis of cattle live-weight changes on tropical grass pasture during the dry and early wet seasons in northern Australia. 1. The nature of weight changes. Journal of Agricultural Science 101:17–24.
- 10. Wharton RH, Utech KBW, Turner HG (1970) Resistance to the cattle tick, Boophilus microplus in a herd of Australian Illawarra Shorthorn cattle: its assessment and heritability. Australian Journal of Agricultural Research 21:163–181.
- 11. Bock RE, Kingston TG, De Vos AJ (1999) Effect of breed of cattle on transmission rate and innate resistance to infection with Babesia bovis and B. bigemina transmitted by Boophilus microplus. Australian Veterinary Journal 77:461–464.
- 12. Peacock FM, Koger M, Olson TA, Crockett JR (1981) Additive genetic and heterosis effects in crosses among cattle breeds of British, European and Zebu origin. Journal of Animal Science 52:1007–1013.
- 13. Seifert GW (1971) Variations between and within breeds of cattle in resistance to field infestations of the cattle tick (Boophilus microplus). Australian Journal of Agricultural Research 22:159–168.
- 14.
Mayr E (1963) Animal Species and Evolution: Belknap Press of Harvard University Press, Cambridge, Massachusetts.
- 15.
Dobzhansky T (1952) Nature and origin of heterosis. In:Gowen JWeditor. Heterosis. Ames, Iowa: Iowa State College Press.
- 16. Barendse W, Harrison BE, Hawken RJ, Ferguson DM, Thompson JM, et al. (2007) Epistasis between calpain 1 and its inhibitor calpastatin within breeds of cattle. Genetics 176:2601–2610.
- 17. Hawken RJ, Zhang YD, Fortes MRS, Collis E, Barris WC, et al. (2012) Genome-wide association studies of female reproduction in tropically adapted beef cattle. Journal of Animal Science 90:1398–1410.
- 18. Turner LB, Harrison BE, Bunch RJ, Porto Neto LR, Li YT, et al. (2010) A genome wide association study of tick burden and milk composition in cattle. Animal Production Science 50:235–245.
- 19. Hindorff LA, Sethupathy P, Junkins HA, Ramos EM, Mehta JP, et al. (2009) Potential etiologic and functional implications of genome-wide association loci for human diseases and traits. Proceedings of the National Academy of Sciences of the United States of America 106:9362–9367.
- 20. Speliotes EK, Willer CJ, Berndt SI, Monda KL, Thorleifsson G, et al. (2010) Association analyses of 249,796 individuals reveal 18 new loci associated with body mass index. Nature Genetics 42:937–U953.
- 21. Allen HL, Estrada K, Lettre G, Berndt SI, Weedon MN, et al. (2010) Hundreds of variants clustered in genomic loci and biological pathways affect human height. Nature 467:832–838.
- 22. Sanders JO (1980) History and development of zebu cattle in the United States. Journal of Animal Science 50:1188–1200.
- 23. Gibbs RA, Taylor JF, Van Tassell CP, Barendse W, Eversole KA, et al. (2009) Genome-wide survey of SNP variation uncovers the genetic structure of cattle breeds. Science 324:528–532.
- 24. Bolormaa S, Pryce JE, Kemper K, Savin K, Hayes BJ, et al. (2013) Accuracy of prediction of genomic breeding values for residual feed intake and carcass and meat quality traits in Bos taurus, Bos indicus, and composite beef cattle. Journal of Animal Science 91:3088–3104.
- 25. Porto-Neto LR, Sonstegard TS, Liu GE, Bickhart DM, Silva MVB, et al. (2013) Genomic divergence of zebu and taurine cattle identified through high-density SNP genotyping. BMC Genomics 14:876.
- 26. Burrow HM (1998) The effects of inbreeding on productive and adaptive traits and temperament of tropical beef cattle. Livestock Production Science 55:227–243.
- 27. Sutherst RW, Bourne AS, Maywald GF, Seifert GW (2006) Prevalence, severity, and heritability of Stephanofilaria lesions on cattle in central and southern Queensland, Australia. Australian Journal of Agricultural Research 57:743–750.
- 28. Prayaga KC, Corbet NJ, Johnston DJ, Wolcott ML, Fordyce G, et al. (2009) Genetics of adaptive traits in heifers and their relationship to growth, pubertal and carcass traits in two tropical beef cattle genotypes. Animal Production Science 49:413–425.
- 29.
Olson TA (1999) Genetics of color variation. In:Fries R, Ruvinsky Aeditors.The genetics of cattle. Wallingford: CAB International. pp.33–53.
- 30. Porto-Neto LR, Kijas JW, Reverter A (2014) The extent of linkage disequilibrium in beef cattle breeds using high-density SNP genotypes. Genetics Selection Evolution 46:22.
- 31. Karim L, Takeda H, Lin L, Druet T, Arias JAC, et al. (2011) Variants modulating the expression of a chromosome domain encompassing PLAG1 influence bovine stature. Nature Genetics 43:405–413.
- 32. Mariasegaram M, Chase CC, Chaparro JX, Olson TA, Brenneman RA, et al. (2007) The slick hair coat locus maps to chromosome 20 in Senepol-derived cattle. Animal Genetics 38:54–59.
- 33. Fortes MRS, Kemper K, Sasazaki S, Reverter A, Pryce JE, et al. (2013) Evidence for pleiotropism and recent selection in the PLAG1 region in Australian Beef cattle. Animal Genetics 44:636–647.
- 34. Hudson NJ, Reverter A, Wang YH, Greenwood PL, Dalrymple BP (2009) Inferring the transcriptional landscape of bovine skeletal muscle by integrating co-expression networks. PLoS ONE 4:e7249.
- 35. Lehnert SA, Byrne KA, Reverter A, Nattrass GS, Greenwood PL, et al. (2006) Gene expression profiling of bovine skeletal muscle in response to and during recovery from chronic and severe undernutrition. Journal of Animal Science 84:3239–3250.
- 36.
Heffernan C (2013) Introduction: special issue on animal health and climate change. Animal Health Research Reviews.
- 37. Running SW (2012) A measurable planetary boundary for the Biosphere. Science 337:1458–1459.
- 38. Barwick SA, Johnston DJ, Burrow HM, Holroyd RG, Fordyce G, et al. (2009) Genetics of heifer performance in ‘wet’ and ‘dry’ seasons and their relationships with steer performance in two tropical beef genotypes. Animal Production Science 49:367–382.
- 39. Barendse W, Reverter A, Bunch RJ, Harrison BE, Barris W, et al. (2007) A validated whole genome association study of efficient food conversion. Genetics 176:1893–1905.
- 40. Matukumalli LK, Lawley CT, Schnabel RD, Taylor JF, Allan MF, et al. (2009) Development and characterization of a high density SNP genotyping assay for cattle. PLoS ONE 4:e5350.
- 41. Zimin AV, Delcher AL, Florea L, Kelley DR, Schatz MC, et al. (2009) A whole-genome assembly of the domestic cow, Bos taurus. Genome Biol 10:R42.
- 42. Browning SR, Browning BL (2010) High-Resolution Detection of Identity by Descent in Unrelated Individuals. American Journal of Human Genetics 86:526–539.
- 43. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, et al. (2007) PLINK: A tool set for whole-genome association and population-based linkage analyses. American Journal of Human Genetics 81:559–575.
- 44. Gibbs RA, Taylor JF, Van Tassell CP, Barendse W, Eversole KA, et al. (2009) Genome-wide survey of SNP variation uncovers the genetic structure of cattle breeds. Science 324:528–532.
- 45.
Akerman J (1982) American Brahman: a history of the American Brahman. Houston, Tex.: American Brahman Breeders Association. 384 p.
- 46.
Felius M (1985) Genus Bos: Cattle Breeds of the World. Rahway, NJ: Merck and Co., Inc.
- 47. Alexander DH, Novembre J, Lange K (2009) Fast model-based estimation of ancestry in unrelated individuals. Genome Research 19:1655–1664.
- 48. Perez-Enciso M, Misztal I (2011) Qxpak.5: Old mixed model solutions for new genomics problems. BMC Bioinformatics 12:202.
- 49.
Falconer DS, Mackay TFC (1996) Introduction to Quantitative Genetics. Harlow, Essex, UK: Longmans Green.
- 50. Bolormaa S, Pryce JE, Reverter A, Zhang Y, Barendse W, et al. (2014) A multi-trait, meta-analysis for detecting pleiotropic polymorphisms for stature, fatness and reproduction in beef cattle. PLoS Genetics 10:e1004198.