Finding genes for economically important traits: Brahman cattle puberty

Additional keywords: Bos indicus, corpus luteum, genome-wide association, scrotal circumference.

Age at puberty is an important component of cattle performance. It determines the beginning of reproductive life for any breeding animal and it influences generation interval, affecting the rate of genetic improvement and herd productivity. Selecting for early pubertal heifers and bulls is a practical approach for reducing the generation interval and potentially increasing fertility (Lesmeister et al. 1973; Foster 1994; Siddiqui et al. 2008; Johnston et al. 2010).

For genetic selection, puberty must be defined as a measurable and inherited trait. Puberty in heifers was defined by plasma progesterone concentration, ultrasound images of corpus luteum (CL) and detection of oestrous (Shamay et al. 2005; Romano et al. 2007; Johnston et al. 2009). Puberty in bulls was defined by scrotal circumference thresholds and was related to sperm concentration, motility and morphology (Lunstra et al. 1978; Bagu et al. 2004; Siddiqui et al. 2008; Corbet et al. 2009). Pubertal traits vary in heritability (h²) from low to moderately high (Cammack et al. 2009).

There is a historical and increasing body of evidence reporting the genetic correlations between female and male puberty in cattle (Martin et al. 1992). Scrotal circumference of bulls correlates with puberty in their female relatives. Correlation estimates range from r = −0.71 to r = −0.15 (Burns et al. 2011). This evidence supports the indirect selection strategy for early pubertal heifers, using correlated traits measured in bulls. Given the size of the genetic correlation between these traits, some genes and genetic markers associated with puberty in heifers are likely to be associated with puberty in bulls, and vice versa.

Bos indicus cattle, such as Brahman cattle, are reportedly older at puberty when compared with most Bos taurus breeds (Lunstra and Cundiff 2003; Lopez et al. 2006). It is of economic relevance to select for early age at puberty in Brahman cattle. To assist selection and identify genes underlying fertility, we performed a GWAS in a sample of bulls with a measurement for puberty, namely the age at a scrotal circumference of 26 cm (AGE26). The present study is the first report of a GWAS regarding scrotal circumference in bulls, while the GWAS on AGECL in heifers represents a re-analysis of a previously reported study (Hawken et al. 2012). The Brahman heifers measured for AGECL were the mothers of the bulls measured for AGE26. Heritability, genetic correlation and genes associated with AGECL and AGE26 are reported and discussed.

Cattle were bred and supplied by the Cooperative Research Centre for Beef Genetic Technologies (Beef CRC) as described previously (Corbet et al. 2009; Johnston et al. 2009). In brief, data from 1007 Brahman heifers and 1118 Brahman bulls were included in the current analysis. The phenotypes of the heifers have been reported previously (Johnston et al. 2009). The bulls included in the present experiment were the offspring of the heifers, which became breeding heifers for the Beef CRC. Their pedigree consists of over 50 grandsire families; 54 bulls that sired heifers, which were mated to 55 industry sires to produce the bulls used in the study.

The traits that defined puberty were AGECL in heifers and AGE26 in bulls. The first CL was detected by ovarian ultrasound examinations conducted every 4–6 weeks, after heifers had reached an average liveweight of 200 kg (Johnston et al. 2009). AGECL was estimated from annotation of the date when the first CL was observed and the date of birth. Scrotal circumference (SC) was measured with a standard metal tape (Fordyce et al. 2006). Between weaning and 24 months of age, eight measurements of SC were taken for each bull, at 3-month intervals. Summary statistics for the age and SC of bulls at each of the eight time points are presented in Table 1. Using these repeated measurements for individual regressions, we interpolated the age when the bull achieved 26 cm of SC (AGE26, expressed in days). Achieving SC of 26 cm was considered a threshold for puberty in Bos indicus bulls. The threshold was different from the previously described 28-cm threshold because Bos indicus typically present a more elongated scrotum and a smaller SC than do Bos taurus (Lunstra et al. 1988; McGowan et al. 2002; Silva et al. 2011). For a visual assessment of age v. SC in this population, see Fig. 1. A more in-depth description and quantitative analysis of all reproductive and growth phenotypes available for this population of bulls was reported elsewhere (Corbet et al. 2009, 2011).

Table 1.  Means ± s.e. for measurements of scrotal circumference (SC) at eight time points

Fig. 1.  Age and scrotal circumference (SC) measured across eight time points. Scatter plot of the age in days (x-axis) versus the SC in cm (y-axis), with time points represented by alternating colours (black and white).

The Illumina BovineSNP50 array (Van Tassell et al. 2008; Matukumalli et al. 2009) was used to genotype samples according to the manufacturer’s protocols (Illumina Inc., San Diego, CA). Positions for each SNP were based on the UMD3 assembly of the bovine genome sequence (available from http://www.livestockgenomics.csiro.au/cgi-bin/gbrowse/btauUMD3/, verified 14 December 2011). Genotypes for 843 heifers had been generated previously using the BovineSNP50 v1 array and were reported in two recent GWAS that analysed growth traits and reproduction traits, including AGECL (Bolormaa et al. 2011; Hawken et al. 2012).

In the present study, 1118 bulls were genotyped using the Illumina BovineSNP50 v2 array. Quality control was similar to the previous study, with repeat samples included and the Bead Studio software from Illumina (www.illumina.com, verified 28 February 2012) used to determine the genotype calls. Animals with call rates <98% were excluded, resulting in 1105 bulls being retained for analyses. SNPs with auto-calling rates <85% and SNPs with a minor allele frequency <0.01 were excluded from further analyses.

For the current study, missing genotypes for bulls and heifers were imputed using the BEAGLE 3.2 program (Browning and Browning 2010). Imputation and quality-control procedures yielded a total of 43 821 SNPs used in the genome-wide association analysis. Genotype calls were coded as 0 for the homozygote of the first allele (A), 1 for the heterozygote, and 2 for the homozygote of the second allele (B). Alleles A and B were defined as per top and bottom rules from Illumina (http://www.illumina.com/documents/products/technotes/technote_topbot.pdf, verified 17 February 2012). Imputation resulted in better coverage of the X chromosome, adding new data to the previous heifer study on AGECL (Hawken et al. 2012). This improved X chromosome coverage for the heifers and all genotypes obtained for bulls are reported here for the first time.

Genetic correlation and h² (for AGE26) were estimated fitting the two traits in a bivariate analysis using a mixed model and three generation pedigree relationships. For AGECL, fixed effects included in the model and h² have been reported previously (Johnston et al. 2009); here, we estimated the h² de novo in the bivariate analysis, using the same fixed effects. Fixed effects included in the model for AGE26 were contemporary group (bulls born in the same year and location – or cohort – and post-weaning location) and birth month. Solutions to the effects in the model as well as variance components were estimated using VCE v.6 software (Groeneveld and García-Cortés; http://vce.tzv.fal.de, verified 20 February 2012).

Genome-wide association studies (GWAS) were performed for AGECL and AGE26 separately and by using two approaches. First, GWAS was performed with the full datasets for bulls (1105) and heifers (843) and, subsequently, each population was split into calibration and a validation datasets, for providing additional evidence of SNP association and estimated effects. The calibration and validation datasets were generated according to birth date, dividing the populations into a younger half (~50%) and an older half (~50%). For AGECL, the calibration dataset was formed by heifers born until November 2001 (n = 480) and older half was born from December 2001 onwards (n = 363). For AGE26, the calibration dataset comprised those bulls born until December 2005 (n = 565) and the bulls in the validation dataset were born from January 2006 onwards (n = 550). Using birth date to define calibration and validation datasets is a tested approach (Hayes et al. 2010). A GWAS for AGECL using the full dataset has been reported previously (Hawken et al. 2012), but the results reserving 363 heifers for a validation exercise are reported only here. No previous study exists to relate the GWAS for AGE26 reported here.

The effect of each SNP was estimated three times, as follows: first, using phenotypic data from the full dataset; second, using phenotypic data from the calibration dataset; and finally, using phenotypic data from the validation dataset. In all three instances, the effect of each SNP was estimated using the mixed model of Eqn 1, as follows:

where y_ij represents the vector of observations from the ith heifer or bull at the jth trait (j = 1 and 2 for AGECL and AGE26, respectively); X is the incidence matrix relating fixed effects in β with observations in y_ij; Z is the incidence matrix relating random additive polygenic effects in u with observations in y_ij; S_k in the vector of genotypes for the kth SNP across all animals; s_jk represents the additive association of the kth SNP on the jth trait; and e_ij is the vector of random residual effects. Fixed effects included in the model were the same as described for the models used to estimate the genetic correlations and h². Tests for the significance of SNP association were conducted using Qxpak5 (Perez-Enciso and Misztal 2004, 2011). Qxpak5 performs a likelihood ratio test, testing the model with SNP versus the model without the SNP, against a chi-square distribution with 1 degree of freedom and this was carried out one SNP at a time. Chromosome X genotypes of males and females were not analysed together, because the traits were analysed separately.

False discovery rates were estimated using Eqn 2, as follows:

where n represents the total number of SNPs included in the study (in the present study, n = 43 821), P is the P-value threshold being used and k is the actual number of associated SNPs in the given P-value threshold.

The percentage of the genetic variance accounted for by the ith SNP was computed according to the formulae of Eqn 3, as follows:

where p_i and q_i are the allele frequencies for the ith SNP estimated across the entire population, a_i is the estimated additive effect of the ith SNP on the trait in question (AGECL or AGE26), and σ_g² is the REML estimate of the (poly-)genetic variance for the trait in question.

Descriptive statistics, h² and associated errors for AGECL and AGE26 are presented in Table 2. At puberty, heifers were on average older (751 ± 142 days) than bulls (555 ± 101 days). The heritability estimated for puberty in heifers was lower (AGECL, h² = 0.56 ± 0.11) than that in bulls (AGE26, h² = 0.78 ± 0.10). The traits AGECL and AGE26 were genetically correlated, with r = 0. 20 ± 0.10 estimated using REML methods.

Table 2.  Summary statistics and heritability estimates for age at the first corpus luteum (AGECL) and age at the scrotal circumference of 26 cm (AGE26)
Std, standard deviation; AGECL phenotype is reported in Johnston et al. (2009)

Figure 2 provides an overview of GWAS results for the full datasets presented as Manhattan plots (i.e. significance on the y-axis v. genome map position on the x-axis) for AGE26 and AGECL. Two regions showed clear association peaks in the GWAS of AGE26, one on the X chromosome and another on BTA14. The GWAS of AGECL also showed a peak on BTA14 for the same region where SNPs were significant for AGE26. There were 134 SNPs associated with AGECL and 146 SNPs associated with AGE26 (P < 0.0001, false discovery rate = 0.03). Of these SNPs, 32 (~22%) were associated with both traits (P < 0.0001, Table 3). These 32 SNPs that were in common for both traits were all located on BTA14, between 21.95 Mb and 28.4 Mb. The most significant SNP for AGECL was Hapmap23509-BTC-073113 at Position 27 198 715 of BTA14 (P = 1.36E–09). The most significant SNP for AGE26 was BTA-30242-no-rs at Position 85 495 447 of Chromosome X (P = 5.35E–13). The effects, P-values and proportion of genetic variance explained for each of the 32 SNPs associated (P < 0.0001) with both AGECL and AGE26 are presented in Table 4. The effects, P-values and proportion of the variance explained for the SNPs underlying the association peak located on X chromosome are also reported (Table 5).

Fig. 2.  Manhattan plots for age at a scrotal circumference of 26 cm (AGE26) and age at the first corpus luteum (AGECL). Plots represent the associations between genome-wide singe-nucleotide polymorphisms (SNPs) and AGE26, as well as, SNPs and AGECL. In both plots, chromosomal positions are in the x-axis and –log (P-values) are in the y-axis.

Table 3.  Number of single-nucleotide polymorphisms (SNPs) associated with age at the first corpus luteum (AGECL), age at the scrotal circumference of 26 cm (AGE26) and both traits at various P-value thresholds
FDR, false discovery rate

Table 4.  Effect, P-values and proportion of the variance explained for single-nucleotide polymorphisms (SNPs) associated (P < 0.0001) with both age at the first corpus luteum (AGECL) and age at the scrotal circumference of 26 cm (AGE26) and located in chromosome 14
SNP positions and distance to nearest gene are reported in base pairs and are based on the UMD3 assembly of the bovine genome sequence. Most significant SNPs in terms of P-values are reported in bold. Note that SNP BTA-97369-no-rs was the only SNP not associated with both calibration and validation dataset for neither AGE26 nor AGECL (see Supplementary Material Tables S1, S2). MAF, minor allele frequency

Table 5.  Effect, P-values and proportion of the variance explained for single-nucleotide polymorphisms (SNPs) underlying the association peak in the X chromosome observed for age at the scrotal circumference of 26 cm (AGE26)
SNP positions and distance to nearest gene are reported in base pairs and based on the UMD3 assembly of the bovine genome sequence. Most significant SNP in terms of P-values are shown in bold. Note that all these SNPs were associated (P < 0.01) in both calibration and validation datasets for AGE26, with one exception, Hapmap32639-BTA-158837 (see Supplementary Material Tables S1, S2). MAF, minor allele frequency

The nearest annotated genes and their distance to the reported SNPs were determined using the UMD3.0 bovine genome assembly (Tables 4, 5). However, to consider only the genes nearest to each SNP could be misleading. It is preferable to consider all genes mapped to the association regions. In total, 23 genes mapped to the association region in Chromosome 14, including ATP6V1H, CA8, CHD7, CSF2RA, FAM110B, IMPAD1, NPBWR1, NSMAF, OPRK1, PCMTD1, PENK, POLR2K, RAB2A, RB1CC1, RGS20, RLBP1 L1, RP1, SDR16C6, SNTG1, SOX17, TGS1, TOX and XKR4. Also, 46 genes mapped to the association region in Chromosome X, including ABCB7, AKAP4, APOL, AR, ARR3, BRWD3, CACNA1F, CCDC120, CHIC1, CHM, CLCN5, CXCR3, CYLC1, DGAT2 L3, DGAT2 L6, EDA, EFNB1, FGF16, FOXP3, GPKOW, H11CXORF26, KIF4A, LRCH2, MGC140080, MGC152340, MIR374A, MIR374B, NHSL2, PAGE4, PJA1, PLS3, POU3F4, RBM10, SLC9A7, SNX12, STARD8, SUV39H1, TAF1, TAF9B, TIMM17B, TMEM28, UPRT, USP11, USP27X, ZDHHC15 and ZNF182.

Manhattan plots for the calibration and validation datasets are provided in the Supplementary Material Fig. S1. There were 111 SNPs that were associated (P < 0.01) with AGE26 in both calibration and validation datasets (Supplementary Material Table S1). There were 98 SNPs that were associated (P < 0.01) with AGECL in both calibration and validation datasets (Supplementary Material Table S2). Chromosome X harboured more SNPs associated with AGE26, in both calibration and validation sets, than did any other chromosome (49%, or 55 of 111 SNPs), followed by BTA 14 (20%, or 23 of 111 SNPs). Most (58%, or 57 of 98 SNPs) of the SNPs associated with AGECL in both calibration and validation sets were located at BTA 14. As a result of our calibration and validation exercise, the same regions and candidate genes as observed in the GWAS with the full dataset remained relevant.

Heritability estimates of cattle puberty in the literature are variable and influenced by population, breed, management and environmental factors, as well as the use of different phenotypic measurements (Martin et al. 1992; Cammack et al. 2009). Heritability for age at puberty in the Beef CRC Brahman heifers (AGECL, h² = 0.57) has been reported previously (Johnston et al. 2009); we confirmed these results by estimating a very similar value (h² = 0.56). These heritabilities are within the range, from 0.20 to 0.67, of other estimates for age at puberty (Arije and Wiltbank 1971; Smith et al. 1976; Laster et al. 1979). Reported heritability estimates for scrotal circumference range from 0.39 to 0.75 (Lunstra et al. 1988; Martinez-Velazquez et al. 2003; Corbet et al. 2009). Our result for puberty in bulls (AGE26) was closer to the higher end of that range (h² = 0.78). Overall and in agreement with our results, the heritability of puberty in cattle ranges from moderate to high, allowing for improvement through genetic selection.

Pubertal traits measured in heifers (AGECL) and bulls (AGE26) had a positive and favourable genetic correlation. This result confirms previous reports that showed that bigger scrotal circumference in bulls was correlated with early puberty in their female relatives (Martin et al. 1992; Burns et al. 2011). Previous evidence and our results point to the practicality of selecting for early pubertal cattle because both heifers and bulls can be selected with beneficial correlated effects.

Considering the genetic correlation, genes and SNPs associated with puberty in heifers were likely to be relevant for puberty in bulls, and vice versa. On BTA14, a large number of SNPs associated with puberty were identified in both bulls and heifers. The SNPs on BTA14 are located close to RP1, XKR4 and TOX. These genes are not known to affect puberty in any species. However, the region on BTA14 includes another 20 annotated genes and presents homology to Human chromosome HSA8q12, which is associated with height (Gudbjartsson et al. 2008). Among those 20 annotated genes, PENK, RPS20, SNORD54 and PLAG1 are plausible functional candidates. PENK has a role in GnRH regulation (Rosie et al. 1992; Taylor et al. 2007). RPS20 and a small RNA (SNORD54) were recently associated with calving ease (Pausch et al. 2011). The alleles that lower calving-ease have a positive effect on growth traits (Pausch et al. 2011). Calving ease is a trait influenced by the size of the calf and so, similarly to stature, it reflects frame size. Cattle with smaller frame size achieve puberty earlier than those with larger frame size (Vargas et al. 1999). Taken together, this evidence could point to RPS20 and SNORD54 as pleiotropic genes affecting puberty and calving-ease mainly through their effect on growth. Another recent study argued that PLAG1 was the relevant gene underlying this region on BTA14 and affecting bovine stature (Karim et al. 2011). Considering the physiological link between growth and puberty, it is possible to hypothesise that the association region in BTA14 may relate to a gene with many pleiotropic effects or a functional polymorphism that affects more than one gene. There is already evidence for a functional SNP in this region affecting the expression levels of multiple genes (Karim et al. 2011).

The two most significant SNPs for AGE26 were located on the X chromosome and not on BTA14. The genes nearest to the most significant SNPs, EDA and DGAT2 L3, do not appear to play a role in reproduction, given the current evidence from literature. A formal discussion of the functional link between these genes and cattle puberty is beyond our present objective. However, the association peak on Chromosome X does include the following three genes that are positional and functional candidate genes: AR (androgen receptor), TAF1 and TAF9B (both TATA box binding protein (TBP)-associated factors). These genes are candidate genes for several reasons. The androgen receptor is known to play a role in sexual development, specifically, it affects spermatogenesis, testis localisation and testis size in mice models (Verhoeven et al. 2010). In pigs, the AR is considered a candidate gene for reproduction and performance traits (Trakooljul et al. 2004). Manipulation of fetal androgen exposure alters the timing of puberty in sheep (Jackson et al. 2008). The genes TAF1 and TAF9B encode transcription factors that form the TFIID complex, a regulator of cell cycle and differentiation. Testis-specific TAF proteins have been reported as relevant for spermatid differentiation in Drosophila (Hiller et al. 2004). Taken together, this evidence and our results support AR, TAF1 and TAF9B as candidate genes for scrotal development and puberty in Brahman bulls.

Results of the present study are consistent with the hypothesis that the genes underlying the associated regions on BTA14 and Chromosome X play significant roles on defining age at puberty in Brahman cattle. Candidate genes mapping to these regions are the subject of ongoing research. Information about specific genes and markers will add value to genomic selection. Therefore, identifying these associated regions contributing to genetic variation in AGECL and AGE26 should assist in the selection for early onset of puberty in Brahman cattle.