Mitochondrial Genome Variation after Hybridization and Differences in the First and Second Generation Hybrids of Bream Fishes

Hybridization plays an important role in fish breeding. Bream fishes contribute a lot to aquaculture in China due to their economically valuable characteristics and the present study included five bream species, Megalobrama amblycephala, Megalobrama skolkovii, Megalobrama pellegrini, Megalobrama terminalis and Parabramis pekinensis. As maternal inheritance of mitochondrial genome (mitogenome) involves species specific regulation, we aimed to investigate in which way the inheritance of mitogenome is affected by hybridization in these fish species. With complete mitogenomes of 7 hybrid groups of bream species being firstly reported in the present study, a comparative analysis of 17 mitogenomes was conducted, including representatives of these 5 bream species, 6 first generation hybrids and 6 second generation hybrids. The results showed that these 17 mitogenomes shared the same gene arrangement, and had similar gene size and base composition. According to the phylogenetic analyses, all mitogenomes of the hybrids were consistent with a maternal inheritance. However, a certain number of variable sites were detected in all F1 hybrid groups compared to their female parents, especially in the group of M. terminalis (♀) × M. amblycephala (♂) (MT×MA), with a total of 86 variable sites between MT×MA and its female parent. Among the mitogenomes genes, the protein-coding gene nd5 displayed the highest variability. The number of variation sites was found to be related to phylogenetic relationship of the parents: the closer they are, the lower amount of variation sites their hybrids have. The second generation hybrids showed less mitogenome variation than that of first generation hybrids. The non-synonymous and synonymous substitution rates (dN/dS) were calculated between all the hybrids with their own female parents and the results indicated that most PCGs were under negative selection.


Introduction
Mitochondria not only provide energy for animal cells [1], but also play a crucial role in the dynamics of molecular evolution, species hybridization, gene introgression and species differentiation [2]. Mitochondrial genome (mitogenome) has been extensively used to evaluate genetic diversity of different populations in fish species and analyze genetic characteristics of the hybrid fishes due to its small molecular weight, maternal inheritance, relatively rapid substitution rate, and lack of recombination [2]. Guo et al. [3] proved that the atpase8 and atpase6 genes are useful genetic markers to monitor the variations in the hybrid progeny among the different carp strains. Avise and Saunders [4] took advantage of polymorphisms in mitochondrial DNA (mtDNA) to analyze hybridization and introgression among sunfish species (Lepomis, Centrarchidae).
Although mitochondrial DNA is usually maternally inherited, some studies also reported the variation of mitogenome sequences between hybrids and their female parents. In mammals, the studies in murine hybrids and the interspecies cross between the domestic cow and the Asian wild gaur proved the evidences of paternal inheritance and recombination of mtDNA [5][6][7]. These results indicated that maternal inheritance of mitochondrial genome may involve species specific regulation and can be disrupted by hybridization. Although many studies demonstrated that mtDNA of hybrid fishes followed strict rules of maternal inheritance [3,[8][9][10][11][12][13], variations of mitochondrial DNA between hybrids and their female parents have also been reported in some fish species. For example, Guo et al. [14] found that complete mtDNA nucleotide identity between the triploid crucian carp and its male parent allotetraploid was higher than that between the triploid crucian carp and its female parent Japanese crucian carp (98% and 93%, respectively).
There are six bream fish species distributed in Chinese natural lakes or rivers, including Megalobrama amblycephala, M. skolkovii, M. terminalis, M. pellegrini and M. elongata belonging to Megalobrama genus, as well as white Amur bream (Parabramis pekinensis) in Parabramis genus, the sister genus of Megalobrama [15,16]. Because of their economically valuable traits, bream fishes have been considered as main aquaculture fish species in China since 1960s [17]. In order to breed superior culture strains, hybridization has been conducted among these fish species [18]. Xie et al. [19] conducted the hybrid breeding of M. hoffmanni (♀) × M. amblycephala (♂) and discovered that the first generation hybrid featured high survival rate and was fertile. It resembled M. hoffmanni in flesh quality and had advantages over M. hoffmanni in resistance to hypoxia. The previous study had also detected that the hybrids of M. amblycephala (♀) × P. pekinensis (♂) exhibited better disease resistance [20].
In this study, we sequenced the complete mitogenomes of 7 hybrid groups of the bream species for the first time. Along with the 10 mitogenomes reported previously by our group, a total of 17 mitogenomes were analyzed in the study, representing 5 bream species, 6 first generation hybrids and 6 second generation hybrids. We explored the way the mitogenomes are inherited during hybridizations among different bream species and how the variability of the studied mitogenomes is distributed in mitochondrial genes of hybrids in comparison to their parental species.

Ethics statement
All experiments were conducted following the ''Guidelines for Experimental Animals" of the Ministry of Science and Technology (Beijing, China). The study was approved by the Institutional Animal Care and Use Ethics Committee of Huazhong Agricultural University. All efforts were made to minimize animal suffering.

Samples and DNA extraction
The parents of each crossing combination were all from wild population and the hybrid breeding artificial reproduction was conducted in 2012 (first generation, F 1 ) and 2015 (second generation, F 2 ) in Fish Breeding Base of College of Fisheries (Hubei Bai Rong Improved Aquatic Seed CO., LTD, Huanggang, 438800), Huazhong Agricultural University, Hubei province of China. Specimens of hybrid offspring were collected randomly from their population. Total genomic DNA was extracted from the fin tissue using a modified ammonium acetate precipitation protocol [21]. The bream fish types and cross combinations are showed in Table 1. Ten individuals from each group were used for the study.

Primer design, PCR amplification and sequencing
Twenty-four pairs of PCR primers (Table 2) were designed to amplify the whole mitogenome sequences based on the conserved sequences from M. amblycephala (GenBank NC_010341.1), using Primer 6.0 and Oligo 7.0 software. The primers were synthesized by Life Technologies Biotechnology Company (Shanghai, China). The amplifications were performed in 25 μL reaction volume containing 1× LA PCR buffer II (Mg 2+ ), 1.25 mM of dNTPs, 0.5 mM of each primer, 1.25 U of LA Taq polymerase (Takara, Dalian, China), approximately 100 ng of template genomic DNA. PCR was performed under the following conditions: denaturation at 95°C for 5 min, followed by 30 cycles at 98°C for 10 s, 52-58°C for 45 s and 72°C for 1 min, as well as further incubation for 10 min at 72°C. Subsequently, the PCR products were purified and directly sequenced by Quintarabio Biotechnology Company (Wuhan, China).

Gene annotation and sequence analysis
After sequencing, the sequence fragments were edited by DNASTAR 5.0 software to obtain the complete mitogenome sequences. Annotation of protein-coding genes (PCGs), ribosomal RNA genes, transfer RNA genes and definition of their respective gene boundaries were performed by MitoAnnotator software (http://mitofish.aori.u-tokyo.ac.jp/annotation/input.html).
The secondary structures of the tRNA were predicted by tRNAscan-SE1.21 software [22]. Sequence alignment was carried out with the Clustal W (http://www.ebi.ac.uk/Tools/msa/ clustalw2/). In order to know the selective pressure acted on each protein-coding gene in these species, we calculated the non-synonymous and synonymous substitution rates (dN/dS) using YN00 in PAML software [23,24] between all the hybrids with their own female parents.

Phylogenetic relationship
According to the phylogenetic tree (Fig 1), all bream fishes and their hybrids were clustered together. The three outgroup species (A. nigrocauda, C. alburnus and C. mongolicus) were grouped into a single clade. The phylogenetic analyses clustered the 4 bream fishes from the genus Megalobrama and 1 from the genus Parabramis into the same clade, which confirmed the close taxonomic relationship between these two genera. Among the Megalobrama species, MS and MP have the closest phylogenetic relationship, with MA being more closely related to them than MT. This result is consistent with the traditional taxonomic relationship of these species reported in previous studies [26,27]. According to the mitogenomes, each hybrid (F 1 and F 2 groups) was more closely related to its own female parent.

Mitogenome organization and composition
Genome length, AT-richness and base composition of the 17 mitogenomes are reported in Table 3. The lengths of the complete genome were all in the range of 16,621 to 16,623 bp, and the genome lengths of the hybrids corresponded with that of the female parents. The overall base composition was slightly AT-rich (Table 3). Each genome contained the same 22 transfer RNA genes, 13 PCGs, 2 ribosomal RNA genes, and 2 main non-coding regions of the control region and the origin of the light strand replication. Most genes of all 17 mitogenomes were encoded on heavy strand (H-strand) except for 1 PCG (the nd6 gene) and 8 tRNA genes (tRNA-Gln, tRNA-Ala, tRNA-Asn, tRNA-Cys, tRNA-Tyr, tRNA-Ser, tRNA-Glu, and tRNA-Pro), which were encoded on light strand (L-strand) (Fig 2). Structurally, all mitogenomes tested in this study shared the same gene arrangement and also displayed conserved genomic arrangement with other teleost species [28][29][30].

Sequence similarity and variation
Sequence comparisons of the mitogenomes between the hybrids and their female parents revealed sequence similarities above 99% in all cases (Table 4). In the F 1 , the highest similarity (99.8797%) was detected between the MA×MP hybrid and its female parent MA, while the lowest similarity (99.4826%) was found between MT×MA hybrid and its female parent MT. A certain number of variable sites were detected between all F 1 hybrids and their respective female parent. The results showed that reciprocal crosses between two breams can generate hybrids with different degrees of mitogenome variation. MT×MA had 86 variable sites, while the reciprocal cross MA×MT only had 23 variable sites. MA×MS had fewer variable sites (21) than MS×MA (38). There is now a broad consensus that reciprocal crosses between species can yield hybrids with different performance [31,32]. In centrarchid fishes, of the 18 species pairs with reciprocal crosses, 17 pairs showed asymmetrical viabilities [31]. The hybrids of Pogonias cromi (♀) × Sciaenops ocellatus (♂) were found to be viable, while the fertilization rate for the reciprocal hybrids from S. ocellatus (♀) × P. cromi (♂) was 0% [33]. Kim et al. [34] had reported that the hatching rates were significantly different in the reciprocal crosses between Japanese flounder (Paralichthys olivaceus) and spotted halibut (Verasper variegatus). Nuclear-cytoplasmic interaction (especially the interaction between mitochondria and nucleus) was regarded as an important cause for asymmetric performance during reciprocal crosses [32,35]. Whether the different variations of mitogenomes detected in the present study could be related to some effects of reciprocal hybrids among bream species, for example to how mito-nuclear interaction works different in the reciprocal hybrids, needs further investigation. The variation of mitogenomes in F 1 hybrids with M. amblycephala as female parent was distinctly lower than those with M. terminalis or M. skolkovii as female parent, which might indicate that M. amblycephala mitogenome interacts better with a hybrid nuclear background when compared to M. terminalis and M. skolkovii. In addition, there were more mutation sites in MT×MA than MS×MA. This is in accordance with their genetic relationship, as MA has  closer phylogenetic relationship with MS than with MT [26,27]. These results may indicate that the genetic distances between parents may have an impact on the mitogenomes of the hybrids, maybe allowing a more stable mitogenome inheritance when the hybrid nucleus derives from closer related species. Regarding F 2 generation hybrids, alignment analysis revealed great similarity among the 6 groups' mitogenomes ( Table 4). The mitogenome sequences of (MA×PP)×MA, (MA×PP) ×(MA×PP), (MA×MT)×(MA×MT), (MA×MT)×(MA×PP) and (MA×PP)×(MA×MT) were completely the same. The highest sequence similarity occurred between the (MA×MT)×MA and its female parent, with 2 variable sites (1 in 16s RNA and 1 non-synonymous mutation in the nd6 gene). Four variable sites were observed in (MA×PP)×(MA×PP) and (MA×MT)×(-MA×MT) compared to their female parents, which were generated by the self-cross of MA×PP and MA×MT, respectively. There was no obvious difference between the hybrids obtained by reciprocal crosses (Table 4). Compared with the variation in F 1 groups, F 2 groups had relatively fewer variable sites. All the results mentioned above reflect that hybridization has a weaker influence on the mitogenome sequences in F 2 than in F 1 , maybe related to the fact that both the parents of the F 2 derived from MA strains.

Protein-coding genes
Of the 213 variable sites found in the F 1 groups, 150 (70.42%) were located in the mitochondrial PCGs. The nd5 gene displayed the highest variability, with 6 non-synonymous mutations and 45 synonymous mutations. On the other hand, the cox2 and atp8 showed the lowest variation rates, with no variation sites found. The highest variability occurred in the nd6 gene in the F 2 groups (Fig 3), while there were no variable sites found in the nd2, cox2, atpase8, atpase6, nd3, nd4l, nd5 and cytb genes.
Like other vertebrate mtDNAs [36,37], most variable sites were at the third codon position, resulting in synonymous mutations (Fig 4). The non-synonymous/synonymous rate ratio (dN/ dS) is widely used as an indicator of selective pressure at the sequence level among different species. It is commonly accepted that dN > dS, dN = dS, and dN < dS generally indicate positive selection, neutral mutation, and negative selection, respectively [23]. In the present study, all dN/dS values of the most PCGs (except the cox2 and atp8 genes) are less than 1 between F 1 hybrids with its own female parent (Table 5)  There were no variation site found in the cox2 and atpase8 genes of F 1 hybrids, and the Dloop region, tRNA, 12S rRNA, nd2, cox2, atpase8, atpase6, nd3, nd4l, nd5 and cytb genes of F 2 hybrids. The variation rates here refer to the average variation rates.
doi:10.1371/journal.pone.0158915.g003 (Table 5). All the 19 variable sites found in 13 PCGs of all the F 2 hybrids were non-synonymous mutations; therefore, the dN/dS was not calculated between the F 2 hybrids and their own female parents. Non-synonymous substitutions are more strongly affected by natural selection than synonymous substitutions and fixations of slightly deleterious mutations are expected to increase the non-synonymous substitution rate [38]. There may be a fixation of slightly deleterious mutations, leading to an increase of non-synonymous substitution in F 2 hybrids.
Some previous studies supported that mitochondrial variants may contribute to phenotypic variation in poultry and livestock. In cattle, mtDNA variation has been associated to economic traits such as milk yield, calving rates, weight and so on [39][40][41]. Fernandez et al. [42] justified the use of the polymorphism in the cytb gene as a marker for maintaining an adequate intramuscular fat level in Iberian pigs. In the present study, the nd5 gene displayed the highest variability in the F 1 groups. Notably, 48 variation sites were found in the nd5 gene between MT×MA and its female parent MT in the present study. The polymorphism of the mtDNA nd5 gene had been reported to be significantly associated with growth traits at 6 months in the  cattle [43]. Sutarno et al. [40] also found a significant association between calving rate and mitochondrial polymorphisms of two regions (the D-loop and the nd5 gene) in two cattle breeds. What is the likely cause of the association of mtDNA nd5 gene polymorphism with phenotypes? Firstly, the mutation rate of mtDNA is higher than that of nuclear DNA [2]. The rapid rate of mutation in mtDNA makes it possible to produce advantageous or disadvantages phenotypes in a relatively short time [43]. Secondly, genetic variation affecting trans-acting nuclear factors can alter the mtDNA sequence and affect phenotypes, so associations between some nuclear genomic regions and traits could be mediated through the effect of mtDNA variation [43,44]. In addition, the enzyme coded by the nd5 gene play a vital role in respiratorychain activities, thus influencing energy supply [45,46]. Whether the high variation rate of the nd5 gene identified in the F 1 groups could be related to phenotypic variation between MT×MA and MA, needs further investigation.
Except for cox1, which begins with GTG, all other PCGs initiate with the classical ATG start codon. In most hybrids except MT×MA and MS×MA, the stop codons of the 13 PCGs include 7 TAA codons and 6 incomplete stop codons, with 3 TA-(the atpase6, cox3 and nd4 genes) and 3 T-(the nd2, cox2 and nd3 genes). In MT×MA and MS×MA, the stop codon of the cytb gene is a single T (T-), and the stop codons for the other 12 PCGs are the same as other species. It seems that this kind of incomplete termination codon could be fixed by post-transcriptional polyadenylation [47], which is common in vertebrate mitogenomes [37,48,49]. There are 2 overlapping reading-frames on the same strand and 1 on the opposite strand. The atp6 and atp8 shared 7 bp, the nd4l overlapped 7 bp with the nd4 and the nd5 shared 4 bp with the nd6, which was coded on the opposite strand (Table 6). Gene overlapping is common in mitogenomes of other vertebrate mitogenomes [48,[50][51][52].

Ribosomal and transfer RNA genes
All mitogenomes contain two rRNA subunits, 12S rRNA and 16S rRNA, which are separated by tRNA-Val as in the other vertebrates [37]. Twenty variable sites were detected in the rRNA genes of F 1 groups, with 2 in 12S rRNA and 18 in 16S rRNA. A total of 26 variable sites were observed only in five tRNA genes (4 in tRNA-Phe, 8 in tRNA-Val, 4 in tRNA-Ile, 1 in tRNA--Tyr and 9 in tRNA-His). In the F 2 generation, 3 variable sites were detected in 16S rRNA and none in 12S rRNA. There were no variable sites found in 22 tRNA genes of all the F 2 hybrids as well. Twenty-one tRNA genes could be folded into the typical cloverleaf secondary structure [53], with the exception of tRNA-Ser, because the tRNA-Ser (AGY) lacks a DHC arm [40]. In MA×MT, MA×PP, MA×MT and MA×MP, a variable site in the anticodon was observed in tRNA-Ile. Among tRNA genes, there were 2 overlapping reading-frames on the opposite strand. Two nucleotides overlapped between tRNA-Ile and tRNA-Gln, and one between tRNA-Thr and tRNA-Pro (Table 6) were found in the mitogenomes of all bream and hybrids.

Non-coding regions
The major non-coding region (D-loop), located between tRNA-Pro and tRNA-Phe genes, was 937 bp in length for most hybrids, except it was 936 bp in MT×MA. As reported in other fish species [54,55], 3 conserved domains were identified by multiple homologous sequence alignment, consisting of extended termination associated sequences (ETAS), a central conserved domain (CCD) and conserved sequence blocks (CSB). The terminal associated sequence (TAS) was observed in the first domain with 4 variable sites. It was regarded as the most variable region in the D-loop, which was consistent with other fish species [54]. The central conserved domain was recognized as the most conservative region in the D-loop, containing CSB-F, CSB-E and CSB-D. The third domain, located at the 3' end of D-loop, was comprised of CSB1, CSB2 and CSB3. The D-loop region is known as an AT-rich region and our data also supported this observation [54]. Based on the structure of the D-loop region reported previously in other fish species [27,[54][55][56][57], the conserved motifs in D-loop region of the 17 groups in this study were identified (Table 7). A relatively lower rate of substitution in the D-loop region was found in our study. The mutation rate of the D-loop region (0.30244%) was lower than those of the nd5 (0.46296%) and nd2 (0.41467%) genes (Fig 3) in the F 1 groups. In the F 2 groups, there were no variable sites found in the D-loop region, the mutation rate of which was lower than those of the nd1 (0.05128%), cox1 (0.03224%), cox3 (0.12739%), nd4 (0.02412%) and nd6 genes (0.15964%) (Fig 3).
A non-coding region of 32 bp, the origin of light strand replication (oril) [58], is encoded on L-strand and located between tRNA-Asn and tRNA-Cys. This region has a palindromic sequence, so it can be folded into a stem-loop secondary structure [27]. The conserved motif 5'-GCCGG-3', which was often observed in the tRNA-Cys gene of many vertebrates' mitogenomes [29,59,60], was also detected in our study.

Conclusions
Our results indicated that, although all the mitogenomes of the hybrids were consistent with maternal inheritance, a certain level of variation was detected in the mitogenomes of F 1 and F 2 hybrids, with MT×MA and MS×MA groups having relatively more variations. The second generation hybrids showed less mitogenome variation than that of first generation hybrids. The most variable gene in F 1 groups was the nd5, while the nd6 displayed the highest variability in F 2 groups. The number of variation sites was found to be related to phylogenetic relationship of the parents. The dN/dS analysis showed that the most PCGs (except the cox2 and atp8 genes) were under negative selection. The information reported in this study could be useful for further studies on mitogenome inheritance and variation in fish hybrids.