Maternal ancestry analyses of red tilapia strains based on D-loop sequences of seven tilapia populations

Sample collection

The sampling scheme and experimental protocols were approved by the Bioethical Committee of Freshwater Fisheries Research Center (FFRC) of the Chinese Academy of Fishery Sciences (CAFS) (BC2013863, 9/2013). Experimental fish were sampled from seven populations, including the three red tilapia strains of Oreochromis spp., Chinese Taiwan (TW), Israel (IL), and Malaysia (MY) strains, and other populations of tilapia, including the GIFT strain of O. niloticus, O. aureus (AR), O. niloticus (NL), O. mossambicus (MS). The TW and IL populations were transferred from Fujian Province, China in 2014, and the MY population was introduced from Malaysia in 2009 (Yang et al., 2015). All tilapia populations were domesticated and bred in an experimental aquaculture farm in Wuxi (Jiangsu Province, China). The methods for handling the animals and the experimental procedures conducted were in accordance with the guidelines for the care and use of animals for scientific purposes set by the Ministry of Science and Technology, Beijing China (No.398, 2006). Forty-eight fin clips were sampled from each population and soaked in absolute ethanol until the DNA was extracted.

DNA extraction and amplification

Genomic DNA was extracted using the phenol-chloroform method (Sambrook & Russell, 2001). The integrity was detected by 1% agarose gel electrophoresis. The purity and concentration of the DNA was detected using the NanoDrop spectrophotometer. The DNA concentration of each sample was adjusted to about 20 ng/µL and kept under −20 °C until ready to use.

Primers of the D-loop were designed according to the complete sequences of tilapia mtDNA (Accession NO: NC_014060) from the National Center for Biotechnical Information (NCBI). The D-loop sequences of 867 bp was amplified using a primer pair (sense primer: 5′-CTACTTCTTCCTCTTCCTTGT-3′, anti-sense primer: 5′-TCCGTCTTAACATCTTCAGT-3′), which was synthesized by Sangon Biotech (Shanghai) Co.Ltd. The PCR amplification was performed on an Eppendorf Mastercycler Pro 384 PCR thermocycler (Eppendorf, Hamburg, Germany). Amplifications were performed in a volume of 50 µL, containing 5 µL 10 × PCR Buffer, 3 µL MgCl₂(0.25 mM), 4 µL dNTPs (2.5 mM), 1 µL Taq polymerase (2.5 U/µL), 1 µL of each primer (10 µM), 2 µL genomic DNA (20 ng/µL), and 33 µL DNase/RNase-free deionized water. PCR amplification was performed under the following conditions: pre-denaturing for 2 min at 94 °C, 35 cycles of denaturing for 40s at 94 °C, annealing for 55s at 55 °C, prolonging 1 min at 72 °C; final prolonging for 10 min at 72 °C; and then held at 12 °C. Subsequently, the reaction products were detected using 1% agarose gel electrophoresis, and the bidirectional sequencing was carried out with the ABI3730XL sequencing instrument of the Shanghai Majorbio Company.

Sequence arrangement and data analysis

The sequences were edited using the BioEdit version 7.0.9 software (Hall, 1998). To ensure accuracy, all DNA fragments were sequenced in two directions, and the assembled sequences were manually checked to prevent the ambiguity of the base or sequencing error. After the completion of the splicing, all sequences were used for homologous alignment and length determination by the BioEdit version 7.0.9 software.

After comparing and sorting the D-loop sequences, the fuzzy sequences were deleted and 332 homologous sequences were obtained upon completion. Genetic variation parameters were calculated by the DnaSP 5.1 software (Librado & Rozas, 2009), including polymorphic (segregating) sites (S), number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), average number of nucleotide differences (k) and Tajima’s D. For phylogenetic analysis, MEGA 5.05 software (Kimura, 1980; Tamura et al., 2011) was utilized to calculate the Kimura 2-parameter (K2P) distance among populations, and to construct unweighted pair group methods with an arithmetic (UPGMA) dendrogram set with 1,000 replications of bootstrapping. Arlequin 3.5 software (Excoffier & Lischer, 2010) was used to analyze the nucleotide composition, F-statistics (F_ST) and the analysis of molecular variance (AMOVA) among seven tilapia populations. The Network 4.6 software (Polzin & Daneshm, 2003) was used to construct the network for haplotypes of D-loop sequences.

Variation and haplotype distribution of D-loop sequences in tilapia

The nucleotide frequencies of the seven tilapia populations were consistent, and clearly, the rate of A + T (the average value is 64.3%) was higher than in C + G (the average value is 35.7%). By sequence alignment, a total of 32 haplotypes were found in the D-loop sequences and deposited in the GenBank database under the accession numbers MH515150–MH515185 (except for MH515152, MH515172, MH515175, and MH515182). Different numbers of haplotypes (from 1 to 20) were detected among the populations (Table 1). Four of these haplotypes were shared haplotypes (Hap_2, Hap_22, Hap_23, and Hap_24), of which 2 haplotypes were composed of NL and GIFT populations (Hap_22, Hap_23), while the others were unique to each population. The 4 dominant haplotypes in all individuals were Hap_2, Hap_22, Hap_24, and Hap_26 accounting for 34.90%, 12.30%, 14.50%, and 14.50%, respectively.

Genetic diversity and genetic distance among seven tilapia populations

The genetic diversity parameters of the tilapia populations based on D-loop sequences are shown in Table 2. A total of 146 polymorphism sites were found, the overall haplotype diversity (Hd) of the tilapia populations was 0.817, and each population ranged from 0 to 0.834. The average number of nucleotide differences (K = 0 − 47.32) and nucleotide diversity (π = 0 − 0.060) were determined. Among them, the AR population had the lowest genetic diversity (Hd = 0, π = 0), the TW and NL populations had a higher haplotype diversity (Hd = 0.834, 0.826, respectively), and the highest nucleotide diversity (π) was detected in the NL population (π = 0.060). Tajima’s test indicated that the Tajima’s D value of the TW, GIFT, and MS populations were negative, and other populations were positive. Among them, NL and MS reached a significant level (P < 0.05), GIFT, IL and MY populations reached extremely significant levels (P < 0.01).

The pairwise genetic distance was calculated using the Kimura 2-parameter (K2P) model among seven tilapia populations (Table 3, below diagonal). The inter-population distances among seven populations ranged from 0.014 to 0.101. The closest inter-population distance (0.014) was between the MS and TW populations, while the furthest inter-population distance (0.101) was between the GIFT and AR populations. In this study, the UPGMA tree based on the K2P genetic distance was shown in Fig. 1. According to the phylogenetic tree, the haplotype was obviously divided into two branches (the AR population was separated from other six populations). In addition, the MS and TW populations were clustered and then clustered with the IL population; GIFT and MY populations were clustered and then clustered with the NL population.

Genetic differentiation among tilapia populations

The results of the analysis of molecular variance (AMOVA) were shown in Table 4. Based on the results of the genetic differentiation analysis, the variance percentage of genetic variation among populations in total variance was 63.57%, and a high degree of inter-population differentiation was observed among populations (Fst = 0.636, P < 0.01). The seven populations were divided into two groups consisting of wild-type or breeding populations, and red tilapia. The genetic differentiation index among populations within groups accounted for 0.633 (P < 0.01). The pairwise F-statistics values (Fst) of the seven tilapia populations (Table 3, upper right corner) showed that the pairwise genetic differentiations among the populations were very significant (P < 0.01), except among the MY and IL populations (P > 0.05).

Haplotype network of D-loop sequences in tilapia populations

The median-joining (MJ) network diagram of D-loop haplotypes was described in Fig. 2. The MJ network presented a star-like profile, which was linked to many haplotypes from different regions, and the shared haplotypes and dominant haplotypes were clearly defined. Obviously, three shared haplotypes were composed of two populations (GIFT and NL, MY and IL), one shared haplotype was composed of six populations (except AR) and a dominant haplotype was composed of an AR population.

Maternal ancestry of red tilapia strains

According to the genetic distances among the tilapia populations (Table 3, below diagonal), it was calculated that TW, IL red tilapia, and MS populations had the closest genetic distance, which were 0.014 and 0.032, respectively. The genetic distance between MY red tilapia and the GIFT population was the closest (0.034). Moreover, the UPGMA tree clearly divided the tested samples of red tilapia into two independent branches. TW and IL were clustered into one branch, and then clustered with MS, while MY divided into another branch and clustered with GIFT and NL.

Genetic diversity and population dynamics

In the present study, the content of A+T (64.3%) in tilapia D-loop sequences was higher than the content of G+C (35.7%), which was in line with the distribution characteristics of the base content in the D-loop (control region) of many fishes (Broughton, Milam & Roe, 2001). In D-loop sequences of seven tilapia populations, 146 polymorphic (segregating) sites (S) and 32 haplotypes were detected, suggested that D-loop sequences could be an effective marker for genetic diversity analysis for tilapia populations. Overall, the tilapia populations showed high haplotype diversity and nucleotide diversity in this study, indicating the populations that contain an abundance of genetic resources for further use in breeding or practice. Specifically, NL had the higher genetic diversity (Hd >0.5, π >0.005), which was consistent with other research using microsatellite markers (Romana-eguia et al., 2004; Yang et al., 2011), and isozymes (Zhao et al., 1997). It was speculated that the original NL population introduced in China was larger and it had a potential for further selective breeding. The genetic diversity of the red tilapia populations (TW, IL, and MY) was higher than the other tilapia populations, probably due to its genetic background of cross-breeding. However, the genetic diversity of the AR population was the lowest (Hd = 0, π = 0), which was found to be similar to the results of the previous reports assessed by RAPD (Xia et al., 1999), microsatellite probes (Wang, Xia & Wu, 2000), and mtDNA restriction enzyme analysis (Cao & Xia, 1997). This may be due to the small population size introduced in China as well as mass breeding over generations, which might result in a decline in the genetic polymorphism of the population. The genetic purity of the AR population was adverse to further selective breeding and it is necessary to introduce the AR population again in order to improve its genetic diversity and reduce its depression as a result of inbreeding. In contrast, the purity of this population also could be used for inhybridization with other populations or strains. The GIFT and MS populations were also detected with low diversity parameters (Hd < 0.5, π < 0.005), indicating that its population may have recently experienced a bottleneck or founder effect produced by minority populations (Grant & Bowen, 1998). The high purity of the AR, GIFT, and MS populations has great significance that the excellent economic traits obtained through long-term multi-generation breeding could insure the stability inherited through the genetic selection process. Therefore, AR, GIFT, and MS were often used as parents to breed red tilapia and stabilize the excellent traits in the red tilapia populations. The Tajima’s D value for some red tilapia (TW) and wild-type or breeding populations (MS, GIFT) were negative, which is probably due to the population expansion caused by larger scale breeding after smaller breeder population produced a family selection in the hatchery.

Genetic relationships of tilapia populations

The red tilapia stocks were reported with different genetic sources that originated from O. mossambicus, O. aureus, or O. niloticus (Wohlfarth et al., 1990; Sandeep et al., 2012). In this study, four relative tilapia populations were used for exploring the maternal ancestors for three strains of red tilapia. The MJ network of D-loop sequences for haplotypes in the tilapia populations, which were divided into three different major regions and red tilapia populations, presented in two regions. Four dominant haplotypes were shared by the seven populations, where Hap_24 was a shared by MY and IL populations, Hap_22 and Hap_23 were shared by GIFT and NL populations, and Hap_2 was shared by six populations except for the AR population. It was suggested that six populations (except AR) may originate from similar maternal ancestors.

The analysis of molecular variance (AMOVA) by grouping and non-grouping showed that the main genetic variation was derived from inter-population, which is similar to the results of Habib et al. (2011), the low variance within the population and high variance among populations was reported among Channa fishes. The fixation index (Fst) was commonly used to examine the genetic variation of populations and the contribution of this variation to genetic differentiation (Holsinger & Bruce, 2009). Significant Fst values (Fst >0.25, P < 0.01) were found in this study, which demonstrated that a higher level of genetic differentiation exists among tilapia populations except for red tilapia. The results indicated that red tilapia populations may have evolved independently after separating from a common ancestor, but that those three strains of red tilapia were closed.

Analysis of maternal ancestry of red tilapia

While the maternal ancestors of the existing three strains of red tilapia are not well documented, their derivation is generally attributed to the crossbreeding of the mutant reddish-orange Oreochromis mossambicus with other species including O. aureus, O. niloticus, and O. hornorum (Wohlfarth et al., 1990; Sandeep et al., 2012). Mitochondrial DNA (mtDNA) has the characteristics of maternal inheritance, simple structure, and rapid evolution. Therefore, the phylogenetic tree constructed by mtDNA can directly reflect the origin of the maternal ancestry (Cann, Stoneking & Wilson, 1994).

Based on the K2P genetic distances among seven tilapia populations, two branches were constructed in the UPGMA dendrogram and it was speculated that these populations might be derived from two different primary maternal ancestors, which was consistent with the results of the MJ network. Specifically, three strains of red tilapia might derive from different maternal origin, MS and GIFT (or NL) populations, respectively. The TW and IL populations were closely related to each other and were clustered with MS, which was confirmed that two strains of red tilapia were produced from local crossbreeding of the rare, mutant-colored (reddish-orange) female O. mossambicus (Wohlfarth et al., 1990). The GIFT strain was selected from four O. niloticus strains imported directly from Africa and four strains widely cultivated in Asia (Eknath et al., 1993). The K2P genetic distance between MY red tilapia and the GIFT population was relatively small (D = 0.034), speculating that MY population was probably bred with the GIFT population or that MY and GIFT populations might come from the similar, artificially selected NL population. In addition, the degree of genetic differentiation between red tilapia (IL, TW, and MY) and the breeding source populations (MS, GIFT or NL) were relatively small, demonstrating that a close genetic relationship was maintained between the breeding varieties and breeding source populations; this highly homology was related to the characteristics of maternal inheritance and non-recombination of mtDNA (Mabuchi, Senou & Nishida, 2010).