Three genetic stocks of frigate tuna Auxis thazard thazard (Lacepede, 1800) along the Indian coast revealed from sequence analyses of mitochondrial DNA D-loop region

Abstract Frigate tuna Auxis thazard thazard is a cosmopolitan species and one of the smallest members of the tribe Thunnini (the true tunas), and currently managed as a single stock. In the present study, genetic variation was surveyed using sequence data of mitochondrial DNA D-loop region to test for the presence of genetic stock structure of frigate tuna along the Indian coast. A total of 364 individuals were sampled from 8 major fishing zones along the Indian coast. Significant genetic heterogeneity was observed for the sequence data (ΦST = 0.0439, P = < 0.001). Analysis of molecular variance (AMOVA) showed significant genetic variation among the three groups analysed (ΦCT = 0.1223, P = < 0.05), which was also supported by spatial AMOVA results. Therefore, the null hypothesis of single panmictic population of frigate tuna along the Indian coast can be rejected. Phylogenetic analysis of nucleotide sequences demonstrated that frigate tuna can be grouped into three different mitochondrial clades (Clades I, II and III). However, there were no significant genealogical branches or clusters of samples corresponding to sampling locality. The results of the present study suggest the possibility of three genetically differentiated units of frigate tuna across the coastal waters of India.


Introduction
Tunas are large, highly migratory pelagic species. Genetic differentiation has generally been supposed to be low among tuna populations within and between the oceans. This is because of high dispersal potential and the apparent lack of migrational barriers at sea. Nevertheless, recent advances in molecular genetic methods have allowed the analysis of large sample sets with higher sensitivity, consequently detecting small to large signals of population subdivision in several species.
Frigate tuna Auxis thazard thazard (Lacepede, 1800) is a coastal tuna species found circumglobally in tropical oceans up to depths of 50 m (Collette & Nauen 1983). The species is mainly confined to continental shelves and has a localized migratory habit (Maguire et al. 2006). From larval records, it is deduced that frigate tuna spawn throughout its distribution range. In correlation with temperature and other environmental changes, spawning season varies with areas, but in some places it may even extend throughout the year. In the southern Indian Ocean, the spawning season extends from August to April; north of the equator it is reported from January to April (Collette & Nauen 1983), while in Indian waters spawning peak is August to November (Pillai & Gopakumar 2003). One of the important features of Auxis spp. is that larvae have the widest range of temperature tolerance (between 21.6 and 30.5˚C) among tuna species studied. The maximum population doubling time is less than 15 months, while maximum reported age is 5 years. It feeds on small fish, squids, planktonic crustaceans (megalops) and stomatopod larvae. Because of their abundance, they are considered an important element of the food web, particularly as prey for other species of commercial interest (Collette & Nauen 1983).
Frigate tuna has been mainly exploited for canned products due to the excellent properties of the meat, with its mild taste and low cholesterol content (Infante et al. 2004). Despite abundance and high commercial importance, there has been no investigation on the genetic structure of frigate tuna. Earlier genetic studies of frigate tuna were mainly confined to species identification using cytochrome b and cytochrome c oxidase I sequence data analysis (Robertson et al. 2007;Melissa et al. 2008); authentication of canned products using DNA markers (Infante et al. 2004;Lin & Hwang 2007) and complete sequencing of mitochondrial DNA (Catanese et al. 2007). In India, the first attempt at compiling and consolidating scientific information on tuna was made at the symposium on scombrids, in 1962 at Mandapam (Sudarsan 1993). Assessments of stock and exploitation rates based on length frequency data were made by Silas et al. (1986) and James et al. (1987). Based on the population parameters and stock estimates, the coastal tunas have been found to be exploited at or above the optimum levels. Along the coastal waters of India the exploitation rate of frigate tuna is 0.72, whereas the optimum exploitation rates was at 0.40, indicating that frigate tuna is exploited relatively at high rates in the coastal waters (Pillai & Gopakumar 2003). Commercial tuna fishing in India started in 1985 with the operation of chartered tuna long-liners and few Indian-owned tuna vessels. As more tuna vessels continue to be introduced, tuna is likely to face problems of over fishing (Menezes et al. 2006). Because many species of pelagic fishes are under intense fishing pressure, it is imperative that their population genetic structure be assessed in order to provide proper management at both regional and international levels. Genetic markers in fisheries biology are used to determine whether samples from natural populations are genetically differentiated from each other. The detection of differentiation would imply that source groups comprise different stocks and should be treated as separate management units.
Mitochondrial DNA (mtDNA) has proven to be a useful genetic marker in population genetic studies and fisheries management due to its high mutation rate, haploid nature and maternal mode of inheritance, which makes it a sensitive indicator of genetic drift resulting from geographical subdivision (Lindak & Paul 1994;Garber et al. 2005). Sequence analyses involving the rapidly mutating mtDNA control region has been used to identify genetic stocks in a growing list of marine fish species (e.g. Pardini et al. 2001;Stepien et al. 2001;Guarniero et al. 2002;Tringali et al. 2003;Garber et al. 2004;Menezes et al. 2012), including members of Scombroidei (Alvarado Bremer et al. 1998;Ely et al. 2002).
In this study, sequence analyses of the mitochondrial D-loop region were employed to examine the levels of genetic diversity among frigate tuna samples collected from the Indian coast and to identify discrete genetic units, if they exist, for management, while testing the null hypothesis of panmixia. Selection of a 500-bp mtDNA D-loop fragment for analysis was based on its documented hypervariability and the availability of specific primers for amplification of this region (Menezes et al. 2006(Menezes et al. , 2012. This is the first study to provide evidence of genetic differentiation in frigate tuna.

DNA isolation
A total of 364 fin clip samples of frigate tuna were collected from 8 locations off the Indian coast (Table I; Figure 1) and preserved in absolute alcohol until DNA extraction. High molecular weight genomic DNA was extracted from the fin clip using the standard TNES-urea-phenol-chloroform protocol (Asahida et al. 1996). The DNA samples were stored at 4 C prior to PCR analysis.

DNA amplification
A fragment of 500 bp containing the first half of the control region (D-loop) was amplified using the primer set designed (Menezes et al. 2006) from a GenBank sequence of Auxis thazard (accession number NC005318). The primer sequences were as follows: 5? CCGGACGTCGGAGGTTAAAAT 3? (forward) and 5? AGGAACCAAATGCCAG-GAATA 3? (reverse). Amplification was carried out by the procedure as described by Menezes et al. (2006). The verification of successful polymerase chain reaction (PCR) amplification was assessed by 1% agarose gel electrophoresis.

DNA sequencing
Amplified PCR products were purified using Axygen PCR Cleanup Kits (Axygen Biosciences, California, USA) as recommended by the manufacturer. Initial sequencing reactions were performed in both directions with the two primers used to amplify the mtDNA D-loop region. As the forward and reverse sequences were perfectly complementary, remaining samples were sequenced only in the forward direction. Sequencing reactions were prepared with the BigDye Terminator Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) according to the supplier's recommendations. Sequencing reactions were resolved on ABI 3130xl Genetic Analyzer (Applied Biosystems, Foster City, CA, USA). Representative sequences have been deposited in GenBank, with accession numbers JN398671-JN399010.

Data analyses
DNA sequences were edited with the program BioEdit (version 7.0.1, Hall 1999) and aligned using the ClustalW algorithm (Thompson et al. 1994) as implemented in MEGA5 (Tamura et al. 2011). DnaSP 4.0 (Rozas et al. 2003) was used to calculate the number of haplotypes, polymorphic sites, nucleotide diversity (p) and haplotype diversity (h) from sequence data. Molecular diversity indices, such as transitions, transversions, substitutions, and indels were obtained using program Arlequin version 3.11 (Excoffier et al. 2005). The aligned sequences were used to analyse the population structure and genetic variation by using Arlequin version 3.11 (Excoffier et al. 2005). An analysis of molecular variance (AMOVA) (Excoffier et al. 1992) was used to test the spatial genetic variation among sampling areas by estimating F ST . Hierarchical AMOVA was used to examine the amount of genetic variability partitioned among groups, among populations within groups, and within populations. The significance of F ST was tested by 1000 permutations for each pairwise comparison. For all F ST analyses, the Tamura & Nei (1993) distance method was used.
Spatial structure was investigated via spatial analysis of molecular variance implemented in the program SAMOVA v. 1.0 (Dupanloup et al. 2002) to identify groups of geographically adjacent popu- lations of frigate tuna that were maximally differentiated based on sequence data. This analysis identifies groups of sample sites that are most similar and that are geographically meaningful. The method is based on a simulated annealing procedure that maximizes the F CT between population groups. Our analyses were based on 100 simulated annealing steps, and a prior definition of the number of groups, K, ranging from 2 to 7. For each analysis, with increasing K, we examined the proportion of genetic variance due to differences between groups (F CT ) to find the value of K for which F CT was largest and statistically significant.
The null hypothesis of neutral evolution was tested using Tajima's D test (Tajima 1989) and Fu's Fs test (Fu 1997) with 1000 permutations as implemented in Arlequin. Tajima's D test compares two estimators of the mutation parameter θ: Watterson's estimator θs and Tajima's estimator θp. A significant D value can be due to factors such as population expansion, bottlenecks, or heterogeneity of mutation rates (Tajima 1989). Fu's Fs test compares the number of haplotypes observed with number of haplotypes expected in a random sample under the assumption of an infinite sites model without recombination (Fu 1997). Fu (1997) found that Fs is sensitive to population demographic expansion, which generally leads to large negative Fs values.
Arlequin was used to calculate the historic demographic parameters θ 0 (population before expansion), θ 1 (population after expansion) and s (relative time since population expansion). The θ (θ 0 and θ 1 ) are product of 2mN 0 and 2mN 1 (where m is mutation rate and N is the effective population size at times 0 and 1). The tau (s) value can be transformed to estimate the actual time (T) since population expansion using formula T = s/2m where m is the mutation rate per site per generation. In the present study, the mutation rate of 3.6Â10 8 mutations per site per year was applied for the control region sequence of frigate tuna as this rate has been reported for the mtDNA control region in teleosts (Donaldson & Wilson 1999). The goodness-of-fit of observed mismatch distributions to the theoretical distribution under a sudden expansion model (Rogers & Harpending 1992) was tested with the sum of squared deviations (SSD) between observed and expected mismatch distributions, and with the raggedness index (Hri) of Harpending (1994), as implemented in Arlequin.
Estimates of the expected number of female migrants between sampling localities per generation (N f m) were calculated by the program Migrate-n version 3.2.16 (Beerli & Palczewski 2010) using the maximum-likelihood approach, employing 10 short chains and one long chain, and a burn-in period of 10,000 generations for each chain. The bestfitting model of substitution was determined with jModeltest 0.1.1 (Posada 2008) using the Akaike information criterion (AIC). The TIM2ef + G model with gamma = 0.6240 (lnL = 7557.1823; AIC = 16478.3646) was chosen as the best-fitting model. Phylogenetic relationships among mitochondrial D-loop haplotypes were assessed using the neighbor-joining (NJ) (Saitou & Nei 1987) method as implemented in MEGA5 (Tamura et al. 2011). The tree was rooted with a control region sequence (GenBank: JF752303.1) from Katsuwonus pelamis as an outgroup. Only nodes with bootstrap support of greater than 50% are shown in the final tree. All positions containing alignment gaps and missing data were eliminated by invoking the pairwisedeletion option for indels. Evaluation of statistical confidence in nodes was tested by 1000 bootstrap replicates (Felsenstein 1985).

Genetic diversity
A fragment of 500 bp containing the first half of the control region (D-loop) was sequenced in 364 frigate tuna. A total of 178 variable sites (131 parsimony informative), constituting 340 haplotypes were detected among mtDNA control region sequences. Frequencies and occurrences of each haplotype are presented as an online supplementary table. Most of the haplotypes were unique to particular individuals. The nucleotide composition was consistent with A-T bias reported in D-loop region of mitochondrial genome (Brown et al. 1986). Mean base composition of D-loop was 36.96% A, 32.91% T, 14.80% C, 15.34% G and observed transitions outnumbered observed transversions by a mean ratio of 1.94. The estimate of haplotype diversity (h) was high for all the eight populations with values ranging from 0.989 to 1.000 (Table I). The nucleotide diversity (p) for sequence data varied greatly among the samples with values ranging from 0.026 to 0.091 (Table I).

Population genetic structure
Analysis of molecular variance (AMOVA) performed on mtDNA sequence data set revealed significant genetic heterogeneity among the eight sampling sites (F ST = 0.0439, P = < 0.001) (Table II). In addition, hierarchical analysis of molecular variance showed significant differentiation among the three groups analysed (F CT = 0.1223, P = < 0.05) (Table II). Eleven pairwise comparisons of sequence data among localities were significant. However, after sequential Bonferroni correction (Rice 1989 (Table III).
SAMOVA was performed to test the significance of the partitioning of genetic variance resulting from different groupings of the sampling localities into geographical groups. For mtDNA, an assumption of three groups (K = 3) was geographically meaningful group with statistically significant F CT value (F CT =0.10154, P =< 0.05) ( Table IV).
Estimates of female migrants among three groups (VE, RA, KA, TU, PO, VI); (KO); (PB) varied from 63.93 to 215.85. The highest levels of inferred migration were from the (PB) to (VE, RA, KA, TU, PO, VI), while the lowest were observed between (KO) and (PB).
Phylogenetic analysis of nucleotide sequences demonstrated that D-loop sequences of frigate tuna were grouped into three different mitochondrial clades (Clades I, II and III) (Figure 2). Haplotypes belonging to Clade I were found in all sampling sites, while haplotypes of Clades II and III were restricted only to two (PO and VI) and four (RA, KO, PB and PO) locations, respectively. Clade I contained 326 of the 340 haplotypes identified, while only five and nine haplotypes were present in Clades II and III. Clade I is weakly supported by a bootstrap value less than 50%. In contrast, Clades II and III are strongly supported with 98% and 100% bootstrap values, respectively. However, haplotypes from individual sites or geographic regions do not form exclusive clades, and sharing of haplotypes occurs across sites and regions.

Historic demography
The observed mismatch distribution for all the samples was slightly multimodal (Figure 3).
Tajima's D test of selective neutrality and Fu's Fs test was highly negative and significant (Table II), suggesting a sudden population expansion, which is supported by the non-significant sum of squared deviation (P = 0.3873) and Harpending's raggedness index (P = 0.6453). Rapid population expansion was also supported by large differences in θ 0 and θ 1 within all collections localities ( Table V). The overall tau (s) value of sequence data was estimated to be 11.992. Following the equation T = s/2m and mutation rate 3.6Â10 8 per site per year (Donaldson & Wilson 1999), it was estimated that expansion occurred approximately 105,027 years ago for all the samples, with the exception of PB which occurred about 668,277 years ago.

Population genetic structure
Hierarchical analysis of molecular variance of mtDNA D-loop sequences showed significant genetic heterogeneity among the three groups (VE, RA, KA, TU, PO, VI), (KO), and (PB) analysed. Results were further corroborated by SAMOVA analysis, which recognized KO and PB as discrete populations, and all remaining sites (VE, RA, KA, TU, PO, VI) together as a third discrete population. Thus, the null hypothesis of single panmictic population of frigate tuna along the Indian coast can be Large marine pelagic fishes may have a homogeneous genetic population structure on a large geographical scale, as they are highly migratory with a very wide reproductive area and have extensive egg and larval dispersal through/by ocean currents. However, coastal pelagic fishes generally have more geographic population subdivisions than do oceanic species (Crosetti et al. 1994;Rossi et al. 1998), probably because of their inshore habitat requirements, shorter migration distances, and vulnerability to climatic fluctuations. Frigate tuna being a coastal species conforms to this pattern in the present study. Prior to the present study, there have been no mtDNA genetic analyses on frigate tuna population structure. However, a recent study undertaken on skipjack tuna (Katsuwonus pelamis) from the Indian coast using sequence analysis of the Dloop region indicated four different stocks (Menezes et al. 2012). The possibility of genetically discrete yellowfin tuna populations has been observed in the northwestern Indian Ocean using mtDNA sequence and microsatellite data (Dammannagoda et al. 2008). Genetic analysis of an mtDNA gene and six microsatellite loci revealed two stocks of skipjack tuna in northwestern Indian Ocean (Dammannagoda et al. 2011). Intra-oceanic population structuring has also been observed in several other tuna and scombroid fishes (Sharp 1978;Ward et al. 1997;Reeb et al. 2000;Chow et al. 2000;Alvarado Bremer et al. 2005). Sequence analysis of mtDNA of bigeye tuna suggests the existence of two mitochondrial lineages in Pacific Ocean (Chiang et al. 2006).
Environmental factors, both past and present, play a large role in organizing population structure by creating or eliminating connections among populations. Historical environmental events relevant to the creation and maintenance of distinct fish popula-tions include glaciations, formation of land bridges, and sea level changes. Present geographical barriers include thermoclines, nutrient-poor currents or ocean gyres (Sinclair & Iles 1989). The Andaman Sea is a part of the Indian Ocean and is located to the southeast of the Bay of Bengal, west of Thailand, south of Myanmar, and east of the Andaman and Nicobar Islands. The Andaman Sea, encircled by islands and land areas as stated above, has been a distinct closed and isolated basin during the Pleistocene, the last glacial period (about 12,000 years ago), when the sea level was estimated as being 100 m lower than at present (Morley & Flenley 1987). Several fish species are considered to have speciated in the closed basin in the Pleistocene and become endemic to the Andaman Sea at present (Randall & Satapoomin 1999;Motomura & Senou 2002). The closed and isolated basin of the Andaman Sea may have constituted a barrier to the gene flow between Port-Blair and rest of the Indian coast frigate tuna populations, and the two populations may only rarely reinforce one another with recruits.
Besides environmental factors, many life-history traits such as spawning play a significant role in determining genetic variability and population structure of marine fish species. The observed population structure of KO could be related to specific spatial spawning regimes. George (1989) assessed the distribution and abundance of Auxis spp. eggs and larvae along the southwest coast of India. The highest density of larvae was located between 10å nd 11˚30?N, indicating that southwest coast of India, including the waters of Kochi, are a spawning site of frigate tuna. Thus, we suggest that the divergence observed in the present study may result from fidelity to spawning sites of individuals from these areas.
Phylogenetic analysis of sequence data revealed three different mitochondrial clades (Clades I, II and III). However, these clades did not have genealogical branches or clusters of samples corresponding to sampling locality. In several other tuna and scombroid fishes, including bigeye tuna (Alvarado Bremer et al. 1998;Chow et al. 2000), blue marlin (Finnerty & Block 1992), sailfish (Alvarado Bremer 1994; Table IV. Population structure based on mtDNA differentiation of frigate tuna samples (in spatial analysis of molecular variance, SAMOVA). The row in bold type indicates the details of geographically meaningful groups with maximum genetic differentiation. The tree is rooted with a control region sequence from Katsuwonus pelamis. Graves & McDowell 1995) and swordfish (Alvarado Bremer et al. 1995Rosel & Block 1996), the heterogeneous distribution of divergent sets of mtDNA lineages has been reported. The asymmetrical pattern of distribution of these clades within a species is due to secondary contact and interbreeding between populations subsequent to their geographical isolation for a long period (Durand et al. 2005). Another explanation for physical mixture of the frigate tuna clades at sampling sites includes the high rate of migration and differential passive transport of larvae belonging to the different clades to fishing grounds by monsoonal currents in the Indian Ocean.

Genetic diversity
The high levels of genetic diversity at both haplotypes and nucleotide levels were observed for frigate tuna samples along the Indian coast. High levels of divergence between haplotypes and nucleotides are associated with a long evolutionary history in a large, stable population, or with secondary contact be-tween different lineages (Grant & Bowen 1998). Nucleotide and haplotypes diversity are generally higher in older populations. Thus, the samples from Port-Blair region may represent older populations, with higher values of haplotype and nucleotide diversity. This statement is also supported by the highest tau value and lowest ti/tv being ratio obtained for this population.

Historical population dynamics
The tau (s) values for the mtDNA sequence data set were all in same range except for PB population, indicating that all the samples (except PB) originated from a single colonization event about 105,027 years before present, while the population of Port-Blair originated from a separate colonization event about 668,277 years before present. The long evolutionary period for the Port-Blair samples may have allowed its divergence as a genetically different population from rest of the samples of the present study.
Overall, large differences were observed between θ 0 and θ 1 , suggesting a rapid population expansion event. Population expansion was also supported by overall negative and significant values of Tajima's D test and Fu's Fs test. In addition, all Harpending's raggedness indices and mismatch distributions were nonsignificant, which should be expected when populations have undergone population expansion. However, the graph of the observed mismatch distributions did not strictly follow the simulated curves under the sudden expansion model, but were slightly multimodal, which underlines the finding of substructures in the data sets (Schneider & Excoffier 1999;Ray et al. 2003). Large differences between the initial θ 0 value of Port-Blair and that of the remaining populations suggest a significantly larger effective population size of frigate tuna in the Andaman Sea. In addition, the neutrality test for Port-Blair samples is not significant, suggesting that the effective size of the frigate population has been large  *, **, *** Significant at P < 0.05, P < 0.01, and P < 0.001, respectively. and stable for a long period. This further strengthens the suggestion of genetically independent populations of frigate tuna in the Andaman Sea.

Conclusion
In conclusion, this is the first study to provide evidence of genetic differentiation of frigate tuna in Indian waters. The results of the present study demonstrate three genetically distinct populations of frigate tuna in Indian waters; one around the Andaman Sea, a second around the southwest coast of India and a third around the rest of the Indian Seas. However, the sample size for Port-Blair was small (N=14) and, because of maternal mode of inheritance, mtDNA variability offers only a partial view of frigate tuna biology and population histories. Therefore, further studies with complementary genetic markers (i.e. microsatellite DNA), additional sampling, and other biological studies, including tagging experiments, are needed to corroborate the results presented here. However, based on these results, we would propose that the Indian frigate tuna stock be managed as at least three distinct units: one in Kochi, the second in the Andaman Sea and the third along the rest of the Indian coast. Failure to manage distinct populations may lead to local overexploitation, extreme reduction of genetic variability (i.e. genetic bottleneck), and ultimately to severe depletions due to loss of fitness and the viability to withstand survival challenges. Given the fact that no conservation strategies have been implemented for this species, the need to establish independent exploitation regimes for each of the population units identified should be evaluated.