NRF1 and ZSCAN10 bind to the promoter region of the SIX1 gene and their effects body measurements in Qinchuan cattle

The SIX1 homeobox gene belongs to the six homeodomain family and is widely thought to play a principal role in mediating of skeletal muscle development. In the present study, we determined that the bovine SIX1 gene was highly expressed in the longissimus thoracis and physiologically immature individuals. DNA sequencing of 428 individual Qinchuan cattle identified nine single nucleotide polymorphisms (SNPs) in the promoter region of the SIX1 gene. Using a series of 5′ deletion promoter plasmid luciferase reporter assays and 5′-rapid amplification of cDNA end analysis (RACE), two of these SNPs were found to be located in the proximal minimal promoter region −216/−28 relative to the transcriptional start site (TSS). Correlation analysis showed the combined haplotypes H1-H2 (-GG-GA-) was significantly greater in the body measurement traits (BMTs) than the others, which was consistent with the results showing that the transcriptional activity of Hap2 was higher than the others in Qinchuan cattle myoblast cells. Furthermore, the electrophoretic mobility shift assays (EMSA) and chromatin immunoprecipitation assay (ChIP) demonstrated that NRF1 and ZSCAN10 binding occurred in the promoter region of diplotypes H1-H2 to regulate SIX1 transcriptional activity. This information may be useful for molecular marker-assisted selection (MAS) in cattle breeding.


Temporal expression of SIX1 in longissimus thoracis during different developmental stages.
The expression of the bovine SIX1 gene in longissimus thoracis at 1, 3, 6 and 12 months of age was significantly higher than that during any other stages ( Fig. 1b and c). However, the expression level was reduced dramatically after 18 months (Fig. 1b and c), but no significant difference was observed in the expression level after 24 months, and the expression level trend remained stable at both the mRNA and protein levels. Based on these results, there appears to be a close relationship between the growth rate traits and SIX1 gene expression in Qinchuan cattle.
Determination of the transcription start site of the SIX1 gene. To analyze the structure and molecular mechanisms of the bovine SIX1 gene, we performed 5′-RACE to identify the transcriptional start site (TSS) of SIX1. Two successive rounds of PCR were performed using antisense primer R1 and nested primer R2 (Table S1). As shown in Fig. 2a, two bands of 555 and 475 bp were amplified. In total, 18 positive clones had three different 5′ ends in the first exons, i.e., at 316, 297 and 236 bp upstream of the TSS (Fig. 2b). Sequence alignment showed that the most upstream of the TSS was completely in accord with the published SIX1 mRNA sequence (XM_588692.7). Therefore, we verified the site and designated as +1.
Seven SNPs were identified in the 1. 8   The samples of the longissimus thoracis were obtained at 1, 3, 6, 12, 18, 24, 36 and 48 months after birth. SIX1 mRNA expression was normalized to the housekeeping gene GAPDH and the expression levels were calculated relative to the gene expression in the liver and 24 months, respectively. The value of each column represents the mean ± standard deviation of three independent experiments. The unpaired Student's t-test was used to detect significant differences. "*"P < 0.05 and "**"P < 0.01.
SCientifiC RepoRts | 7: 7867 | DOI: 10.1038/s41598-017-08384-1 Isolation of the functional proximal minimal promoter of the SIX1 gene. The proximal minimal promoter of the SIX1 gene was isolated to determine the SNPs that are responsible for influencing the function of promoter. Seven reporter constructs with progressive deletions from the 5′ end of the promoter were generated. The luciferase reporter constructs, named pGL3-1802, pGL3-1488, pGL3-1044, pGL3-708, pGL3-483, pGL3-216 and pGL3-28, were transfected into Qinchuan cattle myoblast cells (QCMCs). The luciferase assays revealed a 5-7-fold increase in the promoter activity of pGL-1802/+142 compared to that in the empty vector, indicating a functional promoter in the −1802/+142 region of the SIX1 gene. When the promoter was deleted to position −1044, the promoter activity of pGL-1044/+142 was decreased by 32% compared with the pGL-1488/+142 (Fig. 3a). After inducing the QCMCs via HS after the transfection, the results of the luciferase assays were similar  Luciferase activities in the bovine SIX1 promoter constructs in QCMCs. (a) All plasmids containing 5′ unidirectional deletions of the promoter region of the SIX1 gene (pGL3 −1802, −1488, −1044, −708, −483, −216, −28 and pGL3-Basic) were transfected into QCMCs. After 5 h we replaced the transfection mixture with DMEM with 2% HS (myotubes). The cells were collected for the luciferase assay at 48 h. The results are expressed as the mean ± standard deviation in arbitrary units based on firefly luciferase activity normalized to the Renilla luciferase activity in triplicate transfections. The unpaired Student's t-test was used to detect significant differences. "*"P < 0.05 and "**"P < 0.01. (b) A graphical representation of the bovine SIX1 gene proximal promoter region from +1 to −1802 base pairs, predicting the regions with a high GC content. Folded lines indicate the GC percentage, which is represented on the y-axis, and the x-axis denotes the bp position on the 5′ untranslated region; the bottom of the blue area indicates the relative positions of the CpG islands. Coordinates are given relative to the translational start site (shown as +1). Arrows indicate TSS-1, positions −216 and −28 bp and nine SNP loci in the promoter, of which SNP1 and SNP2 were located in the proximal minimal promoter of the SIX1 gene.
to those during the undifferentiated stage (Fig. 3a). This result demonstrated that positive regulatory elements were located in the −1488/−1044 region. However, when the promoter was further deleted to −28 bp, the promoter activity of pGL-28/+142 was almost abolished in both the undifferentiated and differentiated QCMCs compared to that of pGL-216/+142. These results indicate that the proximal minimal promoter of the SIX1 gene is located within the region −216/−28 relative to TSS-1 and that the undifferentiated cell model is better for determining transcriptional activity.
Genetic polymorphism of in the Qinchuan Cattle SIX1 proximal minimal promoter region. In the present study, the genetic parameters of the bovine SIX1 gene, including the genotype and allele frequencies, were directly calculated for all 428 animals. SNP1 (g. −85 G > T) and SNP2 (g. −63 G > A) were located in the proximal minimal promoter region −216/−28 of the SIX1 gene. We also determined the genotypic and allelic frequencies, genetic diversity parameters of heterozygosity (He), effective allele numbers (Ne), polymorphism information content (PIC) and Hardy-Weinberg equilibrium ( Table 1). The results indicated that for the g. −85 G > T mutation, the TT genotype (14.96%) was less frequent than the wild allele GG (44.16%). The allele frequencies, He, Ne and PIC at the current locus were 0.6460 (G), 0.3540 (T), 0.4574, 1.8428 and 0.3528, respectively. For the g. −63 G > A mutation, the GG genotype was the most prevalent (62.15%) followed by GG (30.84%) and AA (7.01%). The values of the allele frequencies, He, Ne, and PIC at this locus were 0.7757 (G), 0.2243 (A), 0.3480, 1.5537 and 0.2874, respectively. In our study, g. −85 G > T and g. −65 G > A had an intermediate level of polymorphism (high level, PIC > 0.5; moderate level, 0.5 < PIC > 0.25; low level, PIC < 0.25) 17 . The genotypic distributions of SNP1 and SNP2 were in Hardy-Weinberg disequilibrium (chi-square test, χ 2 < χ 2 0.05 ).

Effects of single marker on BMTs.
The association between the two SNPs (g. −85 G > T and g. −63 G > A) and the seven economic traits of body length (BL), withers height (WH), chest depth (CD), chest circumference (CC), back fat thickness (BF), ultrasound loin muscle area (ULA) and intramuscular fat content (IFC) are presented in Table 2. Individuals with genotype the GG had a significantly greater CC and IFC than those with the TT genotype (P < 0.05) for g. −85 G > T. For g. −63 G > A, individuals with the AA genotype had a significantly greater BL, CD, CC, BF and ULA than those with the GA and GG genotypes (P < 0.05) ( Table 2). These results indicated that the two loci g. −85 G > T and g. −63 G > A were associated with the BMTs in Qinchuan cattle.

Effects of diplotypes on BMTs.
To further analyse the associations between the diplotypes of the two the SNPs and the BMTs, six diplotypes were combined by four haplotypes (Hap1, Hap2, Hap3 and Hap4) (  (Table 4). Additionally, the H 3 -H 4 diplotypes had a significantly higher BL, CD, CC and BF (P < 0.05) than the H 1 -H 1 , H 1 -H 3 and H 3 -H 3 diplotypes (Table 4). In contrast, the H 3 -H 4 diplotypes displayed a reduced IFC (P < 0.05).

Potential transcription-factors in the SIX1 haplotypes possessed different transcriptional activities.
To detect the transcriptional activities of the haplotypes, the various haplotypes were cloned, and then luciferase reporter, named pGL3-Hap1, pGL3-Hap2, pGL3-Hap3 and pGL3-Hap4, were constructed. The transcriptional activities of these haplotypes were determined using a dual-luciferase reporter assay system in the QCMCs. The results showed that Hap1 had a 0.83-fold (P < 0.01), 0.04-fold (P > 0.05) and 0.91-fold (P < 0.01) lower activity than Hap2, Hap3 and Hap4, respectively (Fig. 4). A pairwise comparison of the transcriptional activities of Hap1 and Hap2 (different at the SNP2 locus, P Hap1/2 < 0.01), Hap3 and Hap4 (different at the SNP2 locus, P Hap3/4 < 0.01), Hap1 and Hap3 (different at the SNP1 locus, P Hap1/3 > 0.05), Hap2 and Hap4 (different at the SNP1 locus, P Hap2/4 > 0.05) was performed. We hypothesized that SNP2 in the g. −63 G > A mutation potentially resulted in the binding of the transcription-factor Nuclear respiratory factor 1 (NRF1) and Zinc finger and SCAN domain containing 10 (ZSCAN10) proteins, which may play a prominent role in the regulation of transcription activities (Table S2, Fig. 5a).

NRF1 and ZSCAN10 bind to the promoter region of the bovine SIX1 gene in vitro and in vivo.
Electrophoretic mobility shift assays (EMSAs) and a chromatin immunoprecipitation assay (ChIP) were used to determine whether NRF1 and ZSCAN10 bind to the promoter region of the bovine SIX1 gene (Fig. 5a), resulting in increased transcriptional activities of Hap2 and Hap4. As shown in Fig. 5b, the nuclear protein from the QCMCs bound to the 5′-biotin labelled NRF1 probes and formed two main complexes (lane 2, Fig. 5b). The competition assays verified that the mutant probe had a slight effect on the main complexes (lane 3 Fig. 5b). However,  Table 1. Genotype frequencies (%) of the SIX1 gene for the SNPs in the Qinchuan cattle populations. HWE, Hardy-Weinberg equilibrium; χ0.05 2 = 5.991, χ0.01 2 = 9.210. specificity of the NRF1/DNA interaction was prevented by competition from excess non-labelled DNA (lane 4, Fig. 5b). The final lane shows that the complex was super-shifted when it was incubated with the NRF1-antibody (lane 5, Fig. 5b). ZSCAN10 yielded similar results as NRF1 (Fig. 5c). Although the EMSA in the ZSCAN10 experiments did not reveal a super-shifted product at the ZSCAN10 binding sites, however, the brand of the main complex was clearly decreased. The super-shifted may have formed a high molecular weight polymer, which caused a reduced gel mobility shift (lane 5, Fig. 5c). The ChIP results revealed that NRF1 and ZSCAN10 interacted with the binding sites; the relative enrichment levels were ~4.9 and ~6.6-folds over the IgG control respectively ( Fig. 6a  and b), respectively, based on three independent experiments.

Discussion
Identifying QTLs and candidate genes that can be utilized in marker-assisted breeding through the manifestation of economically important traits will facilitate Qinchuan cattle breeding programmes. Variants of candidate genes can be associated with economically important traits, such as growth and carcass traits 18,19 . SIX1 has emerged as a candidate gene that is an important regulator of vertebrate development and the maintenance of differentiated tissues states 20,21 . Our work supports that SIX1 is one of the genes that controls myogenesis and may influence BMTs.
In the present study, the tissue distribution of SIX1 mRNA showed the highest expression in the longissimus thoracis (Fig. 1a), which is consistent with previous findings observed in other species, such as in human 22 , porcine 23 and duck 24 . However, the SIX1 expression level is very low in the kidney in this study, which differs from the results of previous studies in which SIX1 was shown to be a crucial factor for early kidney development 6,7 . The main reason for this discrepancy is that SIX1 expression is regulated temporally in the kidney during early embryonic development 6,7 . Previous studies have mainly studied the regulatory mechanism of the SIX1 gene in early embryogenesis, particularly in the successive steps of myogenesis, but less is known about its expression pattern and regulatory mechanism during later development. In the present study, we detected an up-regulation of SIX1 expression from 1 to 12 months after birth; however, the expression of SIX1 decreased dramatically after 18 months. Increases in skeletal muscle mass are largely determined by hypertrophy of muscle fibers during postnatal growth because the number and size of muscle fibres do not change after birth 25 . An increase in muscle mass solely through muscle fibres hypertrophy could influence the body measurement and meat quality 26 . Qinchuan cattle require on average of 12 to 24 months attaining physiological and skeletal maturity, which appears closely related to the decreasing SIX1 expression during this time ( Fig. 1b and c). Overall, the mechanism of the bovine SIX1 gene may be associated not only with the early events of myogenesis but also with muscle growth and meat quality.
Previous studies have shown that genetic polymorphisms in SIX1 are associated with BMTs. In pigs, Wu et al. 23 reported that a C/T and A/G polymorphism in the SIX1 promoter and intron were significantly associated with the meat color value (MCV1) and meat marbling (MM1) of longissimus dorsi and dressing percentage (DP). In the present study, the correlation analysis showed that cattle with the G allele at g. −85 G > T and the A allele at g. −63 G > A loci had better BMTs. A pairwise comparison revealed that a SNP2 mutation at g. −63 G > A associated with the prominent role of transcription-factors NRF1 and ZSCAN10 in the proximal minimal promoter  region, which play a major role in transcriptional regulation 27,28 , thus contributing to regulate the transcriptional activity of the SIX1 gene. NRF1, which encodes a protein that homodimerizes and functions as a transcription factor, is associated with neurite outgrowth through the regulation of globin gene expression 29 . Previous studies have shown that the activation of NRF1 in fibroblasts induces increases in cytochrome c expression and mitochondrial respiratory capacity. Overexpressing NRF1 in skeletal muscle resulted in an increased expression of myocyte enhancer factor (MEF) 2 A and GLUT4 30 , which were associated with a proportional increase in insulin-stimulated glucose transport in muscles 31 . Moreover, NRF1 can interact with DYNLL1 32 , PPARGC1A 33 and Cytochrome c oxidase (COX) 34 , to initiate nitric oxide synthase 33 and active energy metabolism 34 for individual development. ZSCAN10 (also named as Zfp206) belongs to the C2H2 zinc fingers family and encodes a zinc finger transcription factor that is specifically expressed in embryonic stem cells (ESCs) 35,36 . ZSCAN10 is a regulator of pluripotency and maintains the pluripotent state by interacting with other pluripotency factors, such as Sox2 and Oct4 37 . ZSCAN10-null mice have a reduced weight and mild hypoplasia in the heart, spleen, eyes and long bones 38 . Consistently, ZSCAN10 expression in ESCs mid-gestation embryos and adults suggests that ZSCAN10 plays a role in maintaining progenitor cell subpopulation and regulating individual development. In the present study, we identified that NRF1 and ZSCAN10 binds to Hap2 and Hap4 and induce a significantly higher transcriptional activity than others based on the luciferase reporter, EMSA and ChIP assays. These results demonstrate that NRF1 and ZSCAN10 play important roles in regulating the transcriptional activity of the SIX1 gene and contributing to better body measurements.
In combination, the H 1 -H 2 and H 3 -H 4 diplotypes showed better BMTs than the H 1 -H 1 , H 1 -H 3 and H 3 -H 3 diplotypes. Correlation analysis results were consistent with the interpretation that the transcriptional activities of Hap2 had a higher activity than those Hap1 and Hap3. Notably, the H 3 -H 4 diplotypes showed a lower IFC than the other diplotypes, which is an indicator of lean meat quality characteristics, affects juiciness, flavour, and morphology and is closely related to beef tenderness 15 . There is significant breed variation in IFC, with the Belgian Blue bulls producing the leanest meat, Limousine producing intermediate levels and Aberdeen Angus producing the highest fat content 39,40 . This increased IFC has a positive impact on tenderness, taste and flavour 40 .   The unpaired Student's t-test was used to detect significant differences. "*"P < 0.05 and "**"P < 0.01.
Through aggregate selection, based on our findings, we infer that the H 1 -H 2 diplotypes could be used as a molecular marker of combined genotypes for future selection of BMTs in Qinchuan cattle.
In conclusion, our study revealed that bovine SIX1 is highly expressed in the longissimus thoracis and physiological immature individuals. Two SNPs were located within the proximal minimal promoter region from −216 to −28, inducing the binding of transcription factors NRF1 and ZSCAN10 on the promoter region of diplotypes H 1 -H 2 and regulating SIX1 transcriptional activity. Correlation analysis showed the combined haplotypes H 1 -H 2 (-GG-GA-) had significantly greater BMTs than the others. This information may be useful for MAS in cattle breeding.

Materials and Methods
Ethics Statement. All animal procedures were performed according to guidelines laid down by the China Council on Animal Care, and the protocols were approved by the Experimental Animal Manage Committee (EAMC) of Northwest A&F University.
Animal sources, data collection, and genomic DNA isolation. In total, 428 female cows aged 18 to 24 months were randomly selected from the National Beef Cattle Improvement Center's experimental farm (Yangling, China). The BMTs were measured as described previously 41 , including BL, WH, CD and CC. ULA, BF and IFC were measured via ultrasound using a Sono-grader model 2 (Renco, USA). Genomic DNA was extracted from blood samples, using standard method 42 , and stored at −20 °C until subsequent analyses.
Real-time PCR analysis of expression pattern. Fifteen tissues (heart, liver, spleen, kidney, rumen, reticulum, omasum, abomasa, small intestine, large intestine, subcutaneous fat, longissimus thoracis, soleus, psoas and testicle) were obtained from three adult Qinchuan cattle. The longissimus thoracis samples were collected during eight developmental stages from male Qinchuan cattle, including 1, 3, 6, 12, 18, 24, 36 and 48 months after birth, and three parallel individuals were sampled during each period. Total RNA was extracted from the tissues using a Total RNA Kit (Tiangen, Beijing, China) and then reverse-transcribed using a PrimeScript ™ RT Reagent Kit (Perfect Real Time) (TaKaRa, Dalian, China). The reaction was performed using a SYBR Green PCR Master Mix  Supplementary Table S1. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was invoked as the endogenous control gene. The relative expression levels of the target mRNAs were calculated using the 2 −ΔΔCt method 43 .
Western blotting. Tissues protein were extracted using the T-PER Tissue Protein Extraction Reagent (Pierce, Thermo Fisher Scientific, USA). The total protein samples were quantified using the Pierce BCA Protein Assay Kit (Thermo Scientific), and 50 μg were electrophoresed on a 10% SDS-polyacrylamide gel, and then transferred to nitrocellulose. After blocking in defatted milk powder, the membranes were incubated with the SIX1 antibody (sc-514441, Santa Cruz, USA) and GAPDH antibody (sc-293335, Santa Cruz). The blots were washed and subsequently treated with a peroxidase labelled secondary antibody. The signals were detected as chemical luminescence to expose X-ray films using the ChemiDoc ™ XRS+ System (Bio-Rad, Hercules, CA, USA).

5′-Rapid amplification of cDNA ends (5′-RACE).
To identify the transcriptional start site (TSS) of the bovine SIX1 gene, 5'-RACE was performed on total RNA from the longissimus thoracis muscle using a BD SMARTTM RACE cDNA amplification kit (Clontech Inc., CA, USA) according to the manufacturer's protocol. PCR was performed using a Universal Primer A Mix (UPM, Clontech Inc., CA, USA) and the nested PCR primers (Table S1) located in exon 1 of the SIX1 gene. The conditions and methods used were as previously described 27 . For sequencing PCR products were separated by electrophoresis in 2% agarose gels and subsequently cloned into T-Vector pMD19 (simple) (Takara).

Primer design, PCR amplification, and SNPs detection. Three pairs of PCR primers (primer A, B
and C) were designed to amplify a 1.8 kb genomic region upstream of the bovine SIX1 gene (AC_000167.1 from 73068130 to 73074697). The primer sequences are reported in Table S1. The PCR amplifications were carried out using pooled genomic DNA from 428 Qinchuan cattle as a template 44 Table S1 for 30 s and 72 °C for 60 s. All PCR products were sequenced to verify amplification of the intended target. Finally, the sequences were imported into BioXM software (Version 2.6) for the SNP analysis.
Promoter cloning and generation of the luciferase reporter constructs. Fragment primers for −1802, −1488, −1044, −708, −483, −216, −28 and +144 were designed (Table S1) to amplify unidirectional deletions of the bovine SIX1 promoter. Promoter constructs were generated by PCR using specific primers with the sequence of the KpnI and BglII restriction sites incorporated and the wild individuals as DNA templates. Then, all fragments were cloned into Vector pMD19-T (simple) (TaKaRa), and ligated into the luciferase reporter construct pGL3-basic vector digested with the same restriction enzymes KpnI and BglII (TaKaRa). These plasmids were named pGL3-1802, pGL3-1488, pGL3-1044, pGL3-708, pGL3-483, pGL3-216 and pGL3-28. The PCR Figure 6. ChIP assay of NRF1 and ZSCAN10 binding to the SIX1 promoter in vivo. We analysed the immunoprecipitated products of the NRF1 (a) and ZSCAN10 (b) antibodies via RT-PCR and ChIP-QPCR. We used total chromatin from muscles as the input, and normal rabbit IgG served as the negative control antibody. "**"P < 0.01. Error bars represent the SD (n = 3).
SCientifiC RepoRts | 7: 7867 | DOI:10.1038/s41598-017-08384-1 amplification conditions were similar to the conditions used in the previous step and plasmid DNA was further confirmed by DNA sequencing.
Potential cis-acting elements identification. The Genomatix database (http://www.genomatix.de) was used to search for potential cis-acting elements in the proximal minimal promoter region of the SIX1 gene. These potential cis-acting elements were compared for genotypic differences. Sequences that contained one or two SNPs or being adjacent to SNPs were evaluated.
Cell culture and transfection. QCMCs were isolated from Qinchuan foetal bovine as described previously 45,46 . The QCMCs were maintained in Dulbecco's Modified Eagle Medium (DMEM-F12) and supplemented with 20% newborn calf serum (NBCS, Invitrogen, USA) and antibiotics (100 IU/mL penicillin; 100 µg/ mL streptomycin) at 37 °C and 5% CO 2 in a normal atmosphere incubator. Cells were grown overnight to 80-90% confluence at a density of 1.2 × 10 5 cells in the growth medium without antibiotics in 24-well plates. In each well, 800 ng of the each series of reporter plasmid (pGL3 −1802, −1488, −1044, −708, −483, −216 and −28) were co-transfected with pRL-TK normalizing reporter plasmid (Promega, USA) into QCMCs with 3 μL X-tremeGENE HP DNA transfection reagent (Roche, USA). The pGL3-Basic vector served as a negative control. At 5 h after the transfection, we replaced the media with DMEM with 2% horse serum (HS) (GIBCO, Invitrogen) and incubated for 40 h to induce the differentiation of the QCMC myoblasts into myotubes. We performed all remaining steps as previously described 28 . Firefly luciferase activity and Renilla luciferase activity were measured according to the dual-luciferase reporter assay standard protocol in three independent experiments. The relative luciferase activities were determined using a NanoQuant Plate ™ (TECAN, infinite M200PRO).

Statistical analysis.
Statistical analysis was performed for allelic frequencies, genotype frequencies, He, Hardy-Weinberg equilibriums and PIC parameters according to Nei's methods 17 . Linkage Disequilibrium (LD) and haplotype distributions of the SNPs were analyzed using the expectation maximization algorithm with Haploview software 47 . The association between single SNPs and body measurements traits were analyzed using the general linear models (GLM) procedure in SPSS (version 13.0). The linear model was used: where Y ijkl were the traits measured on each individual cow, μ was the overall population mean for the traits, G i was the fixed genotype effect, A j was the fixed effect of age, A k was the fixed effect due to the age of dam, S l was the fixed effect due to the season of sampling (spring vs. fall) and E ijkl was the standard error.  Table S1) were incubated room temperature for 20 min with a reaction mixture containing 2 μL 10 × binding buffer, 1 μL poly (dI . dC) and 10 μg of nuclear protein extract in a volume of 20 μL. For the competition assay, unlabelled or mutated DNA probes were added to the reaction mixture and incubated for 15 min; then, 200 fmol labeled probes in a volume of 20 μL were added and incubated at room temperature for 20 min. For the super-shift assay, 10 μg of the NRF1 (ab86516, Abcam, USA) or ZSCAN10 (ab45344, Abcam) antibodies were added to the reaction mixture and then incubated on ice for 30 min; then, 200 fmol of the labelled probes in a volume of 20 μL were added and incubated at room temperature for 20 min. Finally, the main complexes were resolved on 6% non-denaturing polyacrylamide gel electrophoresis (PAGE) using 0.5 × TBE buffer for 1 h and image with ChemiDoc ™ XRS+ System molecular imager (Bio-Rad). ChIP assay. The ChIP assays were performed using the SimpleChIP ® Enzymatic Chromatin IP Kit (CST, Massachusetts, USA) according to the manufacturer's protocol. The samples (n = 3) from the Hap2 and Hap4 Qinchuan bovine were used. The protein-DNA complexes were cross-linked with 37% formaldehyde and neutralized with glycine. After digesting the DNA with micrococcal nuclease into fragments of approximately 150-900 bp in length, the fragmented chromatin samples were suspended in the ChIP dilution buffer. The cross-linked chromatin samples were immunoprecipitated with 4 μg of the NRF1 or ZSCAN10 antibodies and normal rabbit IgG overnight at 4 °C. The immunoprecipitated products were isolated with protein G agarose beads, and the bound chromatin was then collected with salt washes. The eluted ChIP in the Elution Buffer was then digested with proteinase K and purified for a quantitative PCR analysis. The ChIP primers used in the RT-PCR experiment are listed in Table S1. Percent input was calculated as follows: % Input = 2 [−ΔCt(Ct[ChIP] − (Ct[Input]−Log2(Input Dilution Factor)))] 48 . We used the immunoprecipitated products from the normal rabbit IgG group as a negative control.