# High Sequence Variations in Mitochondrial DNA Control Region among Worldwide Populations of Flathead MulletMugil cephalus

Authors: Brian Wade Jamandre; Jean-Dominique Durand; Wann-Nian Tzeng
Journal: International Journal of Zoology
 (2014)

---

## Abstract

The sequence and structure of the complete mtDNA control region (CR) ofM. cephalus from African, Pacific, and Atlantic populations are presented in this study to assess its usefulness in phylogeographic studies of this species. The mtDNA CR sequence variations among M. cephalus populations largely exceeded intraspecific polymorphisms that are generally observed in other vertebrates. The length of CR sequence varied among M. cephalus populations due to the presence of indels and variable number of tandem repeats at the 3
′ hypervariable domain. The high evolutionary rate of the CR in this species probably originated from these mutations. However, no excessive homoplasic mutations were noticed. Finally, the star shaped tree inferred from the CR polymorphism stresses a rapid radiation worldwide, in this species. The CR still appears as a good marker for phylogeographic investigations and additional worldwide samples are warranted to further investigate the genetic structure and evolution in M. cephalus.


---

## Body

## 1. Introduction

The advent of next generation sequencing started an era of discovery for a wide array of genetic markers for population genetics, phylogeny, and so forth. However, we have only scratched the surface on the analyses of these new markers and most researchers still resort to using the mitochondrial DNA (mtDNA) for population and biodiversity studies. Furthermore, the nonrecombining nature due to maternal inheritance and fast evolution rate compared to the nuclear genome [1] makes mtDNA a popular marker for studying vertebrate population genetics [2–4]. Within the mitochondrial genome, evolutionary rates are highly variable among genes. For example, in salamanders, nucleotide substitution rates may range from 0.1 to 1 substitution/site/million years (MY) [5]. Among the coding genes, the cytochromeb (cytob) and cytochrome oxidase I (COI) present the highest evolutionary rates; however, the noncoding region of the mitochondrial genome (i.e., the control region (CR)) evolves 2–5 times faster [6]. Because of this exceptionally high mutation rate characteristic, the CR was quickly considered of great utility for addressing intraspecific evolutionary questions [7, 8].In Teleostei, the CR is located betweentRNA
Pro and tRNA
Phe. Structurally, the CR is composed of a central conserved domain flanked by two highly variable (left and right) domains (Figure 1). In these large domains, several short (approximately <30 bp) but conserved sequenced blocks (CSBs) have been described and varied in occurrence between taxa [2]. In the central conserved domain, the CSB-D is involved in the heavy strand replication, including initiation by a 3-strand displacement loop (D-loop) [9]. The CSB-I, CSB-II, and CSB-III are located in the right domain. The functions of these CSBs are not well understood but in mammals the transition to DNA synthesis occurs in the region surrounding CSB-II, which is also evident in flounder fishes [2]. Within the regions of the CR, the tRNA
Pro end (extreme 5′ end) is a popular chosen marker for evolutionary and population genetics studies because of its rapid evolutionary rate compared to other mitochondrial genes and nuclear coding regions [2, 10]. In fishes, this region is also widely used for inter- and intraspecies investigations, but its hypervariability cannot be generalized to all taxa [10]. For example, in some salmonids (brown trout) the tRNA
Pro end of the mtDNA CR seems to evolve at a lower rate [11, 12]. Indeed, differences in sequence patterns and evolutionary rate of the CR are evident but poorly understood [2, 10].Figure 1
The mtDNA CR structure ofMugil cephalus. Arrows indicate locations of primer binding sites. Broken lines indicate location of variable number of tandem repeats (CSB—conserved sequence block).Intraspecific diversity recovered from the CR is quite variable between fish taxa and among them, the cosmopolitan flathead mullet,Mugil cephalus, which occurs between 42°N–42°S latitudes mostly in the tropical and subtropical areas of the world [13], shows exceptionally high CR diversity [14–16]. Between remote populations from Atlantic and Pacific oceans, Rocha-Olivares et al. [15] showed a level of divergence ranging from 38 to 75%. More recently, Jamandre et al. [16] observed CR haplotypes with sequence divergence of 48% between sympatric populations ofM. cephalus in continuous sampling locations along the northwest Pacific. Such exceptional level of sequence divergence brought Rocha-Olivares et al. [15] to question the origin of the CR amplified in their study, as nuclear pseudogenes stemming from the duplication of mtDNA regions can greatly diverge from their mitochondrial copies. However, by sequencingM. cephalus CR they successfully demonstrated the mitochondrial origin of all haplotypes observed in Pacific and Atlantic populations of this species. The levels of divergence observed might be due to an accelerated evolutionary rate in theMugil lineage as suggested earlier by Caldara et al. [17]. In this context, CR would be an interesting marker to investigate the worldwide phylogeographic structure ofM. cephalus because all phylogeographic investigations to date, using different markers (mtDNA RFLP, allozymes), failed to highlight any significant genetic structure consistent with the geography [18, 19]. However, before utilizing it in such way, determining the basic properties of this marker for phylogenetic usefulness is warranted, especially by checking the level of homoplasy that may blur phylogenetic signals.In this study, we describe a remarkably high evolutionary rate of mtDNA CR inM. cephalus by analyzing its sequence structure and mutations from previously reported and new worldwide samples. SinceMugil liza from South America belongs to theM. cephalus clade as demonstrated by Fraga et al. [20], Heras et al. [21], and Durand et al. [22], we also included samples of this species in this study.


## 2. Material and Methods

### 2.1. Data Gathering, Sample Collection, and DNA Extraction

To provide a general and worldwide view of sequence variations and genetic diversity inM. cephalus, we gathered all availableM. cephalus CR sequences in GenBank (http://www.ncbi.nlm.nih.gov/) (Table 1). These sequences or haplotypes belong to the differentM. cephalus lineage previously described by Rocha-Olivares et al. [14, 15], Miya et al. [4], and Jamandre et al. [16]. NewM. cephalus samples were also added and analyzed from Senegal (
n
=
1
) and Togo (
n
=
2
) in West Africa and South Africa (
n
=
1
). It was stressed by Fraga et al. [20], Heras et al. [21], and Durand et al. [22] thatMugil liza belongs to theM. cephalus clade. Thus, a sample ofM. liza from French Guyana (and also to represent the southwest Atlantic) was added in this study to further investigateM. cephalus andM. liza relationships.Table 1
Location, number of samples (N), and the GenBank (http://www.ncbi.nlm.nih.gov/) accession numbers of the mtDNA control region sequences of M. cephalus used in this study.
Location
N
Accession number
Source
Northwest Pacific
126
EU663629-3754
Rocha-Olivares et al. [15]
Japan
1
NC003182
Miya et al. [4]
Hawaii (Central Pacific)
19
AF108232-69
Rocha-Olivares et al. [14, 15]
Gulf of Mexico/Atlantic
96
AF108271-352
Rocha-Olivares et al. [14, 15]
South Africa
1
GU120874
This study
Togo/Senegal (West Africa)
3
GU120871-73
This study
Mugil liza (French Guyana)
1
GU120875
This studyFor each specimen collected, muscle tissues and fin clips were taken and preserved in the 95% ethanol-water solution and brought back to the laboratory. Genomic DNA was then extracted from a small piece of either muscle tissue or fin clips, using a commercially available total genomic DNA extraction kit (Bioman Scientific Co. Ltd., Taiwan).


### 2.2. PCR and Sequencing

The primers used to amplify the complete CR are described in Jamandre et al. [16]. After an alignment of sequences from northwest Pacific samples, internal primers were designed to amplify West and South African samples andM. liza (mulCR-IPF: GCATTCTCTCCACAGGGGCC; mulCR-IPR: GGCCCCTGTGGAGAGAATGC); MulCR-IPF and mulCR-IPR were paired to MulPhe [14] and McepCR-F [16], respectively.PCR amplification was carried out in a 25μL total volume containing 0.5 μL DNA template, 0.25 μL taq polymerase, 1 μL of each primer, 2.5 μL PCR buffer, 0.5 μL 10x DNTP, and 19.25 μL ddH2O. Thermal cycling was performed as follows: 2 mins at 94°C, followed by 35 cycles consisting of 30 sec at 94°C denaturation, 45 sec at 58°C annealing, and 45 sec at 70°C extension, and a final elongation of 10 min at 70°C. PCR products were visualized on a 1% agarose gel stained with ethidium bromide. Afterwards, the products were sequenced with forward and reverse primers using dye-terminator cycle sequencing on an ABI 373A sequencer.


### 2.3. Sequence Analyses

CR sequences were aligned using the software MAFFT version 6 [23] with default parameters and manual improvement. Multiple sequence alignment is a crucial step in the analysis of any genomic data especially for sequences with many overlapping and nonoverlapping gaps. As such, the software MAFFT was used to align theM. cephalus CR data set, since it uses both iterative and progressive methods of alignment, which has been proven to be more accurate and consistent in dealing with gaps, compared to other popular alignment software [24].The mtDNA CR domains (5′ and 3′ hypervariable domain and central conserved region), termination associated sequences (TAS), and conserved sequence blocks (CSBs) were identified by alignment with reference to the previously publishedM. cephalus control region sequences [14, 15] and other fish:Crossostoma lacustre [25],Chanos chanos[26],Dascyllus trimaculatus [27],Notopterus notopterus [28],Perca fluviatilis [3],Padogobius nigricans [29], and Sinipercine fishes [30].Intraspecific nucleotide diversity and variability were assessed from all sequence data available. Sequence divergence was estimated accounting for different rates of transition and transversion, unequal base frequencies, and among site variation rates because evolutionary relationships might be blurred by homoplasic mutations (i.e., multiple substitutions) at a specific site [31]. To assess the level of saturation of CR, the number of nucleotide substitutions per site for transitions and transversions was plotted against genetic divergent distances. In this case the Tamura and Nei’s (TrN-G) model was used. Genetic divergences between groups were calculated using the MEGA 4.1 software [32, 33].The appropriate DNA substitution model was obtained via testing alternative models of evolution using the Modeltest software version 3.7 [34]. Bayesian inference phylogenetic tree was constructed with MrBayes v3.1.2 [33] using the model of evolution indicated by Modeltest. Four Metropolis-coupled Markov Chain Monte Carlo chains with 2
×
10
6 generations were sampled every 1,000th step. 500 trees were discarded as burn-in and a consensus phylogram of the 1,500 recorded trees and posterior probability were calculated.


## 2.1. Data Gathering, Sample Collection, and DNA Extraction

To provide a general and worldwide view of sequence variations and genetic diversity inM. cephalus, we gathered all availableM. cephalus CR sequences in GenBank (http://www.ncbi.nlm.nih.gov/) (Table 1). These sequences or haplotypes belong to the differentM. cephalus lineage previously described by Rocha-Olivares et al. [14, 15], Miya et al. [4], and Jamandre et al. [16]. NewM. cephalus samples were also added and analyzed from Senegal (
n
=
1
) and Togo (
n
=
2
) in West Africa and South Africa (
n
=
1
). It was stressed by Fraga et al. [20], Heras et al. [21], and Durand et al. [22] thatMugil liza belongs to theM. cephalus clade. Thus, a sample ofM. liza from French Guyana (and also to represent the southwest Atlantic) was added in this study to further investigateM. cephalus andM. liza relationships.Table 1
Location, number of samples (N), and the GenBank (http://www.ncbi.nlm.nih.gov/) accession numbers of the mtDNA control region sequences of M. cephalus used in this study.
Location
N
Accession number
Source
Northwest Pacific
126
EU663629-3754
Rocha-Olivares et al. [15]
Japan
1
NC003182
Miya et al. [4]
Hawaii (Central Pacific)
19
AF108232-69
Rocha-Olivares et al. [14, 15]
Gulf of Mexico/Atlantic
96
AF108271-352
Rocha-Olivares et al. [14, 15]
South Africa
1
GU120874
This study
Togo/Senegal (West Africa)
3
GU120871-73
This study
Mugil liza (French Guyana)
1
GU120875
This studyFor each specimen collected, muscle tissues and fin clips were taken and preserved in the 95% ethanol-water solution and brought back to the laboratory. Genomic DNA was then extracted from a small piece of either muscle tissue or fin clips, using a commercially available total genomic DNA extraction kit (Bioman Scientific Co. Ltd., Taiwan).


## 2.2. PCR and Sequencing

The primers used to amplify the complete CR are described in Jamandre et al. [16]. After an alignment of sequences from northwest Pacific samples, internal primers were designed to amplify West and South African samples andM. liza (mulCR-IPF: GCATTCTCTCCACAGGGGCC; mulCR-IPR: GGCCCCTGTGGAGAGAATGC); MulCR-IPF and mulCR-IPR were paired to MulPhe [14] and McepCR-F [16], respectively.PCR amplification was carried out in a 25μL total volume containing 0.5 μL DNA template, 0.25 μL taq polymerase, 1 μL of each primer, 2.5 μL PCR buffer, 0.5 μL 10x DNTP, and 19.25 μL ddH2O. Thermal cycling was performed as follows: 2 mins at 94°C, followed by 35 cycles consisting of 30 sec at 94°C denaturation, 45 sec at 58°C annealing, and 45 sec at 70°C extension, and a final elongation of 10 min at 70°C. PCR products were visualized on a 1% agarose gel stained with ethidium bromide. Afterwards, the products were sequenced with forward and reverse primers using dye-terminator cycle sequencing on an ABI 373A sequencer.


## 2.3. Sequence Analyses

CR sequences were aligned using the software MAFFT version 6 [23] with default parameters and manual improvement. Multiple sequence alignment is a crucial step in the analysis of any genomic data especially for sequences with many overlapping and nonoverlapping gaps. As such, the software MAFFT was used to align theM. cephalus CR data set, since it uses both iterative and progressive methods of alignment, which has been proven to be more accurate and consistent in dealing with gaps, compared to other popular alignment software [24].The mtDNA CR domains (5′ and 3′ hypervariable domain and central conserved region), termination associated sequences (TAS), and conserved sequence blocks (CSBs) were identified by alignment with reference to the previously publishedM. cephalus control region sequences [14, 15] and other fish:Crossostoma lacustre [25],Chanos chanos[26],Dascyllus trimaculatus [27],Notopterus notopterus [28],Perca fluviatilis [3],Padogobius nigricans [29], and Sinipercine fishes [30].Intraspecific nucleotide diversity and variability were assessed from all sequence data available. Sequence divergence was estimated accounting for different rates of transition and transversion, unequal base frequencies, and among site variation rates because evolutionary relationships might be blurred by homoplasic mutations (i.e., multiple substitutions) at a specific site [31]. To assess the level of saturation of CR, the number of nucleotide substitutions per site for transitions and transversions was plotted against genetic divergent distances. In this case the Tamura and Nei’s (TrN-G) model was used. Genetic divergences between groups were calculated using the MEGA 4.1 software [32, 33].The appropriate DNA substitution model was obtained via testing alternative models of evolution using the Modeltest software version 3.7 [34]. Bayesian inference phylogenetic tree was constructed with MrBayes v3.1.2 [33] using the model of evolution indicated by Modeltest. Four Metropolis-coupled Markov Chain Monte Carlo chains with 2
×
10
6 generations were sampled every 1,000th step. 500 trees were discarded as burn-in and a consensus phylogram of the 1,500 recorded trees and posterior probability were calculated.


## 3. Results

### 3.1. Control Region Sequence Structure and Variability

The observed lengths ofM. cephalus control region varied among sampling locations. The complete CR sequences ofM. cephalus from northwest Pacific ranged from 909 to 1015 bp, whereas in central Pacific (Hawaii), it ranged from 919 to 930 bp. The Atlantic and Gulf of MexicoM. cephalus CR sequences ranged from 884 to 894 bp while the West and South African individuals ranged from 899 to 912 bp and 911 bp, respectively. The CR length of theMugil liza individual is 906 bp.By comparingM. cephalus CR sequences to CR of other fishes species (see Table 2), three domains were highlighted, namely, the left and right hypervariable domain and the central conserved domain (Figure 1). Several conserved sequence blocks (CSBs) were located at the central conserved region and at the right hypervariable domain (Figure 2). CSB-I, -II, and -III are commonly found in vertebrate CR. These were also present in all sequences ofM. cephalus. Putative CSB-D was also found in the sequences (Figure 2). Pyrimidine tract was also identified at the central conserved domain.Table 2
Comparison of conserved sequence blocks (CSBs) identified in the control region ofMugil cephalus to other teleost fishes.
Species
CSB-D
CSB-I
CSB-II
CSB-III
Mugil cephalus
NW Pacific
Lineage 1/
TTCCTGGCATTTGGTTCC
AACATA
AAAACCCCCCCTACCCCCCCTAA
TCCTGAAAACCCCCCCGGAAACA
Lineage 2/
TTCCTGACATTTGGTTCT
AATATA
AAAACCCCCCTACCCCCCTAAA
TCCTGAAAACCCCCCGGAAACAA
Hawaii
TTCCTGGCATTTGATTCC
AACATA
AAACCCCCCCTACCCCCCCAAA
TCCTGAAAACCCCCCGGAAACAA
Atlantic
TCTTTGATGTTCAGGTTTT
AACATA
AAAACCCCCCCTACCCCCCCAAA
TCCTGAAAACCCCCCGGAAACAA
West Africa
TTTCTGGCATTTGGCTCC
AACATA
AAAACCCCCCCTACCCCCCCTAA
TCCTGAAAACCCCCCGGAAACAA
South Africa
TTCCTGGCATTTGGTTCC
AACATA
AAAACCCCCCCTACCCCCCCTAA
TCCTGAAAACCCCCCGGAAACAA
Mugil liza
TTCCTGGCATTTGGTTCC
AACATA
AAAACCCCCCCTACCCCCCCAAA
TCCTGAAAACCCCCCGGAAACAA
Chanos chanos
TTCAAGAGCATCTGGTTCC
—
AAACCCCCCCCTCCCCCCA
TGTTAAACCCCCAAAACCA
Crossostoma
lacustre
TTACTGGCATCTGGTTCC
CTATGTATGTACAATGAGCATAA
ACAAACCCCCCTACCCCCCT
TGCTCAAACCCCGAAACCA
Dascyllus trimaculatus
TTCCTGGCATTTGGTTCC
CACTACTGTCTTCCCGGACATACA
TAAAACCCCCCTACCCCCCT
TGAAACCCCCCGAAAACA
Notopterus notopterus
TATTCCTGGCATCTGGCTCC
—
AACCCCCCTACCCCCCA
TCAAACCCCTAAACCA
Perca fluviatilis
TTCCTGGCATTTGGTTCC
ATCTTAGGATATCAAGAGCATAA
TAAACCCCCCCTACCCCCCC
TGAAAACCCCCCCGGAAACA
Padogobius nigricans
TTTCTGGCATTTGGTTCC
ATAACTGATATCATGAGCATAA
AATCCCCCCTACCCCCCT
TGAAAACCCCCCGGAAACA
Sinipercine
spp.
TATTCCTGGCATTTGGTTCCT
ATTGGATATCAAGAGCATAA
TAAACCCCCCTACCCCCC
TGAAAACCCCCCGGAAACA
—: not identified.Multiple alignment of consensus sequences ofMugil cephalus from L1 (Lineage 1) and L2 (Lineage 2) of the NW (northwest) Pacific, Hawaii/C (central) Pacific, Gulf of Mexico/ Atlantic, Togo and Senegal/W (west) Africa, and South Africa. Structural features (CSB—conserved sequence block and Pyrimidine-tract) are also indicated.
(a)
(b)Insertions and deletions (indels) were found in both the left and right hypervariable domain and generally associated with repetitive sequences (Figure2). The percent composition of nucleotide varied from 35.96 to 37.50% (A), 29.80 to 34.57% (T), 16.88 to 20.64% (C), and 11.29 to 13.27% (G), which infers thatM. cephalus CR is AT rich and poor in G.The transition/transversion (Ts/Tv) rate ratios arek
1
=
20.443 (purines) and k
2
=
14.496 (pyrimidines). An overall Ts/Tv ratio of 5.131 was estimated from the dataset (Figure 3). A linear relationship between the number of nucleotide substitutions per site and the sequence divergence was observed for both, transversion and transitions mutations (Figure 3). Slope of the transitions /divergence regression indicates a higher transition rate than transversion rate.Figure 3
Bayesian tree inferred from complete CR sequences. Posterior probability values (%) are indicated close to the nodes of the branches.


### 3.2. Phylogenetic Relationship and Sequence Divergence

NW PacificM. cephalus were subdivided into lineages 1 and 2 as described in Jamandre et al. [16]. The CR sequence from the whole mtDNA genome (a sample from Japan) was incorporated and assigned to the lineage 2 of NW PacificM. cephalus. Hawaii was assigned as the central Pacific population. NW Atlantic and Gulf of MexicoM. cephalus were grouped together as no genetic divergence was found among these locations as described in Rocha-Olivares et al. [14, 15]. Togo and Senegal individuals were grouped as West Africa and South Africa has its own lineage.M. liza was assigned as a SW AtlanticM. cephalus group. Consensus sequence was obtained at each group, and these were eventually used for further phylogenetic and sequence divergence analyses.The most appropriate model for the CR data set was the GTR + G model with a gamma distribution shape parameterα
=
0.7997. Phylogenetic relationships ofM. cephalus revealed no geographic structuring. Unrooted Bayesian tree showed short branch lengths stemming from the ancestral nodes (Table 2, Figure 4). The star-like tree had nodes with moderate (71%) to highly significant (100%) posterior probability values. Gulf of Mexico and Atlantic samples gained the longest branch lengths implying high divergence from all of the other populations.M. liza individuals from French Guyana did not present a divergence level higher than the divergence observed withinM. cephalus sensu stricto.Figure 4
Transition (Ts) and transversion (Tv) over genetic distance (Tamura-Nei TrN-G) plot inferred from the whole CR ofMugil cephalus.


## 3.1. Control Region Sequence Structure and Variability

The observed lengths ofM. cephalus control region varied among sampling locations. The complete CR sequences ofM. cephalus from northwest Pacific ranged from 909 to 1015 bp, whereas in central Pacific (Hawaii), it ranged from 919 to 930 bp. The Atlantic and Gulf of MexicoM. cephalus CR sequences ranged from 884 to 894 bp while the West and South African individuals ranged from 899 to 912 bp and 911 bp, respectively. The CR length of theMugil liza individual is 906 bp.By comparingM. cephalus CR sequences to CR of other fishes species (see Table 2), three domains were highlighted, namely, the left and right hypervariable domain and the central conserved domain (Figure 1). Several conserved sequence blocks (CSBs) were located at the central conserved region and at the right hypervariable domain (Figure 2). CSB-I, -II, and -III are commonly found in vertebrate CR. These were also present in all sequences ofM. cephalus. Putative CSB-D was also found in the sequences (Figure 2). Pyrimidine tract was also identified at the central conserved domain.Table 2
Comparison of conserved sequence blocks (CSBs) identified in the control region ofMugil cephalus to other teleost fishes.
Species
CSB-D
CSB-I
CSB-II
CSB-III
Mugil cephalus
NW Pacific
Lineage 1/
TTCCTGGCATTTGGTTCC
AACATA
AAAACCCCCCCTACCCCCCCTAA
TCCTGAAAACCCCCCCGGAAACA
Lineage 2/
TTCCTGACATTTGGTTCT
AATATA
AAAACCCCCCTACCCCCCTAAA
TCCTGAAAACCCCCCGGAAACAA
Hawaii
TTCCTGGCATTTGATTCC
AACATA
AAACCCCCCCTACCCCCCCAAA
TCCTGAAAACCCCCCGGAAACAA
Atlantic
TCTTTGATGTTCAGGTTTT
AACATA
AAAACCCCCCCTACCCCCCCAAA
TCCTGAAAACCCCCCGGAAACAA
West Africa
TTTCTGGCATTTGGCTCC
AACATA
AAAACCCCCCCTACCCCCCCTAA
TCCTGAAAACCCCCCGGAAACAA
South Africa
TTCCTGGCATTTGGTTCC
AACATA
AAAACCCCCCCTACCCCCCCTAA
TCCTGAAAACCCCCCGGAAACAA
Mugil liza
TTCCTGGCATTTGGTTCC
AACATA
AAAACCCCCCCTACCCCCCCAAA
TCCTGAAAACCCCCCGGAAACAA
Chanos chanos
TTCAAGAGCATCTGGTTCC
—
AAACCCCCCCCTCCCCCCA
TGTTAAACCCCCAAAACCA
Crossostoma
lacustre
TTACTGGCATCTGGTTCC
CTATGTATGTACAATGAGCATAA
ACAAACCCCCCTACCCCCCT
TGCTCAAACCCCGAAACCA
Dascyllus trimaculatus
TTCCTGGCATTTGGTTCC
CACTACTGTCTTCCCGGACATACA
TAAAACCCCCCTACCCCCCT
TGAAACCCCCCGAAAACA
Notopterus notopterus
TATTCCTGGCATCTGGCTCC
—
AACCCCCCTACCCCCCA
TCAAACCCCTAAACCA
Perca fluviatilis
TTCCTGGCATTTGGTTCC
ATCTTAGGATATCAAGAGCATAA
TAAACCCCCCCTACCCCCCC
TGAAAACCCCCCCGGAAACA
Padogobius nigricans
TTTCTGGCATTTGGTTCC
ATAACTGATATCATGAGCATAA
AATCCCCCCTACCCCCCT
TGAAAACCCCCCGGAAACA
Sinipercine
spp.
TATTCCTGGCATTTGGTTCCT
ATTGGATATCAAGAGCATAA
TAAACCCCCCTACCCCCC
TGAAAACCCCCCGGAAACA
—: not identified.Multiple alignment of consensus sequences ofMugil cephalus from L1 (Lineage 1) and L2 (Lineage 2) of the NW (northwest) Pacific, Hawaii/C (central) Pacific, Gulf of Mexico/ Atlantic, Togo and Senegal/W (west) Africa, and South Africa. Structural features (CSB—conserved sequence block and Pyrimidine-tract) are also indicated.
(a)
(b)Insertions and deletions (indels) were found in both the left and right hypervariable domain and generally associated with repetitive sequences (Figure2). The percent composition of nucleotide varied from 35.96 to 37.50% (A), 29.80 to 34.57% (T), 16.88 to 20.64% (C), and 11.29 to 13.27% (G), which infers thatM. cephalus CR is AT rich and poor in G.The transition/transversion (Ts/Tv) rate ratios arek
1
=
20.443 (purines) and k
2
=
14.496 (pyrimidines). An overall Ts/Tv ratio of 5.131 was estimated from the dataset (Figure 3). A linear relationship between the number of nucleotide substitutions per site and the sequence divergence was observed for both, transversion and transitions mutations (Figure 3). Slope of the transitions /divergence regression indicates a higher transition rate than transversion rate.Figure 3
Bayesian tree inferred from complete CR sequences. Posterior probability values (%) are indicated close to the nodes of the branches.


## 3.2. Phylogenetic Relationship and Sequence Divergence

NW PacificM. cephalus were subdivided into lineages 1 and 2 as described in Jamandre et al. [16]. The CR sequence from the whole mtDNA genome (a sample from Japan) was incorporated and assigned to the lineage 2 of NW PacificM. cephalus. Hawaii was assigned as the central Pacific population. NW Atlantic and Gulf of MexicoM. cephalus were grouped together as no genetic divergence was found among these locations as described in Rocha-Olivares et al. [14, 15]. Togo and Senegal individuals were grouped as West Africa and South Africa has its own lineage.M. liza was assigned as a SW AtlanticM. cephalus group. Consensus sequence was obtained at each group, and these were eventually used for further phylogenetic and sequence divergence analyses.The most appropriate model for the CR data set was the GTR + G model with a gamma distribution shape parameterα
=
0.7997. Phylogenetic relationships ofM. cephalus revealed no geographic structuring. Unrooted Bayesian tree showed short branch lengths stemming from the ancestral nodes (Table 2, Figure 4). The star-like tree had nodes with moderate (71%) to highly significant (100%) posterior probability values. Gulf of Mexico and Atlantic samples gained the longest branch lengths implying high divergence from all of the other populations.M. liza individuals from French Guyana did not present a divergence level higher than the divergence observed withinM. cephalus sensu stricto.Figure 4
Transition (Ts) and transversion (Tv) over genetic distance (Tamura-Nei TrN-G) plot inferred from the whole CR ofMugil cephalus.


## 4. Discussion

We compiled and analyzed new sequences and extracted from GenBank database of NCBI all available complete CR sequences ofM. cephalus from different populations around the world. The total length of the CR varies greatly among population/lineage. The longest CR sequences were found in individuals belonging to the NW Pacific (909–1015 bp), whereas the shortest were observed in individual samples in the Atlantic and the Gulf of Mexico (884–894 bp). Such size variation is consistent with observation from other teleost fish (e.g., Cichlid—888 bp; Atlantic cod—996 bp; sturgeon—976; salmon—1071; yellowtail ~ 1600; from (2)). However, to date, there are no other reports of fish species with such highly intraspecific sequence length difference. The observed CR length variations inM. cephalus are either due to the presence of variable number of tandem repeats (VNTR) at the right hypervariable domain (Figure 2) and indels. Indels inM. cephalus were already described by Rocha-Olivares et al. [14]; however they only gave general descriptions on indels patterns and phylogenetic significance. In the sequence alignment presented here, there are two sets of permanent deletions unique to lineage 1 of northwest Pacific samples: 22 bp deletion (starting at position 857 bp in the alignment; Figure 2) and 7 bp deletion 45 bp following the former. However, the lineage 1 of northwest Pacific has tandem repeats (ATn) located at 30 bp after CSB-I, which cannot be found in the other populations. Second, West African samples were found to have 10 to 14 bp deletions at the central conserved domain/region.Lastly, both the left and right hypervariable domain contain considerable amount of indels that greatly affected the sequence length of CR among populations ofM. cephalus. An 88 bp indel due to repeated sequence (TATAATATn) was found again in the 3′ or left hypervariable domain in lineage 1 of NW Pacific when the consensus sequences among populations were aligned (Figure 2). These indels seem to be unique to lineage 1 of NW Pacific individuals. Indels in CR are common in teleost fish [35]. However large deletions are uncommon and only reported in some fish: an 86 bp indel in the 3′ region ofOncorhynchus mykiss [36], 40 bp indel in central conserved region of the rockfishSebastes capensis [37], and 48 bp indel in milkfishChanos chanos [26]. These studies focused more on the interspecies sequence variation with the exception of the milkfish. The large deletions found in this study are intraspecific comparison and seem to be very rare in fish. This is the first study to show large deletions in CR within theM. cephalus species and the family Mugilidae. These indels were hypothesized not to affect the organism due to the noncoding nature of CR [1]; however these indels might be of importance for population identification. The CSBs observed inM. cephalus are commonly found in other teleosts with slight length variations. Comparing the CSB-I to other fish (Table 2), it is obvious that the CSB-I inM. cephalus is shorter than the other fishes. CSB-I appears inM. cephalus to be more conserved than CSB-II and CSB-III. The two latter CSBs show a few nucleotide indels in the alignment.The phylogenetic tree inferred from CR shows thatM. cephalus may have experienced quick radiation due to short branch lengths stemming from the ancestral node. Such branching pattern can indicate either homoplasic mutations that can limit resolution of ancient nodes or a quick ancestral radiation. Rapid evolutionary rate as reported in CR of Chaetodontidae [10] may result in an excess of homoplasic mutations ~ mutations that reached a “ceiling.” Because no transition or transversion saturation was detected inM. cephalus CR sequences, this assumption was not considered. However, comparing between samples from different locations, the long branch lengths can be attributed to deep divergences between these samples. Livi et al. [38] described that these divergences observed among worldwide samples ofM. cephalus might have been due to both demographic (bottleneck events, etc.) and geographic and temporal isolations.The star-like shaped phylogenetic tree was already observed in cichlids of the African lakes, which was interpreted as a rapid radiation [39–41]. Such an evolutionary event can be attributed to both allopatric and sympatric speciation, which is due to geological events, reproductive isolation, and other selections caused by intrinsic and extrinsic forces [41]. In the case of a South American silversideOdontesthes perugiae species complex (a strictly freshwater species), rapid evolution was assumed to be due to its marine-estuarine origin and geological events [38, 42].Lastly, the star shaped pattern of the phylogenetic tree inferred no geographic patterning according to the sequence analysis. Previous phylogenetic studies using other regions of mtDNA have also shown no geographic pattern among worldwide population ofM. cephalusgenetic structure [18, 19, 21]. In conformity with other studies and the comprehensive Mugilidae molecular phylogeny by Durand et al. [22], we clearly confirmed the presence ofM. liza in theM. cephalus clade because there are no other or any long branches present to differentiate this species from otherM. cephalus haplotypes. The high CR sequence divergence leads Rocha-Olivares et al. [14, 15] to reconsider the real taxonomical status ofM. cephalus which was also questioned by Heras et al. [21]. This indicates the presence ofM. liza in theM. cephalus clade.


## 5. Conclusion

This study presents a very high sequence and genetic divergence within groups of the same species. Applying the definition of the biological species concept is hard for globally distributed species, mainly due to difficulties in the determination of interbreeding levels between populations or species. In congruence to the growing evidence and as indicated in this study, on a worldwide scale,M. cephalus should be regarded as a complex species.



---

*Source: 564105-2014-04-02.xml*