Abstract

Accurate identification of units for conservation is particularly challenging for marine species as obvious barriers to gene flow are generally lacking. Bryde’s whales (Balaenoptera spp.) are subject to multiple human-mediated stressors, including fisheries bycatch, ship strikes, and scientific whaling by Japan. For effective management, a clear understanding of how populations of each Bryde’s whale species/subspecies are genetically structured across their range is required. We conducted a population-level analysis of mtDNA control region sequences with 56 new samples from Oman, Maldives, and Bangladesh, plus published sequences from off Java and the Northwest Pacific. Nine diagnostic characters in the mitochondrial control region and a maximum parsimony phylogenetic analysis identified 2 genetically recognized subspecies of Bryde’s whale: the larger, offshore form, Balaenoptera edeni brydei, and the smaller, coastal form, Balaenoptera edeni edeni. Genetic diversity and differentiation indices, combined with a reconstructed maximum parsimony haplotype network, indicate strong differences in the genetic diversity and population structure within each subspecies. Discrete population units are identified for B. e. brydei in the Maldives, Java, and the Northwest Pacific and for B. e. edeni between the Northern Indian Ocean (Oman and Bangladesh) and the coastal waters of Japan.

Barriers to gene flow for cetaceans are rarely evident in marine environments (Mendez et al. 2010) meaning that discrimination of lower level conservation units is challenging, as geographic distribution is not an appropriate proxy for isolation. Genetics can be a powerful tool for discriminating among incipient species and geographic forms, as well as distinct demographically independent populations that are experiencing levels of gene flow too high for local adaptation to occur (Taylor 2005). Notable examples of population-level delineations in baleen whales using genetics include humpback whales (Megaptera novaeangliae) in the North Pacific (Baker et al. 1998), North Atlantic (Stevick et al. 2006), Arabian Sea, and South Atlantic and Indian Oceans (Rosenbaum et al. 2009), and blue whales (Balaenoptera musculus) in the Southern Hemisphere (LeDuc et al. 2007).

Despite these advances, the taxonomy and population structure of many cetaceans remain unresolved. The implications of the existence of undetected conservation units at species and distinct population levels are disquieting, especially for taxonomic groups hunted under scientific permit from the International Whaling Commission or those recovering from commercial whaling (Clapham et al. 2008). There is also the potential for the specialized habitat requirements of distinct lineages to be obscured by being aggregated within larger taxonomic groups. This issue is particularly consequential when lower level conservation units inhabit areas that can be potentially affected by human activities, such as fisheries interactions and hydrocarbon exploration and development.

Well over a century has passed since the Bryde’s whale (Balaenoptera edeni) was first described, but the phylogeny of this species complex is still unresolved (Perrin and Brownell 2007). Although the nomenclature is unsettled because the species genetics of the holotype specimen of B. edeni has not yet been determined, 2 subspecies are provisionally recognized by their genetics: a larger pelagic form, Balaenoptera edeni brydei, with a circumglobal distribution in tropical and subtropical waters of the Pacific, Atlantic, and Indian Oceans, and a smaller nearshore form, Balaenoptera edeni edeni, in the Indo-Pacific region (Committee on Taxonomy 2011). However, others have recognized 2 species rather than subspecies (B. brydei and B. edeni; Wada et al. 2003; Sasaki et al. 2006; Kanda et al. 2007; Kato and Perrin 2009). For the purposes of maintaining consistency with current nomenclature (Committee on Taxonomy 2011), we refer to the subspecies B. e. brydei and B. e. edeni, or “large form” and “small form.”

A single-species designation of Bryde’s whales was broadly accepted until the 1990s. Recently, however, it was discovered that populations in several parts of the range exhibit differences in body size, including a larger offshore form (i.e., B. e. brydei) and 1 or more smaller, predominantly coastal forms (i.e., B. e. edeni; Perrin et al. 1996; Best 1997, 2001; Perrin and Brownell 2007; Penry et al. 2011). A new species, Balaenoptera omurai, representing a separate ancient lineage within the Balaenopteridae clade (Sasaki et al. 2006), was also recently described in the Indo-Pacific region (Wada et al. 2003). As Bryde’s whales were previously subjected to commercial exploitation and remain a target of scientific whaling by Japan (Kanda et al. 2007), the ability to distinguish Bryde’s whale taxa and elucidate their respective genetic population structure is required to avoid overexploitation, develop effective conservation plans, and prevent the loss of irreplaceable evolutionary lineages.

Here, we build upon previous research by combining new genetic samples of Bryde’s whales from 3 previously unsampled locations across the Northern Indian Ocean (NIO) with previously published data on samples from the Central Indo-Pacific region and Northwest Pacific Ocean (Yoshida and Kato 1999; Kanda et al. 2007). Through the integration of these data sets, we provide additional insights into the Bryde’s whale phylogeny that supports the existing classification of the 2 taxonomic units (here treated as subspecies): B. e. brydei (large form) and B. e. edeni (small form). We then make population-level inferences across the region, which provide an important baseline for understanding the genetic diversity and spatial structure of Bryde’s whale populations; information that is vital for effective conservation and management.

Materials and Methods

Samples and Molecular Methods

A total of 409 samples originating from across the Western and Central Indo-Pacific and the Northwest Pacific Ocean were used for this study, including those previously published (Yoshida and Kato 1999; Kanda et al. 2007). The study region is defined following the Marine Ecoregions of the World (MEOW) schema developed by Spalding et al. (2007) and encompasses the Western Indo-Pacific Realm eastwards from the Somali/Arabian Province, the Central Indo-Pacific Realm, and the Warm Temperate Northwest Pacific Province nested within the Temperate North Pacific Realm (Figure 1). Our dataset includes 56 newly collected samples from Bangladesh, the Maldives, and Oman (see Supplementary Material online for details). Thirty samples were from biopsies of whales in Bangladesh (BAN). Of these, 29 were sampled from the rim of the Swatch-of-No-Ground (SoNG) submarine canyon and 1 originated from a stranding at Cox’s Bazaar in southeast Bangladesh. Previously unpublished data for the mtDNA control region were obtained for 8 whales sampled off the Maldives (MAL) and 18 individuals stranded or struck by ships along the coast of Oman (OMA). These new genetic data were combined with mitochondrial haplotypes from the south of Java (JAV, n = 27), the coastal waters of Japan (COJ, n = 16), and with a large dataset from the Northwest Pacific (NWP, n = 310; ACCN: EF068013-048, EF068060-063, Kanda et al. 2007; ACCN: AF146378-388, Yoshida and Kato 1999).

Figure 1.

Schematic illustrating the extent of the study region and approximate sampling locations shaded in gray: OMA, Oman; MAL, Maldives; BAN, Bangladesh; JAV, south of Java; COJ, coast of Japan; NWP, Northwest Pacific. New samples were collected from Oman, the Maldives, and Bangladesh. Existing samples had been previously collected from south of Java, coast of Japan, and Northwest Pacific (Yoshida and Kato 1999; Kanda et al. 2007). The eastern portion of the schematic (JAV, COJ, and NWP) was adapted from Yoshida and Kato (1999) and Kanda et al. (2007).

Figure 1.

Schematic illustrating the extent of the study region and approximate sampling locations shaded in gray: OMA, Oman; MAL, Maldives; BAN, Bangladesh; JAV, south of Java; COJ, coast of Japan; NWP, Northwest Pacific. New samples were collected from Oman, the Maldives, and Bangladesh. Existing samples had been previously collected from south of Java, coast of Japan, and Northwest Pacific (Yoshida and Kato 1999; Kanda et al. 2007). The eastern portion of the schematic (JAV, COJ, and NWP) was adapted from Yoshida and Kato (1999) and Kanda et al. (2007).

Total genomic DNA was extracted following procedures outlined in the QIAamp Tissue Kit (QiaGen, Venlo, Netherlands). A 407bp consensus fragment of the mtDNA control region was amplified using polymerase chain reaction (PCR), with primers Dlp 1.5 and Dlp 5 (Baker et al. 1993). Reactions in 25 µL total volume, containing 21.0 µL H20, 1.0 µL of each primer at 10 µM concentration, 1 Illustra (tm) PuReTaq (tm) Ready-To-Go (tm) PCR Bead (GE Healthcare, Pittsburgh, Pennsylvania), and 2.0 µL DNA template, were conducted using an Eppendorf Gradient Mastercycler (94 °C for 4min, followed by 30 cycles of 94 °C for 45 s, 54 °C for 45 s, and 72 °C for 45 s, and a final extension step at 72 °C for 10min). Amplified PCR products were purified using Agencourt AMPure XP (Agencourt Bioscience, Beverly, Massachusetts) and sequenced with dye-labeled (BigDye ver 3.1 (tm); Applied Biosystems, Inc. Life Technologies Corporation, Carlsbad, California) terminators in both directions. Sequence data were collected using a 3730xl DNA Analyzer (Applied Biosystems, Inc.). Geneious ver 5.3.5 (Biomatters Ltd., Auckland, New Zealand) was used to edit and create consensus sequences for the forward and reverse reads.

Analytical Approach

In fulfillment of data archiving guidelines (Baker 2013), we have deposited the primary data underlying these analyses with Dryad.

Identification of Species and Subspecies

To identify which Bryde’s whale species or subspecies were present in our sample, we selected mtDNA control region reference sequences for B. e. brydei (large form; ACCN: AB201259, AP006469), B. e. edeni (small form; ACCN: AB201258), and B. omurai (ACCN: AB201256-7) based on the phylogenetic analysis by Sasaki et al. (2006). Sasaki et al. (2006) attempted to phylogenetically verify specimens used in previous studies and obtained new specimens for each taxon that adhered to the classification defined by Wada et al. (2003). As all 3 of these taxa were phylogenetically distinct, and while their correspondence to “small” and “large” forms requires further work, we consider these sequences to represent the most reliable and consistent reference for defining the species and subspecies in this study. We recognize that this situation could change if the genetic identity of the B. edeni holotype is ever determined. Balaenoptera physalus (ACCN: NC_001321.1) was selected as the outgroup for the phylogenetic analysis given its basal evolutionary relationship to the Bryde’s whale complex (Sasaki et al. 2006).

We identified the species and subspecies in our sample using characteristic attribute (CA) diagnosis (Sarkar et al. 2002; Lowenstein et al. 2009). We accepted the phylogeny of the species/subspecies as described by Sasaki et al. (2006), aligned the reference sequences using ClustalW (Higgins et al. 1994) under default settings in MEGA5 (Tamura et al. 2011), and trimmed to the consensus 299bp mitochondrial control region (bp position 15545–15843 in the mtDNA genome of B. e. edeni [ACCN: AB201258]). To construct a character-based key, we visually inspected the reference sequences for variable sites that could serve as diagnostics for the 3 taxa (sensu Lowenstein et al. 2009). We then aligned our unknown sequences to the chosen reference sequences and used the CAs to identify the species and subspecies present in the unknown sample.

Phylogenetic Analysis

Sequences were collapsed to haplotypes using DnaSP ver 5 (Librado and Rosaz 2009). Unique haplotypes were combined with the outgroup B. physalus and a single B. omurai reference sequence (ACCN: AB201256) and aligned using ClustalW (Higgins et al. 1994) under default settings in MEGA5 (Tamura et al. 2011). The resulting 304bp alignment (including gaps) was used to estimate lineage relationships using maximum parsimony (Fitch 1971). Maximum parsimony analysis was conducted in PAUP ver 4.0b10 (Swofford 2002) with 1000 bootstrap replicates using a heuristic search strategy with tree-bisection-reconnection branch swapping, random taxon addition with 100 repetitions and 1 tree held at each step, and a maximum of 1000 trees saved per replicate in order to decrease the time needed to run large bootstrap replicates (Sessa et al. 2012). The resulting bootstrap 50% majority rule consensus tree was edited using Figtree ver 1.3.1 (Rambaut 2009).

Genetic Diversity Indices

For the statistical analyses, haplotypes of B. e. brydei and B. e. edeni were treated separately based on the outcome of the taxon identification and phylogeny, and only sampling regions where n > 5 were included to enable statistical inference. Samples were grouped based on their geographic sampling site using DnaSP ver 5 (Librado and Rosaz 2009). Genetic diversity indices (number of haplotypes, haplotype diversity, nucleotide diversity with Jukes-Cantor correction, and average number of pairwise nucleotide differences among sequences) were calculated in DnaSP for the total sample and for each geographic region (when n > 5). To further explore the genetic diversity of the newly sequenced samples of B. e. edeni from Oman (n = 16) and Bangladesh (n = 29), diversity indices were separately calculated for the consensus 407bp fragment of the mitochondrial control region (bp position 15500–15906 in the mtDNA genome of B. e. edeni [ACCN: AB201258]).

Population-Level Genetic Structure

The Java and Northwest Pacific haplotype frequencies (Yoshida and Kato 1999; Kanda et al. 2007) were combined with the new haplotype frequencies from the NIO. Tests of genetic differentiation between sampling locations (when n > 5) were conducted in Arlequin ver 3.5 (Excoffier and Lischer 2010) for B. e. brydei and B. e. edeni, respectively. A heterogeneity test for haplotype frequencies was calculated using Fisher’s exact test of population differentiation (implemented with 10 000 Markov chain steps and 1000 dememorization steps) at the 0.05 significance level. Pairwise genetic differentiation between sampling sites was calculated using haplotype frequencies (FST) with 1000 permutations at the 0.05 significance level (Weir and Cockerham 1984) in Arlequin ver 3.5 (Excoffier and Lischer 2010). Pairwise genetic distances were calculated in PAUP ver 4.0b10 (Swofford, 2002) assuming the HKY85 model of nucleotide substitution as selected according to the corrected Akaike information criterion implemented in jModelTest ver 2.1 (Darriba et al. 2012). Levels of genetic divergence between samples were then calculated with the fixation index (ΦST) (Excoffier et al. 1992) in Arlequin ver 3.5 using the distance matrix computed in PAUP. Significance of ΦST for all possible pairwise population comparisons was assessed using 1000 permutations at the 0.05 significance level.

Haplotype Network

The dataset for the haplotype network comprised the consensus 299bp control region sequences for B. e. brydei and B. e. edeni, representing 48 haplotypes and 348 samples (including sampling regions with n < 5). The alignment was converted to Roehl data format (.RDF) using DnaSP. Median-joining haplotype networks (Bandelt et al. 1999), both with and without maximum parsimony postprocessing (Mardulyn 2012), were calculated using NETWORK ver 4.6.0.0 (Fluxus Technology Ltd 1999–2010) with ε = 0 and all variable sites weighted equally. Median-joining networks have been recommended over maximum parsimony approaches in intraspecific studies as they capture a greater degree of ambiguity, thus enabling more realistic interpretations (Cassens et al. 2005).

Results

Presence of Bryde’s Whale Species and Subspecies

Phylogenetic reconstruction of available references sequences for B. e. brydei, B. e. edeni, B. omurai, relative to the outgroup B. physalus, identified 9 CAs that were diagnostic of the 4 taxa within the 299bp consensus region. Sequences from our 56 samples matched closely those of B. e. brydei or B. e. edeni, sharing all CAs with one or the other of these taxa. None of the samples matched the known mtDNA sequence of B. omurai or any other species. These taxon-specific (species or subspecies) clades were supported by the maximum parsimony bootstrap 50% majority rule consensus tree based on 41 parsimony informative characters (Figure 2). Bootstrap values for the 2 clades were high (100% for both clades; Figure 2) and support previous work that has identified the 2 subspecies as sister taxa (Sasaki et al. 2006). Samples identified as B. e. brydei and B. e. edeni were therefore treated separately for subsequent diversity and population-level analyses.

Figure 2.

Phylogenetic reconstruction of mtDNA control region haplotypes of Bryde’s whales sampled from across the Western and Central Indo-Pacific and Northwest Pacific Ocean. The bootstrap 50% majority rule consensus parsimony tree is shown with bootstrap values supporting phylogenetic differentiation of haplotypes identified as Balaenoptera edeni brydei and Balaenopteraedeni edeni. The 9 characteristic attributes (CAs) used to identify the taxa are shown to the immediate right of the tree. Nucleotide positions correspond to the B. e. brydei mitochondrial genome positions 15477–16410 (ACCN: AB201259). Positions 15609, 15616, and 15769 diagnose the B. e. brydei subspecies. Positions 15592, 15681, 15722, and 15726 diagnose the B. e. edeni subspecies. * represents conserved nucleotides in relation to the outgroup, Balaenoptera physalus. H, haplotype number; N, sample size; and sampling location (i.e., OMA, Oman; MAL, Maldives; BAN, Bangladesh; JAV, south of Java; COJ, coast of Japan; NWP, Northwest Pacific), are shown adjacent to the termini in the table to the far right. See Supplementary Material online for details of haplotype accession numbers.

Figure 2.

Phylogenetic reconstruction of mtDNA control region haplotypes of Bryde’s whales sampled from across the Western and Central Indo-Pacific and Northwest Pacific Ocean. The bootstrap 50% majority rule consensus parsimony tree is shown with bootstrap values supporting phylogenetic differentiation of haplotypes identified as Balaenoptera edeni brydei and Balaenopteraedeni edeni. The 9 characteristic attributes (CAs) used to identify the taxa are shown to the immediate right of the tree. Nucleotide positions correspond to the B. e. brydei mitochondrial genome positions 15477–16410 (ACCN: AB201259). Positions 15609, 15616, and 15769 diagnose the B. e. brydei subspecies. Positions 15592, 15681, 15722, and 15726 diagnose the B. e. edeni subspecies. * represents conserved nucleotides in relation to the outgroup, Balaenoptera physalus. H, haplotype number; N, sample size; and sampling location (i.e., OMA, Oman; MAL, Maldives; BAN, Bangladesh; JAV, south of Java; COJ, coast of Japan; NWP, Northwest Pacific), are shown adjacent to the termini in the table to the far right. See Supplementary Material online for details of haplotype accession numbers.

Genetic Diversity

The genetic analysis of the mtDNA control region resulted in the identification of 45 unique haplotypes (H1-H45) for B. e. brydei that were derived from 348 sequences with 34 polymorphic sites (2 singletons, 32 parsimony informative) in the 297bp control region (following the removal of gaps and missing data). For sampling locations where n> 5, B. e. brydei (n = 348) was identified at 3 sampling locations: the Maldives (n = 8), south of Java (n = 27), and offshore in the Northwest Pacific (n = 310). In addition, 2 individuals were sampled on the coast of Oman, and 1 individual was sampled from a ship strike offshore of Bangladesh. Genetic diversity (Table 1) was relatively high (Hd: 0.845; π(JC): 0.01319; k: 3.821) and was generally comparable between samples; the Maldives exhibited a relatively lower k value likely related to small sample size, and the south of Java sample exhibited a relatively lower Hd value.

Table 1

Genetic diversity indices for Balaenoptera edeni brydei and Balaenoptera edeni edeni haplotypes for the 299bp consensus region of the total sample and for individual sampling locations where n > 5 (OMA, Oman; MAL, Maldives; BAN, Bangladesh; JAV, south of Java; COJ, coast of Japan; NWP, Northwest Pacific)

Species Sample N S H Hd π(JC) k 
B. e. brydei All 348 34 44 0.844 0.013 3.752 
MAL* 0.750 0.005 1.536 
JAV 27 12 0.396 0.007 2.108 
NWP 310 33 37 0.810 0.012 3.079 
B. e. edeni All 61 0.391 0.004 1.095 
BAN* 29 0.000 0.000 0.000 
OMA* 16 0.000 0.000 0.000 
COJ 16 0.342 0.002 0.575 
Species Sample N S H Hd π(JC) k 
B. e. brydei All 348 34 44 0.844 0.013 3.752 
MAL* 0.750 0.005 1.536 
JAV 27 12 0.396 0.007 2.108 
NWP 310 33 37 0.810 0.012 3.079 
B. e. edeni All 61 0.391 0.004 1.095 
BAN* 29 0.000 0.000 0.000 
OMA* 16 0.000 0.000 0.000 
COJ 16 0.342 0.002 0.575 

New samples are indicated by *. N, number of sequences; S, number of segregating sites; H, number of haplotypes; Hd, haplotype diversity; π(JC), nucleotide diversity with Jukes-Cantor correction; k, average number of pairwise nucleotide differences among sequences.

In contrast, B. e. edeni showed remarkably low genetic diversity with only 3 haplotypes derived in the 299bp control region from 61 sequences (3 parsimony informative sites) (Hd: 0.391; π(JC): 0.00371; k: 1.095; Table 1). For sampling locations where n > 5, B. e. edeni (n = 61) was identified at 3 sampling locations: Bangladesh (n = 29), Oman (n = 16), and coast of Japan (n = 16). Notably, no genetic diversity was found among the Bangladesh and Oman samples as all 45 individuals shared a single haplotype for the 299bp fragment. Three haplotypes were identified in the coast of Japan sample, 1 of which was identical to the haplotype identified in Bangladesh and Oman. When diversity analyses were conducted on the larger 407bp consensus fragment of the mtDNA control region of the new Oman and Bangladesh samples, we identified 1 additional B. e. edeni haplotype in the Oman sample (H49; data not shown).

Overall, 4 new haplotypes were identified for B. e. brydei (H01, H06, H07, and H44; ACCN: JX090150-52, KC261305), and 1 new haplotype was identified for B. e. edeni (H49; ACCN: KC561138). The remaining haplotypes for B. e. brydei and B. e. edeni have been previously found and presented in other studies (Yoshida and Kato 1999; Kanda et al. 2007; see Supplementary Material online).

Population Structure

Median-joining networks showed comparable results irrespective of whether or not maximum parsimony (MP) postprocessing was included. As expected, the median-joining network without MP postprocessing captured a larger number of inferred nodes and reticulations (Cassens et al. 2005; Mardulyn 2012). However, as the fundamental relationships between haplotypes were not affected, only the more parsimonious network with MP postprocessing is shown (Supplementary Material online).

For the 44 haplotypes identified as B. e. brydei, 2 main clusters are apparent: the NIO (Oman, Maldives, and Bangladesh) and the Northwest Pacific. Haplotypes from off Java are represented across the network (Figure 2; Supplementary Material online). Two clusters, NIO and coastal Japan, respectively, are also evident for B. e. edeni. However, a single individual from the coast of Japan was found to share a NIO haplotype (H46). B. e. brydei comprised 11.1% of the total sample in Oman, 100% of the samples in the Maldives, 4.4% of the Bangladesh sample (the single individual sampled from an offshore ship strike), and 100% of the samples from off Java and the Northwest Pacific. In contrast, B. e. edeni was only sampled close to the coastline, comprising 88.9% of the Oman sample, 96.6% in Bangladesh, and 100% in the coastal Japan (Figure 2; Supplementary Material online).

For B. e. brydei, pairwise FST and ΦST values (Table 2) were highly significant between all sampling sites (P < 0.001), indicating that populations in the Maldives, off Java, and the Northwest Pacific can be considered genetically distinct populations. In contrast, for B. e. edeni, pairwise FST and ΦST results showed no significant genetic differentiation between Bangladesh and Oman (FST: 0.000, P > 0.05; ΦST:0.000, P > 0.05). However, highly significant differentiation was found between the coast of Japan and Bangladesh (FST: 0.866, P < 0.001; ΦST:0.923, P < 0.001), as well as Oman (FST: 0.817, P < 0.001; ΦST:0.893, P < 0.001). One haplotype (H46) was shared between all 3 sampling locations and is possibly indicative of some unquantifiable degree of gene flow across the region or the retention of ancestral polymorphism (Figure 2; Supplementary Material online).

Table 2

Pairwise FST and ϕST values for Balaenoptera edeni brydei and Balaenoptera edeni edeni for each sampling location where n > 5 (OMA, Oman; MAL, Maldives; BAN, Bangladesh; JAV, south of Java; COJ, coast of Japan; NWP, Northwest Pacific)

B. e. brydei  MAL JAV NWP 
MAL — 0.479*** 0.211*** 
JAV 0.561*** — 0.334*** 
NWP 0.564*** 0.452*** — 
B. e. edeni  BAN OMA COJ 
BAN — 0.000 0.866*** 
OMA 0.000 — 0.818*** 
COJ 0.923*** 0.893*** — 
B. e. brydei  MAL JAV NWP 
MAL — 0.479*** 0.211*** 
JAV 0.561*** — 0.334*** 
NWP 0.564*** 0.452*** — 
B. e. edeni  BAN OMA COJ 
BAN — 0.000 0.866*** 
OMA 0.000 — 0.818*** 
COJ 0.923*** 0.893*** — 

FST values are shown above the diagonal, ϕST results are shown below the diagonal.

Significance values are indicated as ***P < 0.001 assessed using 1000 permutations at the 0.05 significance level in Arlequin ver 3.5 (Excoffier and Lischer 2010).

Discussion

Our phylogenetic analyses of the mtDNA control region are consistent with previous taxonomic groupings recognized for the subspecies B. e. brydei and B. e. edeni. Our results provide novel insights into the breadth of the distribution of these subspecies across the Western and Central Indo-Pacific and the warm temperate Northwest Pacific and further elucidate genetic patterns at the population level. The striking differences between the 2 forms indicated by these analyses, and when considered alongside previously identified morphological and behavioral differences, support the designation of each form as a separate species or subspecies.

Taxon Identification and Divergence

Using phylogenetic analyses, we confirmed evolutionary divergence in the mitochondrial DNA of Bryde’s whale subspecies within our sample: the offshore, large form, B. e. brydei, and the coastal, small form, B. e. edeni, as previously reported by Kanda et al. (2007) and Sasaki et al. (2006). Due to the limited information available for these taxa, we rely solely on the best available genetic data to define the species and subspecies in our study. Reference sequences should ideally be based on verified voucher specimens that offer corollary information (e.g., morphological data) for taxon designation (Reeves et al. 2004), and we recognize that this is a limitation of our study; independent classification using morphological data is required for formal taxonomic classification (Reeves et al. 2004; DeSalle et al. 2005).

Individual genetic loci, like morphological characters, do not necessarily reflect the true phylogenetic history; the gene tree is not always consistent or congruent with the species tree (Page and Charleston 1997). This has been previously demonstrated in Bryde’s whales by Sasaki et al. (2006) who found inconsistencies in the phylogenetic relationships between B. e. brydei, B. e. edeni, and B. borealis dependent upon the mitochondrial molecular marker employed. Therefore, the phylogeny we identified is likely to be, at least in part, a function of the single mtDNA marker used in the analysis. Future analyses utilizing larger fragments of the mitochondrial genome alongside additional nuclear markers are likely to further resolve the phylogenetic relationships of the Bryde’s whale species complex (Morin et al. 2010; Sharma et al. 2012).

Morphological, behavioral, and geographic information indicate strong differences between the 2 subspecies. This differentiation is not only of ecological and evolutionary interest, but is also of critical importance for informing the conservation and management of these whales. Size differences and temporal reproductive phase shifts have been recorded in historical whaling data (Mikhalev 2000). Observations of habitat partitioning (i.e., coastal vs. offshore) between the 2 subspecies (Best 2001) indicate the existence of an ecological barrier to gene flow, which may have acted as the mechanism for divergence. These findings are corroborated by field observations of a putative population of coastal small-form whales off South Africa (Best 2001; Penry et al. 2011); however, the genetic identity of this group still needs to be confirmed. This study provides further evidence by showing that B. e. brydei appears to have a more cosmopolitan distribution in both coastal and offshore areas, likely due to greater mobility and offshore habitat use. In contrast, B. e. edeni was only sampled close to the coast of Japan indicating that the coastal waters of the Northwest Pacific may represent their eastern and northern range extent in the North Pacific.

The original 9 specimens of B. omurai were from the Solomon Sea (n = 6) in 1976, off the Cocos Islands (n = 2) in 1978, and Tsunoshima (34°21′N, 130°52′E), Sea of Japan, Japan (n = 1) in 1998 (Wada et al. 2003). More recently, additional specimens have been reported from southern Japan, Taiwan, Philippines, and Thailand (the westernmost specimen of B. omurai from the Andaman Sea). However, the identification of these new specimens of B. omurai is based solely on their morphology and not genetics (Yamada et al. 2006, 2008). Omura’s whale and B. edeni, therefore, appear to be sympatric in parts of their range off southern Japan, Taiwan, and off Thailand in the Andaman Sea. This sympatry may also occur in the waters around Cocos Islands in the eastern central Indian Ocean where 2 specimens of B. omurai were taken under a special research permit in the late 1970s (Wada and Numachi 1991). The exact details of any habitat sympatry are unknown because all the whales from Japan, Taiwan, and Thailand are based on stranded specimens. The lack of B. omurai in our sample adds support to the western limit of this species in the Eastern Indian Ocean being the Andaman Sea, off the western coast of the Malay Peninsula (Yamada et al. 2008, Yamada 2009).

Population-Level Diversity and Structure

The genetic structure observed for B. e. brydei indicates 3 discrete populations experiencing very little gene flow in the Maldives, off Java, and the Northwest Pacific. We note, however, that the small sample size for the Maldives (n = 8) limits the statistical inference that can be made regarding this potential “population” and precludes a definitive conclusion. Given the potential consequences of not recognizing a genetically differentiated group in a species subject to continued hunting, we chose to include the Maldives as a separate population unit in this study as a precautionary measure with the view to informing management.

The population identity of the whales off Java is not clear, as 3 of the 5 haplotypes were also identified within the Maldives sample (n = 1) and the Northwest Pacific sample (n = 2), indicating contemporary or historic gene flow. Interestingly, the Java population also exhibits much lower genetic diversity (Hd = 0.396) than either the Maldives (Hd = 0.750) or the Northwest Pacific (Hd = 0.810), suggesting that the population may be small and subject to the effects of genetic drift, perhaps due to the lower ocean productivity found in this region (Longhurst 2007). Two whales sampled in Oman were identified as B. e. brydei, suggesting that another discrete population may exist in the Arabian Sea, or that the population identified in the Maldives may have a broader geographical range than detected by this analysis. Historical Soviet whaling records report large aggregations of both large- and small-form Bryde’s whales in the Gulf of Aden (Mikhalev 2000), indicating that this may indeed be an important part of the range for both of these taxa. Increased genetic sampling in this region will be crucial in delineating population boundaries for management purposes.

In marked contrast to B. e. brydei, extremely low degrees of genetic diversity (Hd = 0.391) and population structure were found for B. e. edeni across the NIO, at a scale not before seen in baleen whales (Rosenbaum et al. 2000; Patenaude et al. 2007). Only a single haplotype (Figure 2; Supplementary Material online; H46) was shared between the 45 individuals sampled in Bangladesh (Hd = 0.000) and Oman (Hd = 0.000) when the 299bp consensus sequence was examined. As only 1 additional haplotype was identified in Oman (when the larger 407bp fragment of the control region was considered), these low levels of diversity are likely not fully explained by the limited length of the marker used in our study. Notably, no further diversity was found within the Bangladesh sample, indicating that levels of genetic diversity can still be considered unusually low for this subspecies.

We observed strong population structure between the NIO populations of both B. e. brydei and B. e. edeni compared with those in the Northwest Pacific and in the coastal waters of Japan (Table 2; FST and ΦST have significance values of P < 0.001 for all comparisons.). This is consistent with the biogeographic barrier imposed by the peninsulas and islands of Thailand, Malaysia, and Indonesia. However, the shared haplotype between Java and the Northwest Pacific for B. e. brydei (Figure 2; Supplementary Material online; H39), and between the NIO and coast of Japan for B. e. edeni (Figure 2; Supplementary Material online; H46), provides evidence of interoceanic exchange, at least historically, within populations of both taxa. Given our small sample size, it can be assumed that we underestimate the actual rates of genetic exchange between the NIO and the Northwest Pacific, thus implying a porous barrier to long-range movements.

Conclusion

Evidence from phylogenetic analyses, and corroborating morphological and behavioral studies, supports the presence of 2 taxonomic units of Bryde’s whale across the Western and Central Indo-Pacific and the Northwest Pacific Ocean. The distinctiveness of the 2 subspecies confirms the need to designate each taxon as a separate conservation unit with specific management recommendations for each. Bryde’s whales are vulnerable to fisheries bycatch and ship strikes across the study region (Bijukumar et al. 2012) and are currently subject to scientific whaling by Japan in the western North Pacific. There is also the potential impact of hydrocarbon exploration and development in coastal waters. Given these stressors, there is a clear need to implement effective management measures that are fully informed by better defining conservation units at the species and population level using molecular information.

Strong genetic differences were found at the population level within B. e. brydei and B. e. edeni. We found significant differentiation among populations of B. e. brydei in the Maldives, Java, and in the offshore Northwest Pacific and B. e. edeni off Oman and Bangladesh in the NIO and in the coastal waters of southern Japan. We therefore suggest that each population be considered an independent conservation unit for management purposes. The Arabian Sea may also represent an important priority for management given bycatch and ship strikes of these whales in the region, and the catches of 849 Bryde’s whales during the mid-1960s, which based on their total lengths would likely be B. e. brydei (Mikhalev 2000). This is a priority for future research as it cannot yet be determined if the whale populations in the Arabian Sea are independent of the Maldives unit identified in this study. Additional genetic sampling is therefore urgently needed in the Arabian Sea and the Maldives, as well as coastal Southeast Asia, particularly along the Malay Peninsula and in the Gulf of Thailand (Perrin and Brownell 2007).

In addition, biparentally inherited, neutral microsatellite markers and single-nucleotide polymorphisms identified by high throughput sequencing techniques represent powerful future tools to complement population-level mtDNA analyses. Longer mtDNA sequences are likely to provide greater resolution of haplotypes and more informative estimates of genetic diversity and population differences, as indicated by our identification of an additional B. e. edeni haplotype in Oman. It will also be important to collect additional morphological information to validate the findings of phylogenetic and population genetic studies. Photographic documentation of individuals during biopsy sampling and the collection of morphological information from future ship strikes in the Indian Ocean represent 2 opportunistic methods to gather additional data. The application of these new data will enable the finer-scale, spatiotemporal analyses essential for ensuring appropriate management and persistence of these whales (Gaines et al. 2005; Dale and Von Schantz 2007).

Supplementary Material

Supplementary material can be found at http://www.jhered.oxfordjournals.org/.

Funding

Funding to HCR was provided by numerous individuals and foundations.

Acknowledgments

We are grateful to the staff of WCS and AMNH, WCS staff in Bangladesh, colleagues from the Environment Society of Oman, and Robert L. Pitman and the director and staff of the Marine Research Centre for assistance with collection of samples from the Maldives. Skin tissue biopsy collection by primary authors was carried out under approvals from relevant permitting agencies. We are particularly grateful to all the staff of the Sackler Institute for Comparative Genomics, especially George Amato, Rob DeSalle, and Rebecca Hersch. We thank all the summer interns, students, volunteers, and research assistants who have contributed to field and laboratory work. Logistical support in Oman and Bangladesh was provided by a number of dedicated colleagues. We would like to thank Naohisa Kanda for the provision of haplotype frequencies for JAV and NWP samples used from Kanda et al. (2007). We are grateful to William F. Perrin and Ana Rita Amaral for feedback on the final version of the manuscript and would like to thank Martin Mendez and Sergios Kolokotronis for their contributions to the early stages of the analyses presented here.

References

Baker
CS
.
2013
.
Journal of heredity adopts joint data archiving policy
.
J Hered
 .
104
:
1
.
Baker
CS
Medrano-Gonzalez
L
Calambokidis
J
Perry
A
Pichler
F
Rosenbaum
H
Straley
JM
Urban-Ramirez
J
Yamaguchi
M
von Ziegesar
O
.
1998
.
Population structure of nuclear and mitochondrial DNA variation among humpback whales in the North Pacific
.
Mol Ecol
 .
7
:
695
707
.
Baker
CS
Perry
A
Bannister
JL
Weinrich
MT
Abernethy
RB
Calambokidis
J
Lien
J
Lambertsen
RH
Ramírez
JU
Vasquez
O
.
1993
.
Abundant mitochondrial DNA variation and world-wide population structure in humpback whales
.
Proc Natl Acad Sci USA
 .
90
:
8239
8243
.
Bandelt
HJ
Forster
P
Röhl
A
.
1999
.
Median-joining networks for inferring intraspecific phylogenies
.
Mol Biol Evol
 .
16
:
37
48
.
Best
PB
.
1997
.
Two allopatric forms of Bryde’s whale off south Africa
.
Rep Int Whal Comm
 .
39
:
363
369
.
Best
PB
.
2001
.
Distribution and population separation of Bryde’s whale Balaenoptera edeni off southern Africa
.
Mar Ecol Prog Ser
 .
220
:
277
289
.
Bijukumar
A
Jijith
SS
Suresh Kumar
U
George
S
.
2012
.
DNA barcoding of the Bryde’s whale Balaenoptera edeni Anderson (Cetacea: Balaenopteridae) washed ashore along Kerala coast, India
.
J Threatened Taxa
 .
4
:
2436
2443
.
Cassens
I
Mardulyn
P
Milinkovitch
MC
.
2005
.
Evaluating intraspecific “network” construction methods using simulated sequence data: do existing algorithms outperform the global maximum parsimony approach?
Syst Biol
 .
54
:
363
372
.
Clapham
PJ
Agular
A
Hatch
LT
.
2008
.
Determining spatial and temporal scales for management: lessons from whaling
.
Mar Mamm Sci
 .
24
:
183
201
.
Committee on Taxonomy
.
2011
.
List of marine mammal species and subspecies
. [cited 2012 May 16]. Available from: http://www.marinemammalscience.org
Dale
J
Von Schantz
M
.
2007
.
From genes to genomes: concepts and applications of DNA technology
 .
West Sussex (UK)
:
John Wiley & Sons Ltd
.
Darriba
D
Taboada
GL
Doallo
R
Posada
D
.
2012
.
jModelTest 2: more models, new heuristics and parallel computing
.
Nat Methods
 .
9
:
772
.
DeSalle
R
Egan
MG
Siddall
M
.
2005
.
The unholy trinity: taxonomy, species delimitation and DNA barcoding
.
Philos Trans R Soc Lond B Biol Sci
 .
360
:
1905
1916
.
Excoffier
L
Lischer
HE
.
2010
.
Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows
.
Mol Ecol Resour
 .
10
:
564
567
.
Excoffier
L
Smouse
PE
Quattro
JM
.
1992
.
Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data
.
Genetics
 .
131
:
479
491
.
Fitch
WM
.
1971
.
Toward defining the course of evolution: minimum change for specific tree topology
.
Systematic Zool
 .
20
:
406
416
.
Fluxus Technology Ltd
.
1999–2010
.
NETWORK ver 4.6.0.0
. [cited 2012 February 10]. Available from: http://www.fluxus-engineering.com
Gaines
CA
Hare
MP
Beck
SE
Rosenbaum
HC
.
2005
.
Nuclear markers confirm taxonomic status and relationships among highly endangered and closely related right whale species
.
Proc Biol Sci
 .
272
:
533
542
.
Higgins
D
Thompson
J
Gibson
T
Thompson
JD
Higgins
DG
Gibson
TJ
.
1994
.
CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting penalties and weight matrix choice
.
Nucleic Acids Res
 .
22
:
4673
4680
.
Kanda
N
Goto
M
Kato
H
McPhee
MV
Pastene
LA
.
2007
.
Population genetic structure of Bryde’s whales (Balaenoptera brydei) at the inter-oceanic and trans-equatorial levels
.
Conserv Genet
 .
8
:
853
864
.
Kato
H
Perrin
WF
.
2009
.
Bryde’s whales Balaenoptera edeni/brydei
. In:
Perrin
WF
Würsig
B
Thewissen
JGM
, editors.
Encyclopedia of marine mammals
 .
Amsterdam
:
Academic Press
. p.
158
163
.
LeDuc
RG
Dizon
AE
Goto
M
Pastene
LA
Kato
H
Nishiwaki
S
LeDuc
CA
Brownell
RL
Jr
.
2007
.
Patterns of genetic variation in Southern Hemisphere blue whales and the use of assignment test to detect mixing on the feeding grounds
.
J Cetacean Res Manage
 .
9
:
73
80
.
Librado
P
Rozas
J
.
2009
.
DnaSP v5: a software for comprehensive analysis of DNA polymorphism data
.
Bioinformatics
 .
25
:
1451
1452
.
Longhurst
A
.
2007
.
Ecological geography of the sea
 .
6th ed
.
Oxford
:
Elsevier
.
Lowenstein
JH
Amato
G
Kolokotronis
SO
.
2009
.
The real maccoyii: identifying tuna sushi with DNA barcodes–contrasting characteristic attributes and genetic distances
.
PLoS One
 .
4
:
e7866
.
Mardulyn
P
.
2012
.
Trees and/or networks to display intraspecific DNA sequence variation?
Mol Ecol
 .
21
:
3385
3390
.
Mendez
M
Rosenbaum
HC
Subramaniam
A
Yackulic
C
Bordino
P
.
2010
.
Isolation by environmental distance in mobile marine species: molecular ecology of franciscana dolphins at their southern range
.
Mol Ecol
 .
19
:
2212
2228
.
Mikhalev
YA
.
2000
.
Whaling in the Arabian Sea by the whaling fleets Slava and Sovetskaya Ukraina
. In:
Zemsky
VA
Yablokov
AV
, editors.
Soviet whaling data [1949–1979]
 .
Moscow (Russia)
:
Center for Russian Environmental Policy, Marine Mammal Council
. p.
141
181
.
Morin
PA
Archer
FI
Foote
AD
Vilstrup
J
Allen
EE
Wade
P
Durban
J
Parsons
K
Pitman
R
Li
L
et al
.
2010
.
Complete mitochondrial genome phylogeographic analysis of killer whales (Orcinus orca) indicates multiple species
.
Genome Res
 .
20
:
908
916
.
Page
RD
Charleston
MA
.
1997
.
From gene to organismal phylogeny: reconciled trees and the gene tree/species tree problem
.
Mol Phylogenet Evol
 .
7
:
231
240
.
Patenaude
NJ
Portway
VA
Schaeff
CM
Bannister
JL
Best
PB
Payne
RS
Rowntree
VJ
Rivarola
M
Baker
CS
.
2007
.
Mitochondrial DNA diversity and population structure among southern right whales (Eubalaena australis)
.
J Hered
 .
98
:
147
157
.
Penry
GS
Cockroft
VG
Hammond
P
.
2011
.
Seasonal fluctuations in occurrence of inshore Bryde’s whales in Plettenberg Bay, South Africa, with notes on feeding and multispecies associations
.
Afr J Mar Sci.
 
33
:
403
414
.
Perrin
WF
Brownell
RL
Jr
.
2007
.
Proposed updates to the List of Recognised Species of Cetaceans
.
Rep Int Whal Comm
 .
SC/59/O15
:
1
4
.
Perrin
WF
Dolar
MLL
Ortega
E
.
1996
.
Osteological comparison of Bryde’s whales from the Philippines with specimens from other regions
.
Rep Int Whal Comm
 .
46
:
409
413
.
Rambaut
A
.
2009
.
Figtree: Tree Figure Drawing Tool, Version 1.3.1
 .
Edinburgh (UK)
:
Institut of Evolutionary Biology, University of Edinburgh
.
Reeves
RR
Perrin
WF
Taylor
BL
Baker
CS
Mesnick
SL
.
2004
.
Report of the workshop on shortcomings of cetacean taxonomy in relation to needs of conservation and management
 , April 30–May 2, 2004, La Jolla, California, NOAA Tech Mem. NOAA-NMFS-SWFSC-363:59–60.
Silver Spring (MD)
:
NOAA
.
Rosenbaum
HC
Egan
MG
Clapham
PJ
Brownell
RL
Jr
Malik
S
Brown
MW
White
BN
Walsh
P
Desalle
R
.
2000
.
Utility of North Atlantic right whale museum specimens for assessing changes in genetic diversity
.
Conserv Biol
 .
14
:
1837
1842
.
Rosenbaum
HC
Pomilla
C
Mendez
M
Leslie
MS
Best
PB
Findlay
KP
Minton
G
Ersts
PJ
Collins
T
Engel
MH
et al
.
2009
.
Population structure of humpback whales from their breeding grounds in the South Atlantic and Indian Oceans
.
PLoS One
 .
4
:
e7318
.
Sarkar
IN
Thornton
JW
Planet
PJ
Figurski
DH
Schierwater
B
DeSalle
R
.
2002
.
An automated phylogenetic key for classifying homeoboxes
.
Mol Phylogenet Evol
 .
24
:
388
399
.
Sasaki
T
Nikaido
M
Wada
S
Yamada
TK
Cao
Y
Hasegawa
M
Okada
N
.
2006
.
Balaenoptera omurai is a newly discovered baleen whale that represents an ancient evolutionary lineage
.
Mol Phylogenet Evol
 .
41
:
40
52
.
Sessa
EB
Zimmer
EA
Givnish
TJ
.
2012
.
Phylogeny, divergence times, and historical biogeography of New World Dryopteris (Dryopteridaceae)
.
Am J Bot
 .
99
:
730
750
.
Sharma
R
Goossens
B
Kun-Rodrigues
C
Teixeira
T
Othman
N
Boone
JQ
Jue
NK
Obergfell
C
O’Neill
RJ
Chikhi
L
.
2012
.
Two different high throughput sequencing approaches identify thousands of de novo genomic markers for the genetically depleted Bornean elephant
.
PLoS One
 .
7
:
e49533
.
Spalding
M
Fox
HE
Allen
GR
Davidson
N
Ferdaña
ZA
Finlayson
M
Halpern
BS
Jorge
MA
Lombana
AL
Lourie
SA
et al
.
2007
.
Marine ecoregions of the world: A bioregionalization of coastal and shelf areas
.
Bioscience
 .
57
:
573
583
.
Stevick
PT
Allen
J
Clapham
PJ
Katona
SK
Larsen
F
Lien
J
Mattila
DK
Palsbøll
PJ
Sears
R
Sigurjónsson
J
et al
.
2006
.
Population spatial structuring on the feeding grounds in North Atlantic humpback whales (Megaptera novaeangliae)
.
J Zool
 .
270
:
244
255
.
Swofford
DL
.
2002
.
PAUP*. Phylogenetic Analysis Using Parsimony (*and other methods)
 .
Sunderland (MA)
:
Sinauer Associates
.
Tamura
K
Peerson
D
Peterson
N
Stecher
G
Nei
M
Kumar
S
.
2011
.
MEGA5: Molecular Evolutionary Genetics Analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods
.
Mol Biol Evol
 .
28
:
2731
2739
.
Taylor
B
.
2005
.
Identifying units to conserve
. In:
Reynolds
JE
Perrin
WF
Reeves
RR
Montgomery
S
Ragen
TJ
, editors.
Marine mammal research, conservation beyond crisis
 .
Baltimore (MD)
:
John Hopkins University Press
. p.
149
164
.
Wada
S
Numachi
K
.
1991
.
Allozyme analyses of genetic differentiation among populations and species of the Balaenoptera
.
Rep Int Whal Comm
 .
13
(special issue):
125
154
.
Wada
S
Oishi
M
Yamada
TK
.
2003
.
A newly discovered species of living baleen whale
.
Nature
 .
426
:
278
281
.
Weir
BS
Cockerham
CC
.
1984
.
Estimating F-statistics for the analysis of population structure
.
Evolution
 .
38
:
1358
1370
.
Yamada
TK
.
2009
.
Balaenoptera omurai
. In:
Ohdachi
SD
Ishibashi
Y
Iwasa
MA
Saitoh
T
, editors.
The wild mammals of Japan
 .
Kyoto (Japan)
:
Shoukadoh
.
Yamada
TK
Chou
L-S
Chantrapornsyl
S
Adulyanukosol
K
Chakravarti
SK
Oishi
M
Wada
S
Yao
C-J
Kakuda
T
Tajima
Y
et al
.
2006
.
Middle sized Balaenopterid whale specimens (Cetacea: Balaenopteridae) preserved at several institutions in Taiwan, Thailand, and India
.
Mem Natl Sci Mus Tokyo
 .
44
:
1
10
.
Yamada
TK
Kakuda
T
Tajima
Y
.
2008
.
Middle sized balaenopterid whale specimens in the Philippines and Indonesia
.
Mem Natl Sci Mus Tokyo
 .
45
:
75
83
.
Yoshida
H
Kato
H
.
1999
.
Phylogenetic relationships of Bryde’s whales in the Western North Pacific and adjacent waters inferred from mitochondrial DNA sequences
.
Mar Mamm Sci
 .
15
:
1269
1286
.

Author notes

Corresponding Editor: C. Scott Baker