Mitochondrial DNA T7719G in tRNA-Lys gene affects litter size in Small-tailed Han sheep

In farm animals, mitochondrial DNA (mtDNA) effect on economic performance remains hot-topic for breeding and genetic selection. Here, 53 maternal lineages of Small-tailed Han sheep were used to investigate the association of mitochondrial DNA variations and the lambing litter size. Sequence sweeping of the mitochondrial coding regions discovered 31 non-synonymous mutations, and the association study revealed that T7719G in mtDNA tRNA-Lys gene was associated with litter size (P < 0.05), manifesting 0.29 lambs per litter between the G and T carriers. Furthermore, using the mixed linear model, we assayed the potential association of the ovine litter size and haplogroups and multiple-level mtDNA haplotypes, including general haplotypes, assembled haplotypes of electron transport chain contained sequences (H-ETC), mitochondrial respiratory complex contained sequences (H-MRC) and mitochondrial genes (H-gene, including polypeptide-coding genes, rRNA genes and tRNA genes). The strategy for assembled mitochondrial haplotypes was proposed for the first time in mtDNA association analyses on economic traits, although none of the significant relations could be concluded (P > 0.05). In addition, the nuclear major gene BMPR1B was significantly correlated with litter size in the flock (P < 0.05), however, did not interact with mtDNA T7719G mutation (P > 0.05). Our results highlight mutations of ovine mitochondrial coding genes, suggesting T7719G in tRNA-Lys gene be a potentially useful marker for selection of sheep litter size.


Background
Mitochondria are responsible for ATP production in the electron transport chain (ETC) in cells. The ETC consists of five mitochondrial respiratory complexes (MRCs), of which complex I, complex III, complex IV and complex V are constructed by both mitochondrial and nuclear encoded proteins, and complex II is entirely encoded by nuclear genes. Mitochondrial genome codes 13 polypeptides, 2 rRNAs and 22 tRNAs [1]. Since the first report of mitochondrial DNA (mtDNA) effect on milk production traits in dairy cattle [2], mtDNA effects on various economic traits have been widely studied in livestock, including pigs [3,4], dairy cattle [5,6], beef cattle [7,8], sheep [9,10] and chickens [11,12]. Litter size is one of the vital economic traits for animal breeding and production, which has been studied for decades. For sheep, the nuclear gene, BMPR1B was identified as one of causative genes for sheep prolificacy [13], and has been widely used in sheep breeding. For the mtDNA effect on litter size, researchers reported the association with ewe litter size among haplogroups in an Afec-Assaf flock, but did not found the interaction with BMPR1B effect [10]. However, previous studies reported poor mtDNA effects [14][15][16][17], which made it necessary to uncover the genetic contribution of mtDNA for ewe litter size. In this study, Small-tailed Han sheep, a prolific breed of China, were used to investigate the association of litter size with mtDNA coding genes, in which nonsynonymous mutations were considered as possible functional SNPs. Besides the non-synonymous mutation and haplogroup, assembled haplotypes of ETC-contained mtDNA sequences, MRC-contained mtDNA sequences, and mitochondrial genes were analyzed, respectively. The strategy for assembled mitochondrial haplotypes was proposed for the first time in association analyses. In addition, BMPR1B effect and interaction with mtDNA mutations were analyzed.

Animals
In total, 117 lambing Small-tailed Han sheep from 53 maternal lineages (families divided by female ancestors) of the same flock were performed blood sample collection, and recorded one or more times of litter size (Additional file 1: Table S1). All sheep were kept indoors year-round, and fed with mixed silage and hay to meet their nutritional requirements. Litter size (number of lambs born) was recorded, and full pedigree information was collected for all animals as reproductive management of the flock where ewes were mated with the selected rams after spontaneous estrus. All the ewes were genotyped for the BMPR1B at >2-month age.
Genotyping of BMPR1B and sequencing of mitochondrial coding genes Genomic DNA was extracted using the standard phenol/ chloroform method [18]. The BMPR1B was genotyped using PCR-RFLP assay with the Ava II restriction enzyme [19]. Mitochondrial complete coding sequences were amplified by 17 primer pairs [20], and PCR products were sequenced in the Sanger method.

Haplotype and haplogroup constitution
All non-synonymous mutations were used to constitute the haplotype and haplogroup. Considering the possible action for mitochondrial function, haplotypes were furthermore assembled by mtDNA sequence of ETC, MRCs and genes respectively. Haplotypes were determined by the online software FaBox [21], and haplogroups were constituted based on network analysis by Network 4.6.1.4 [22].

Association analysis
Association analyses were carried out in the following mixed model by MIXED procedure in SAS software version 9.2 (SAS Institute Inc., Cary, North Carolina, USA).
In the model, the effects of lambing year-season (ys), parity number (parity), service ram (ram), BMPR1B genotype (BMPR1B), mtDNA mutations (mutations, including the effects of mtDNA non-synonymous mutations, haplotypes and haplogroups), the interaction between BMPR1B and mtDNA mutations (BMPR1B × mutations), the polygenic effect (ID), the permanent environmental effect (EP), and the random residual (e) were included. The response variable was the ewe litter size (ls). Each cell of these effects contained observations. The polygenic effect corrected the genetic background by the additive genetic relationship matrix, i.e. the pedigree information. The permanent environmental effect dealed with the repeated measurement data.  When a set of statistical inferences were simultaneously considered, multiple comparisons were conducted by the FDR using the R project. "ns" represents "not significant", and "*" represents "significant" at the significant level of 0.05 Inference about the interaction effect was made. If non-significant, that effect was dropped from the model and inference was made about the main effects. If significant, the cell means for the interaction became of interest.
When a set of statistical inferences were simultaneously considered, multiple comparisons were conducted by the false discovery rate (FDR) in the R project (R version 3.2.5) [23].

Mutations in mitochondrial coding genes
In total, 95 mutations in mitochondrial coding genes were discovered (Additional file 1: Table S2), including 64 synonymous and 31 non-synonymous mutations (19 missense mutations in protein coding genes, 8 mutations in rRNAs, and 4 mutations in tRNAs respectively), which were illustrated in Table 1.

Effects of mitochondrial haplotype and haplogroup on ovine litter size
The 31 mitochondrial non-synonymous mutations assigned the Small-tailed Han sheep flock to 44 haplotypes (Additional file 1: Table S3), which were clustered into 4 haplogroups (Fig. 1). Using the mixed linear model, no significant effect on any haplotype or haplogroup was associated with litter size of Small-tailed Han sheep (P > 0.05), and the interaction between BMPR1B and haplotype or haplogroup also did not significantly affect litter size (P > 0.05), while the BMPR1B was positively associated with litter size (P < 0.05).

Effects of missense mutations and haplotypes in mitochondrial protein coding genes on ovine litter size
Totally 35 H-ETCs were constituted by 19 missense mutations in protein coding sequences, while 15 missense mutations in MRCI assembled 33 H-MRCIs, and both 2 mutations in MRCIV and MRCV constituted 3 H-MRCIVs and 3 H-MRCVs respectively. Intensively, H-genes were assembled by mutations of individual genes ( Table 2 and Additional file 1: Table S4, S5 and S6). With mixed linear model analyses, no significant association was detected between litter size and the 19 non-synonymous mutations (P > 0.05) ( Table 1), nor was any haplotype (P > 0.05) ( Table 2). BMPR1B was strongly associated with litter size (P < 0.05) ( Table 3), but the interaction of BMPR1B and mtDNA mutations was not remarkable (P > 0.05).

Effects of mutations and haplotypes in mitochondrial rRNA genes on ovine litter size
The 3 mutations in 12SrRNA sorted the sheep flock into 4 haplotypes (H-12SrRNA), and the 5 mutations in 16SrRNA constituted 15 haplotypes (H-16SrRNA) ( Table 2 and Additional file 1: Table S7). Association analyses revealed that neither mutation nor haplotype in rRNA genes affected litter size (P > 0.05) (Tables 1 and 2). The BMPR1B genotype was remarkably associated with litter size (P < 0.05) (Table 3), however, was inconspicuously correlated to rRNA mutations (P > 0.05).

Effects of mutations in mitochondrial tRNA genes on ovine litter size
There were 4 mutations in tRNA genes, and only one variation was observed in each of them (Table 1), Fig. 1 Neighbor-joining tree of different maternal pedigrees for haplogroup assignments by the 31 non-synonymous mutations therefore no haplotype was constituted. Notably, T7719G in tRNA-Lys affected litter size by 0.29 lambs per litter between the G (1.79) and T (1.50) carriers (P < 0.05). Even though the BMPR1B was in prominent correlation with litter size (P < 0.05) ( Table 3), there was no interaction with T7719G (P > 0.05) (Additional file 1: Table S8).

Discussion
In human diseases, mutations in mitochondrial coding genes led to changes of oxidative phosphorylation enzyme complexes [24]. It is rational to image the possible effects of mtDNA on livestock traits. In this study, the Small-tailed Han sheep were used to explore the correlation of mtDNA and ewe litter size.
In order to investigate the overall mtDNA effects, a novel strategy of assembling mitochondrial haplotypes based on biology functions is proposed. Extensively, there are four levels of mitochondrial haplotypes for mitochondrial biological functions. The general haplotype is assembled by all non-synonymous mutations to reflect the integrated characteristics of mtDNA coding In the mixed linear model, unnecessary factors interfere with the estimation of interested factors, as the degree of freedom is wasted. Therefore, based on the thoughts of multiple models [10] and variable selection in a mixed linear model [25], a two-step method of environmental variable selection was put forward in the association analyses, in which nonsignificant environmental factors were ignored to improve accuracy of genetic estimation. In this study, the environmental factors included lambing month, parity of dam, and service ram, and the genetic factors, which we were interested in, included BMPR1B genotype, mtDNA mutations and their interaction. We started a full model (Model I) with all the environmental and genetic factors to test if environmental factors were noises. Results revealed that the effects of parity of dam and service ram were not significant on litter size. Subsequently, these two factors were excluded in the aims to construct an optimized model (Model II), which was used to test the significance of genetic factors, especially the effects of mtDNA mutations on litter size.
Many mitochondrial tRNA mutations in human were reported to be associated with wide range of pathological  When a set of statistical inferences were simultaneously considered, multiple comparisons were conducted by the FDR using the R project. "ns" represents "not significant", and "*" represents "significant" at the significant level of 0.05 conditions [26], for example, deafness was correlated with mitochondrial tRNA-Asp A7551G mutation [27], hypertension was associated with mitochondrial tRNA-Ile A4263G mutation [28], and mitochondrial tRNA mutations might decrease in carcinoma hepatocyte [29]. In sheep, A7755G in tRNA-Lys was reported to affect litter size in an Afec-Assaf flock [10], while our study revealed that T7719G in tRNA-Lys was significantly associated with ewe litter size in Small-tailed Han sheep, with the means of 1.79 and 1.50 for the G and T carriers, manifesting a difference of 0.29 lambs per litter. Furthermore, the mutation was predicted to produce a transversion of U to G at DHU loop on 2D cloverleaf of the tRNA-Lys structure (Fig. 2).

Conclusions
As a conclusion, we summarize the research procedures for mtDNA effects. Firstly, sweeping mitochondrial variations on maternal lineages; secondly, constituting the haplotype, haplogroup and assembled haplotypes; lastly, analyzing the association between mtDNA mutations (individual mutations, haplotypes and haplogroups) and interested traits. The present study discovered the mtDNA T7719G was linked with ewe litter size. For traits of low heritability, the marker assisted selection could increase the accuracy of breeding and selection. For the further study, the mtDNA T7719G should be put into post-association validation, and may become a genetic marker in the sheep breeding programs.

Additional file
Additional file 1: Table S1. Detail information of used sheep. Table S2.
Information of variations in mtDNA coding region of Small-tailed Han sheep. Table S3. The haplotype constituted by the 31 mitochondrial nonsynonymous mutations in mitochondrial coding regions (H). Table S4.