Collection of phenotypic data
The original data consisted of 79,617 reproductive females and 785,993 mating records from July 1993 to March 2011. Pedigree information from reproductive females was traced back two generations. Embryonic mortality was treated as a categorical trait of reproductive females. Embryonic mortality was defined by cows that received a 2nd round of AI at 18–29 days (D), 30–60 D, 61–90 D, 91–140 D, and 141 D–parturition after the 1st AI. The threshold statistical model for embryonic mortality included fixed effects of parity, breeding farm, combination of mating year and month, and age at mating as covariance:
$$ {U}_{ijkl}=\mu + PARIT{Y}_i+ FAR{M}_j+Y{M}_k+b\ast AG{E}_{ijkl}+{a}_{ijkl}+{e}_{ijkl} $$
$$ {Y}_{ijklm}=0\left({U}_{ijklm}\le t\right),{Y}_{ijklm}=1\left({U}_{ijklm}>t\right) $$
where U
ijkl
is an unobservable underlying continuous variable (liability), μ is the overall mean, PARITY
i
is the fixed effects of parity, FARM is the fixed effects of farm, YM
k
is the fixed effects of month and year of AI, AGE
ijkl
is age at AI of observation
ijk
, b is the partial regression coefficient of age at AI, a
ijkl
is the direct additive genetic effects as random effects, e
ijkl
is the residual effect. Y
ijkl
is the observation
ijkl
for embryonic mortality; the trait was treated as a categorical trait: if the cow was a non-recipient of a 2nd AI in the defined period, it was allocated the value 0, if the cow was a recipient of a 2nd AI in the defined period, it was allocated the value 1. The observed binary response takes the value 1 and 0 if the liability is above a fixed threshold (t) and below t, respectively. Genetic parameters were estimated via Gibbs sampling procedure. The THRGIBBS1F90 program [13] was used to fit a threshold animal model for embryonic mortality. This program uses a Bayesian threshold model with Gibbs sampling to evaluate the posterior density of (co)variance estimates. The Gibbs analysis was run as a single chain of 110,000 cycles. After burn-in of the first 10,000 samples, every 100th sample was stored to estimate posterior means and standard deviations using the POSTGIBBSF90 program [13]. This program summarizes Gibbs samples obtained from THRGIBBS1F90 [12].
SNP genotyping and CNVR detection
Whole blood was collected from each cow, and genomic DNA was isolated using the Easy-DNA kit (Invitrogen, Cat. #K1800-01). All sample were genotyped using the using the Illumina BovineHD BeadChip Array (Illumina, Cat. #WG-450-1002), which contains 735,293 autosomal SNPs [42], according to the manufacturer’s instructions. SNP clustering and genotype calling was performed using GenomeStudio V2011 (Illumina, version 1.9.4), and all markers passed quality control (call rate >98%). The UMD3.1 assembly was used to map SNPs position [43]. CNVR were detected as previouly reported [9], detected using PennCNV software (version June 2011) [44, 45], which incorporates factors including log R ratio (LRR), B allele frequency (BAF), marker distance, and the population frequency of B allele (PFB) into a hidden Markov model.
Genome-wide CNV association study for embryonic mortality
Among 861 CNVRs, we selected 116 deleted-type CNVRs (Additional file 1: Table S1), which were a minor allele frequency > 0.01. The association analysis was performed on 791 samples using GEMMA software based on a linear mixed model with genomic relationships [21], which used a genetic relationship matrix estimated from SNP genotypes from the BovineHD BeadChip Array to model the correlation between the phenotypes of the sample subjects. Additionally, single-marker association analyses were performed for 791 samples using the SNPs in the BovineHD BeadChip Array.
qPCR validation of CNVR
Real-time qPCR was performed for CNVR validation using a 7900HT Real-Time PCR system (Applied Biosystems). Primers and probes were designed for CNVR_322 (Additional file 1: Table S8). Amplification reactions (20 μl⋅well−1) were carried out in triplicate with 20 ng of genomic DNA, 1× Absolute QPCR ROX Mix (Thermo Scientific, Cat. #AB-1138/B), 400 nM of each primer, and 200 nM of each probe. The basic transcription factor 3 gene (BTF3), which served as an internal qPCR standard for both copies at a locus (2n) [46], was co-amplified with the primers (Additional file 1: Table S9). Three replicate reactions were performed for each primer pair, and a comparative CT method was used to calculate the copy number [46]. ∆ CT was calculated by subtracting the BTF3 CT value from the sample CT value for each replicate. The average ∆ CT value for the three replicates was calculated. To determine the ∆∆ CT, the average ∆ CT of a calibrator animal, which had two copies of the DNA segment, was used. Finally, the copy number was given using the formula 2 × 2 -∆∆ CT.
Identification of the deletion breakpoint
We designed the forward primer series for 374,581 to 379,470 bp and the reverse primer series for 412,204 to 414,031 bp (Additional file 1: Table S4). The 800-bp PCR products, detected using CF_4 (377,383 to 377,402 bp) and TR_1 (412,204 to 412,223 bp), from the genomic DNA of CNVR_322-detected animals were subcloned into a pCR2.1-TOPO vector (Invitrogen, Cat. #K4500-01). Sixty-three bp on both sides of 33,934 bp flanking regions (378,065 to 378,127 bp and 411,999 to 412,061 bp were based on UMD3.1 assembly) in wild-type animals were amplified using primer pairs (Additional file 1: Table S4) and subcloned into a pCR2.1-TOPO vector. The plasmids were sequenced using M13_R and M13 (-20) primers by the BigDye Terminator v.3.1 Cycle Sequencing Kits (Applied Biosystems) and purified using CleanSEQ (Agencourt, Cat. #A29154) were sequenced, followed by electrophoresis using an ABI 3730 sequencer (Applied Biosystems).
DNA test for CNVR_322 with multiplex PCR
Multiplex PCR was performed using F_1, R_1 and R_2 primer (Additional file 1: Table S10). The PCR conditions were as follows: initial denaturation was carried out at 98 °C for 5 min, and followed by 35 cycles, each consisting of denaturation at 98 °C for 20 s, annealing at 60 °C and extension at 72 °C for 1 min. The amplified products were separated by electrophoresis on 1.5% agarose gels with 100 bp ladder markers (NEB, Cat. #N3231L).
Expression analysis using qPCR
For real time qPCR, we extracted total RNA from cow and mice tissues and dermal fibroblasts using RNeasy mini kits (QIAGEN, Cat. #74104) and treated it with DNase I. The cDNA was synthesized from 50 ng RNA using ReverTra Ace-α (TOYOBO, Cat. #FSK-101) with random primers, according to the manufacturer’s instructions. Bovine ANXA10 and mouse Anxa10 was amplified with the primers in described in Additional file 1: Table S11. Real-time PCR was performed on a 7900HT Real-Time PCR system (Applied Biosystems) using the comparative Ct method with glyceraldehyde-3-phosphate dehydrogenase (GAPD) as the internal control.
Effect of the 34 kb deletion on ANXA10 transcripts
The abomasum was collected from carrier and control steers at a slaughterhouse (Hyogo prefecture). A portion of ANXA10 cDNA, across exons 1 to 2 or exons 1 to 7, was amplified with primers as described in Additional file 1: Table S12. PCR products were electrophoresed by 1.5% agarose gel. The PCR product for exons 1 to exon 7 was subcloned into pCR2.1-TOPO vector (Invitrogen, Cat. #K4500-01) and sequenced using M13_R and M13 (-20) primer.
Western blotting
To express bovine and mouse ANXA10 proteins, the coding region of bovine ANXA10 (NM_001192058, 167–1141 bp) was amplified using PrimeSTAR Max DNA polymerase (Takara, Cat. #R045A) from cDNA of the abomasum using a forward primer (5′- GCTCTAGAatgttttgcggagactatgttcaggg -3′; uppercase letters indicate the XbaI linker) and a reverse primer (5′- CGGAATTCctaAGCGTAATCAGGAACGTCGTAAGGGTAgtagtcatctgcgtcacccgc -3′; uppercase letters indicate the EcoRI linker, and underlined letters indicate the C-terminal hemagglutinin [HA] tag for ANXA10, respectively). The PCR product was cloned into the XbaI and EcoRI sites of the pCAGGS vector [47]. The coding region of mouse Anxa10 (NM_001136089, 166–1,140 bp) was PCR amplified from cDNA of secretroy stomach using a forward primer (5′- GCTCTAGAatgttttgcggggaatatgtccaagg -3′; uppercase letters indicate the XhaI linker) and a reverse primer (5′- CGGAATTCctaAGCGTAATCAGGAACGTCGTAAGGGTAgtagtcttccacatcaccagc -3′; uppercase letters indicate the EcoRI linker, and underlined letters indicate the C-terminal hemagglutinin [HA] tag for ANXA10, respectively). The PCR product was cloned into the XhaI and EcoRI sites of the pCAGGS vector [47]. The sequence and orientation of the insert were confirmed by sequencing. The pCAGGS Vector was used for mock transfections. For cell culture, HeLa S3 cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM; Sigma, Cat. #D5796) with 10% fetal calf serum (FCS; Sigma, Cat. #F-2442) supplemented with 2 mM L-glutamine (Gibco, Cat. #25030-081) and 100 units/ml penicillin and 100 μg/ml streptomycin (Gibco, Cat. #15140-122). Using Lipofectamine 2000 (Invitrogen, Cat. #11668-019), we transfected 2 × 105 cells per well in a 6-well plate with a mixture of 2 μg of the vector. The transfected cells and tissues from cattle were homogenized with RIPA buffer (25 mM Tris–HCl, pH 7.6, 150 mM NaCl, 1% NP-40 [sigma, Cat. #21-3277-2], 1% sodium deoxycholate [sigma, Cat. #D5670], 0.1% sodium dodecyl sulfate [Wako, Cat. #191-07145], protease inhibitor cocktail [Roche, Cat. #11 873 580 001], 1 mM dithiothreitol [Wako, Cat. #191-07145]) and were sonicated with Bioruptor (Cosmo bio, Cat. #UCD-250) for four cycles at a high setting, 30 s ON and 30 s OFF on ice. The expression of bovine ANXA10 was confirmed by western blotting with anti-ANXA10 antibody to full length (R&D, Cat. #AF3544, 0.5 μg/ml) and to the C-terminal of ANXA10 (abcam, Cat. #ab135985, 1 μg/ml). The expression of mouse ANXA10 was confirmed by western blotting with anti-ANXA10 antibody (R&D, Cat. #AF3544, 0.5 μg/ml and Santa Cruz, Cat. #sc-135234, 0.2 μg/ml). Immunoreactivity was detected with a horseradish peroxidase-conjugated anti-goat IgG antibody (Jackson ImmunoResearch, Cat. #205-035-108) or anti-rabbit IgG antibody (Jackson ImmunoResearch, Cat. #711-036-152) and the ECL Prime Western Blotting Detection Reagent (GE Healthcare, Cat. #RPN2232). Chemiluminescence was detected with an ImageQuant LAS 4000 (GE Healthcare) and quantified using the ImageQuant TL Analysis Toolbox.
Replication study and frequency of the CNVR_322 in Japanese Black cattle population
For the replication study, we used 2693 cows, selected from 79,617 reproductive females that were evaluated for embryonic mortality. The 34 kb deletion was genotyped from PCR products using multiplex primers as described in Additional file 1: Table S10. The effects of the CNVR_322 were estimated as the least square means of generalized linear model (GLM) analyses. The statistical model for GLM analysis included fixed effects for the farm, birth year and CNVR_322. The genetic variance explained by the CNVR_322 was calculated based on estimates of the CNVR_322 effect and the frequency of the CNVR_322 [48]. Total genetic variance was estimated by MTDF-REML. The effect size of a CNVR_322 was estimated as the proportion of genetic variance explained by the CNVR_322.
To survey the risk-allele frequency in Japanese Black cattle, we genotyped 6364 cows that were not evaluated for embryonic mortality between 1990 and 2014 (mating records of the cows were not available) and 2391 steers from two central slaughterhouses (Tokyo Metropolitan Central Wholesale Market, Tokyo, and Nanko Wholesale Market, Osaka, Japan) that received animals from nationwide locations in Japan between 2003 and 2008.
Association analyses between CNVR_322 and meat traits in Japanese Black cattle
Data on five carcass traits were collected from 1156 animals, which were genotyped using 50 K SNP array [49]. The 34 kb deletion was genotyped from PCR products using multiplex primers as described in Additional file 1: Table S10. Association analyses were performed using Wald test with a linear mixed model accounting for family relatedness using a genomic relationship matrix.
Generation of Anxa10 null mice using CRISPR/Cas9
To design a single guide (sg) RNA for Anxa10, we searched PAM sequences using CRISPRdirect [50] and off-target sequences using COSMID [51]. We selected 69–71 bp of PAM sequences from the start of methionine in Anxa10 (NM_001136089, Additional file 4a) because the sgRNA-targeting site did not have off-target sequences on the genome. The forward (5′-CACCggagtgctccgcctagcatt-3′; uppercase indicates the BpiI linker) and reverse primers (5′-AAACaatgctaggcggagcactcc-3′; uppercase indicates the BpiI linker) were annealed and then cloned into BpiI (Thermo Scientific, Cat. #ER1012)-digested pX330-U6-Chimeric_BB-CBh-hSpCas9 plasmid (pX330-sgRNA) (Addgene, Cat. #42230) [52]. The sequence of the insert and the direction were confirmed by sequencing using primer (5′-tggactatcatatgcttacc-3′). To confirm the activity of targeted endonucleases, targeted regions of Anxa10 were cloned into pCAG-EGxxFP plasmids (Addgene, Cat. #50716) [37] to monitor the single strand annealing activity [37]. The targeted regions of Anxa10 were amplified using PrimeSTAR Max DNA polymerase (Takara, Cat. #R045A) from genomic DNA of C57BL/6NJ mouse using a forward primer (5′-CGGGATCCtgggggtaaaaagtggtgaa-3′; uppercase indicates the BamHI linker) and a reverse primer (5′-CGGAATTCcaatggtcagtgtgggtcag-3′; uppercase indicates the EcoRI linker). Using Lipofectamine 2000 (Invitrogen, Cat. #11668-019), we transfected 2 × 105 Cos-7 cells per well in a 6-well plate with a mixture of 0.5 μg of the pX330-sgRNA and 0.5 μg of pCAG-EGxxFP to SSA activity. After 24 h transfection, cells were examined with a confocal microscope (FV1000, Olympus Optical) and EGFP-positive cells were counted using ImageJ 1.46 [53].
To inject the plasmids pronuclear, BDF1 (CLEA-Japan) female mice were superovulated using 5U/body PMSG (ASKA) and 5U/body hCG (ASKA) and mated with BDF1 males. Fertilized eggs were collected from the oviduct and cultured in KSOM medium (Millipore, Cat. # MR-020P-5 F). The pronuclear stage eggs were injected with 5 ng/μl pX330-sgRNA (10 mM Tris–HCl, pH 7.5, 0.1 mM EDTA). The eggs were cultured in KSOM overnight then transferred into the oviducts of pseudo-pregnant ICR female mice (CLEA-Japan). The tail genomic DNA were amplified by TaKaRa Ex Taq HS DNA polymerase (TaKaRa, Cat. #RR006) using forward primer (5′-tttcatggaaattcatacaatctctt-3′) and reverse primer (5′-accttggagcattgagtcat-3′). The PCR conditions were as follows: initial denaturation was carried out at 98 °C for 5 min, and followed by 35 cycles, each consisting of denaturation at 98 °C for 20 s, annealing at 60 °C for 30 s and extension at 72 °C for 1 min. The PCR product was directly sequenced using reverse primer by BigDye Terminator v.3.1 Cycle Sequencing Kits (Applied Biosystems) and purified with CleanSEQ (Agencourt, Cat. #A29154) were sequenced, followed by electrophoresis using an ABI 3730 sequencer (Applied Biosystems).
Histological examination and immunohistochemistry
For histopathological examination, the tissue were fixed with 4% paraformaldehyde (PFA) and embedded in paraffin (JUNSEI, Cat. #58290-1501) using standard procedures. Sections (3 μm) were stained with Hematoxylin & Eosin.
For ANXA10 immunostaining, cattle abomasum samples were collected into cryotube and snap frozen with liquid nitrogen. Tissues were trimmed using cryostat blade (Leica, Cat. #14035838926) and then embedded in Tissue-Tek O.C.T. compound (Sakura, Cat. #4583) and frozen. Section (15 μm) was fixed in 4% PFA at 4 °C for 15 min. Mice stomach tissue were fixed with 4% PFA, and the tissues were rinsed with PBS, 30% sucrose/PBS and embedded in O.C.T compound. Sections (20 μm) were washed with 0.05% Triton x100 (Wako, Cat. #168-11805) in PBS. For ANXA10 immunostaining, the sections were immunostained with anti-Annxin A10 antibody (R&D, Cat. #AF3544, 5 μg/ml) in Canget signal immunostain solution A (TOYOBO, Cat. #NKB-401) for 1 h. For platelet thrombi, the sections were immunostained with anti-integrin beta 3 [EPR2417Y] (abcam, Cat. #ab75872, 1:250) [40].