A genome-wide association study maps congenital gonadal hypoplasia to BTA29

Gonadal hypoplasia in Northern Finncattle and Swedish Mountain cattle is hypothesized to result from pleiotropic effects of a white colour gene because most of the affected animals are predominantly white [6,23]. To account for both the different coat colour patterns and the affection status of the animals, four different case-control cohorts were established (Table 1) and allelic associations were performed. P-values below 7.71 x 10^-8 (Bonferroni-corrected threshold for multiple testing) were considered to be significantly associated. A genome-wide association study (GWAS) with 20 unaffected predominantly black and brown animals and 21 affected predominantly white animals revealed 18 significantly associated SNPs on BTA29 between 17.69 Mb to 20.50 Mb (Figure 2A, Table S2). The SNP showing most significant association was BovineHD2900005672 (19,661,149 bp, P=2.19 x 10^-13). The second GWAS, including all available unaffected animals (N=73) and all affected animals, yielded five highly significantly associated SNPs on BTA29 and one highly significantly associated SNP on BTA23 at 52,435,290 bp, (P=9.57 x 10^-9) (Figure 2B). This SNP was considered to be a false positive result because there were no other significantly associated SNPs on BTA23. Moreover, re-genotyping the SNP using Sanger-sequencing revealed discrepant genotypes compared with the genotypes obtained with the Illumina BovineHD Bead chip in some samples. This indicated technical problems. Comparison of predominantly white -40 unaffected and 21 affected – animals revealed a suggestive association (Table S2). Similarly, GWAS in the fourth group (D) using 40 unaffected predominantly white animals and 20 unaffected predominantly coloured animals revealed a suggestive association (Table S2). However the SNPs of the last two GWAS analysis were not significantly associated after Bonferroni-correction (Figure 2C and 2D). Additionally, association analyses were repeated using a genotypic test. The results were in line with the allelic test, underlining the importance of the 17.69 Mb to 20.50 Mb area in BTA29 (Table S3). Collectively, the analysis of the four studied cohorts indicates that the BTA29 region is associated with gonadal hypoplasia and not only reflects different coat colour patterns. The genomic inflation factors for the groups ranged from 0.94 to 1.02 and indicate the absence of potential spurious associations.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Association of 647,971 SNPs in four different case control scenarios.

The Manhattan plots represent the -log₁₀(P) values of association for 647,971 SNPs in four different case-control designs (Table 1). The red dots represent significantly associated SNPs (P < 7.71 x 10^-8).

https://doi.org/10.1371/journal.pone.0075659.g002

CNV analysis identifies a translocated genomic region with the KIT gene

Durkin et al. [15] showed that colour sidedness in various cattle breeds is determined by two CNVs that duplicate a small part of the BTA6 and BTA29. Because of the strong link between white coloration and gonadal hypoplasia in Northern Finncattle and Swedish Mountain cattle, we carried out CNV analysis. A total of 2101 autosomal CNVs were identified in 94 animals. Among them a CNV segment on BTA6, extending from 71.60 Mb to 72.13 Mb and containing the KIT gene, was significantly associated with gonadal hypoplasia (P=5.65 x 10^-6) (Figure 3, Figure 4A). The CNV was present in all 21 affected animals and in 52 of 73 animals of the control group. Also, a CNV segment on BTA29, extending from 19.86 Mb to 20.32 Mb, appeared to be frequent in the Northern Finncattle. The CNV segment on BTA29 was identified in 14 affected and 44 unaffected animals, but was not associated with gonadal hypoplasia (Figure S2). This segment is in the immediate vicinity of the most significantly associated gonadal hypoplasia SNPs (Figure 4B). To ensure that CNV segment on BTA6 was not associated to gonadal hypoplasia nominal p-values of the GWAS for the SNPs located vicinity of the segment were studied (Tables S4 and S5). Only with the Design D in the genotypic test a few SNPs with suggestive association was found. This design included only unaffected predominantly white and coloured animals so the found associations were related to coat colour not gonadal hypoplasia.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Association of 2101 autosomal CNVs with the affection status of 94 animals.

The presence of CNV segments was compared in 21 cases and 73 controls using Fisher exact tests. The dots represent SNPs within the CNV segments. Red dots represent SNPs in significantly overrepresented CNVs in cases vs. controls (P < 2.38 x 10^-5).

https://doi.org/10.1371/journal.pone.0075659.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Schematic view of two CNV segments on BTA6 and BTA29.

The figures display 5-SNP-sliding-window log R ratios in 75 unaffected (grey) and 21 affected (blue) animals on BTA6 (A) and BTA29 (B). The grey shaded box represents the extent of the two CNV segments. The gene content was assessed based on the University of Maryland assembly.

https://doi.org/10.1371/journal.pone.0075659.g004

Durkin et al. [15] also showed that the two colour sidedness loci result from a translocation of a chromosomal segment, encompassing the KIT gene, from BTA6 to BTA29 and a part of this fragment back to BTA6 near the KIT gene. To investigate whether similar duplications are present in Northern Finncattle and the Swedish Mountain breed, we used 38 animals with copy number (CN) four on BTA6 as the case group and 21 animals with CN two as a control group to identify a potential translocation. The CN status was determined based on the array genotypes. We found a strong association on BTA29 (19,661,149 bp, P = 1.66 x 10^-28), indicating a translocation of the BTA6 segment to BTA29 (Figure S3). We next analysed 58 and 36 animals with and without CNV on BTA29, respectively, and identified a strong association on BTA6 (72,198,048 bp, P=1.69 x 10^-14, Figure S4), indicating a translocation from the BTA29 segment to BTA6. These positions agree with those described by Durkin et al. [15] indicating that the duplications identified are the same.

Cytogenetic studies confirm translocations between BTA6 and BTA29

Three Northern Finncattle animals were selected for cytogenetic investigations based on the CN status of the BTA6-segment. One of these animals was affected by gonadal hypoplasia. A Western Finncattle and an Eastern Finncattle animal were also included in the comparative analysis. These breeds are not affected by the gonadal hypoplasia. Traditionally, Western Finncattle is a solid brown breed and Eastern Finncattle is a brown breed with the colour-sided phenotype. Bacterial artificial chromosome (BAC) clones specific to the CN areas of chromosomes BTA6 (red FISH signal) and BTA29 (green FISH signal) were used as probes in FISH experiments. Overlapping red and green signals appear yellow. In the Western Finncattle animal two red (BTA6 homologs) and two green (BTA29 homologs) signals were recorded. The Eastern Finncattle animal showed the following pattern: two yellow signals on BTA6 and two green signals on BTA29. All three animals from the Northern Finncattle breed showed red and yellow signals on BTA6 while the patterns on BTA29 varied. The first unaffected Northern Finncattle showed green and yellow signals on BTA29 whereas the second unaffected Northern Finncattle showed only green signals. The affected Northern Finncattle animal had yellow signals on BTA29 (Figure S5). Taken together, all animals except the solid brown Western Finncattle animal had one or several translocations associated with colour sidedness. The affected Northern Finncattle was homozygous for the Cs₂₉ allele. This supports our GWAS and CNV results that showed strong association between the Cs₂₉ allele and gonadal hypoplasia.

CNV confirmation using quantitative PCR

The number of the Cs₂₉ alleles was analysed by quantitative PCR in nine affected and 21 unaffected animals. The results agree with the CNV results based on the array genotypes. However, using qPCR we could not distinguish between heterozygous and homozygous animals.

Confirmation of CNV using conventional PCR

To characterize the translocation present in Northern Finncattle and Swedish Mountain cattle we used PCR primers according to Durkin et al. [15] and one additional primer pair that flanked the insertion site of wild type BTA29. Analysed animals included also two animals that were omitted from the GWAS analyses because of the genotyping failure. Colour of these two animals was predominantly brown and white. The results showed that all affected animals were homozygous for the Cs₂₉ allele and 15 animals had the Cs₆ allele. Of the unaffected animals, 18 were homozygous and 37 were heterozygous for the Cs₂₉ allele and 20 did not have the Cs₂₉ allele. The Cs₆ allele was detected from 44 control animals (Table 2). The ratio of the Cs₆ allele is slightly higher in affected than in unaffected animals: 15/21 (71%) and 44/75 (59%) respectively. More importantly the occurrence of Cs₆ in control animals that were homozygous for the Cs₂₉ allele was 14/18 (78%) which is almost the same as in affected animals (Table 2). Like the case animals, all homozygous control animals were predominantly white. The ratio of the Cs₆ allele in white control animals that were heterozygous for the Cs₂₉ allele was 11/18 (61%) which is also quite close to the Cs₆ ratio of the case animals (Table S6). These results suggests that the homozygous Cs₂₉ allele might have pleiotropic effects on both coat colour and gonadal hypoplasia but the Cs₆ allele affects only the coat colour. The results are in line with the CNV results based on array genotypes and qPCR. Additionally, we determined the DNA sequence of all PCR products amplified from the selected 12 animals (results not shown). The sequences were comparable with those obtained by Durkin et al. [15].

A genome-wide association study comparing the case and control animals homozygous for the Cs29 allele revealed no association

Given that the Cs₂₉ allele was associated with gonadal hypoplasia we performed an additional conditional GWAS where the case animals were compared with 18 control animals that were homozygous for the Cs₂₉ allele. The only suggestively associated allele was the SNP on BTA23 at 52,435,290 bp which was shown to have discrepant genotypes (Figure S6).

Ethics statement

The blood sampling and clinical examinations were carried out following standard veterinary protocols in Finland. The Animal Ethics Committee of the State Provincial Office of Southern Finland approved all animal work (ESAVI-2010-03428/Ym-23).

Gonad examination and sampling

Clinical examinations were done during farm visits by experienced veterinarians. Bulls and bull calves were palpated for testicle size and symmetry. Female animals older than 16 months, excluding animals more than five months pregnant, were studied by ovarian palpation per rectum. After clinical examination, 9 ml of EDTA venous blood were collected from the animals suitable for this study. Semen samples of artificial insemination bulls included in the study were also collected. From the post-mortem study animals’ gonads were examined visually, palpated and weighed. Histological and tissue samples for DNA extraction were taken from the collected gonads. Histological samples were subjected to standard Bouin’s fixation and embedded in paraffin. Sections (5 µm) were cut and stained with haematoxylin-eosin (HE).

Animals considered affected by gonadal hypoplasia

Cows and heifers with one ovary of extremely small size were considered to be affected by gonadal hypoplasia. Usually the hypoplastic ovary was undetectable by palpation per rectum. Bulls and bull calves were considered to be affected by unilateral gonadal hypoplasia if their testicles clearly differed in size, i.e. one testicle being more than two times larger than the other testicle in young calves and for bulls’ one testicle being more than three times larger than the other testicle. Bilateral gonadal hypoplasia was diagnosed only in one bull selected for semen collection. Both testicles were very small and there were no sperm cells in the ejaculate. The testicle histology only showed Sertoli cells and no spermatogonia in seminiferous tubules in the hypoplastic testicles. When possible, the clinically diagnosed affected animals were re-examined during a second farm visit or the gonads were studied after slaughter.

DNA extraction

Genomic DNA from blood samples was extracted by automatic isolation (Magnetic Separation Module I, Chemagen) and genomic DNA from tissue and semen samples was extracted using a commercially available Qiagen Kit (QIAamp DNA Mini Kit). Extractions of the blood and tissue samples were made according to the manufacturer’s instructions. The extraction of DNA from semen samples was made according to DNA Purification from Tissues-protocol in QIAamp DNA Mini Kit handbook with some modifications. From 200 to 500 µl of frozen semen was centrifuged for 5 min at 100 × g. The supernatant was moved and the pellet was washed with 200 µl of phosphate buffered saline (PBS). The Qiagen buffer ALT was added up to 300 µl and 20 µl of Proteinase K and dithiothreitol was also added. The mix was incubated for 1 h at 56 °C. During incubation the sample was pulse vortexed four times for 15 sec. 300 µl of Qiagen buffer AL was added and the pulse vortex repeated. The sample was incubated for 10 min at 56 °C and thereafter 150 µl of 96% alcohol was added. The sample was pulse vortexed and incubated for 3 min at room temperature. The whole mixture was applied to the QIAamp Mini spin column and centrifuged at 6,000 × g for 1 min. The filtrate was discarded and the column was washed twice with 500 µl of Qiagen buffer AW1 and once with 500 µl of Qiagen buffer AW2 (centrifugation at 6,000 × g for 1 min). The column was dried with centrifugation at 20,000 × g for 3 min. The DNA was eluted with 50 µl of distilled water, incubated for 1 min at room temperature and centrifuged at 20,000 × g for 1 min.

Animals selected for genotyping

DNA samples from 96 animals were included in this study. Animals with total unilateral or bilateral hypoplasia were included in the case group that comprised 21 animals (10 males and 11 females). Of these animals, one was affected with bilateral, 18 with left-sided and two with right-sided gonadal hypoplasia. The control group included 75 unaffected animals. Among the studied animals 91 were Northern Finncattle and five were Swedish Mountain cattle (two cases and three controls).

High-density genotyping and quality control

Ninety-six animals (21 affected / 75 unaffected) were genotyped with the Illumina BovineHD Bead chip interrogating genotypes of 777,962 SNPs. Genotype calling was performed using default parameters of Illumina’s BeadStudio. Quality control was carried out with PLINK v1.07 [28]. We excluded 1224, 343 and 1735 SNPs with Y-chromosomal, mitochondrial and unknown chromosomal position, respectively, for further analysis. The genotypes of two unaffected animals were omitted because genotyping failed in more than 10% of the SNPs. We further excluded 6229 SNPs because genotyping failed in more 10% of the individuals and 121,657 monomorphic SNPs. The final dataset comprised 94 animals and 647,971 SNPs with an average call-rate of 99.67%. The chromosomal position of the SNPs was determined based on the UMD3.1 assembly of the bovine genome [29].

Genome-wide association study

Fisher exact tests of allelic and genotypic associations were performed to compare genotypes in cases vs. controls at each SNP in turn using PLINK [28]. We considered SNPs with P < 7.71 x 10^-8 as significantly associated (Bonferroni-corrected threshold for multiple testing). Quantile-quantile plots were inspected and genomic inflation factors were calculated according to Devlin and Roeder [30] to assess the extent of false positive association signals.

Detection of copy number variants

Genotype signal intensities obtained from genotyping with the Illumina BovineHD Bead chip (see above) were analysed with PennCNV [31] to identify copy number variations (CNV). Briefly, the implemented CNV-detection algorithm considers both the log R ratio (LRR) and the B allele frequency (BAF), as well as the allele frequency and the distance of adjacent SNPs. Two individuals and 6229 SNPs with poor genotyping quality (see above) were not considered for the identification of CNVs. The presence of CNVs containing a minimum number of 10 SNPs corresponding to a minimum length of approximately 35 Kb was compared in cases vs. controls using Fisher exact tests.

Cytogenetic analysis

Metaphase spreads were obtained from short-term lymphocyte cultures, established for five animals. The FISH study was carried out according to a protocol described by Durkin et al. [15]. Two BAC clones (RP42-160M9, RP42-156I13) covering the region of interest on chromosome BTA6 (region 72,566,605–72,817,995bp, UMD 3.0) and two BAC clones (RP42-37P11, RP42-116G8) covering the region of interest on BTA29 (region 20,772,406–21,035,251bp, UMD 3.0) were derived from the RPCI-42 Bovine BAC Library (http://bacpac.chori.org/home.htm). BAC DNA was isolated using an alkaline lysis method and labelled by random priming. The isolated DNAs from two BAC clones specific to BTA6 were mixed equally and labelled using biotin-11-dUTP, while DNA from two BAC clones specific for BTA29 was also mixed and labelled with digoxigenin-11-dUTP. The labelled probes with an excess of bovine Cot-1 DNA were separately denatured for 10 min at 70°C and applied on denatured chromosome slides. Hybridization was carried out overnight at 37°C. After slide washing, biotin-labelled probes were detected using streptavidin-Cy3 (Amersham, 1:200, red colour) and digoxigenin-labelled probes were detected with antidigoxigenin-fluorescein Fab fragments (Roche, 1:200, green colour). The slides were counterstained with Vectashield containing DAPI (Vector Laboratories) and examined with an epifluorescence Nikon E600 Eclipse microscope equipped with a cooled digital CCD camera and Lucia software.

QPCR

QPCR was used to validate CNV identified after bioinformatics analysis. Altogether 30 samples were analysed: nine cases and 21 controls. A simple method based on Weksber et al. [32] and Lachman et al. [33] was applied for analysis. Two primer pairs designed using Primer 3 [34] were located in the KIT region: precisely in exon 8 (GGGCCAGTGGATGTACAGAT) and 9 (TGCAAAGTTAAAAGAGGCAGA); the second pair in exon 18 (CACATTTGAAAGTGATGTCTGG) and 19 (AGAACTTAGAATCGACTGGCATT). 10 µl of QPCR reaction was composed from the Fast SYBR® Green Master Mix (Life Technologies) and 2.5 pmol of each primer. All amplifications were carried out in Applied Biosystems 7500 Fast Real-Time PCR System (Life Technologies) according to manufacturer’s recommendations. H6PD was used as a reference gene; both primers were located in the first exon (AAGGTCCTGGAGTCCCTGTC, GTAGAAAATTCGGCCGGTCT).

PCR and sequencing

All animals were analysed with standard PCR carried out using breakpoint primers designed by Durkin et al. (Table 2) [15] and one additional primer pair PSK_α-β2_R (TGGGTAGACAGGTTTGTTTCC and TCTTGACCACTTGCATTGGA) that flanked the insertion site of the wild type BTA29. A PCR reaction of 20 µl volume containing 20 ng genomic DNA, 1x Qiagen PCR buffer, 1.5 mM MgCl2, 200 µM of each nucleotide, 5 pmol of each forward and reverse primer (Sigma) and 0.5 units of Taq Polymerase (Qiagen) was performed under the following conditions: initial denaturing at 95 °C for 3 min, followed by 35 cycles at 94 °C for 30 sec, 60 °C for 1 min, 72 °C for one minute and the final extension at 72 °C for 3 min. PCR products were visualized together with GeneRulerTM100bp DNA Ladder (Fermentas, Thermo Scientific) on 1.5% agarose gel.

Five case samples, two control samples with a duplication of the BTA6 segment, two control samples with both studied duplications and three control samples without the duplications were subjected to sequencing. A 10-20 ng of purified PCR product was mixed with 0.5 µl BigDye® Terminator v3.1 (Life Technologies) and 0.5 µl of the forward or reverse PCR primer (2.5 pmol). The sequencing reaction was performed under the following conditions: initial denaturing at 96 °C for 10 sec, followed by 35 cycles, including 10 sec in 96 °C, 5 sec in 50 °C, and 4 min in 60 °C. The gel filtration of the sequencing reaction was applied using the MultiScreen filtration plate (MAHVN4510; Millipore) and Sephadex G-50 Fine (Sigma) followed by capillary electrophoresis carried out on 3130xl Genetic Analyzer (Life Technologies). Base calling, sequence alignment and polymorphism detection were made using the Phred/Phrap/Polyphred software [35-37]. Sequences were inspected using Consed [38].

Seven cases and ten controls were sequenced (see above) with primers CAATTTTAAGCATGTGCTGAGG and ACAGCCTCTGGTCTGTCTGG to validate the genotype calls for the SNP on BTA23 at 52,435,290 bp.