Biodiversity of the Betta smaragdina (Teleostei: Perciformes) in the northeast region of Thailand as determined by mitochondrial COI and nuclear ITS1 gene sequences.

In Thailand, there are currently five recognized species members of the bubble-nesting Betta genus, namely Betta splendens, B. smaragdina, B. imbellis, B. mahachaiensis and B. siamorientalis. In 2010, we indicated the possibility, based on COI barcoding evidence, that there might be two additional species, albeit cryptic, related to the type-locality B. smaragdina in some provinces in the northeast of Thailand. In the present study, after a more extensive survey of the northeast, and phylogenetic analyses based on COI and ITS1 sequences, the B. smaragdina group may be composed of at least 3 cryptic species members. The phylogenetic positions of these B. smaragdina group members in the bubble-nesting bettas' tree together with those of their congeners have been consolidated by better DNA sequence quality and phylogenetic analyses. With a better supported tree, the species statuses of B. siamorientalis and the Cambodian B. smaragdina-like fish, B. stiktos, are also confirmed.

Introduction

Fighting fish belonging to the genus betta are native to most ASEAN countries except the Philippines. Identifying new species of the betta fighting fish, the bubble nesters and the mouth brooders, in these Southeast Asian countries is still an active undertaking with most reports basing their species differentiation on morphological criteria (Schindler and Schmidt, 2008, Tan, 2009a, Tan, 2009b, Kowasupat et al., 2012a, Kowasupat et al., 2012b, Schindler and Linke, 2013, Tan, 2013).

In the literature up to 2010, it was established, by morphological criteria, that there were three species of bubble-nesting fighting fish in Thailand viz B. splendens Regan, 1909 (the Siamese fighting fish), B. imbellis Ladiges, 1975 and Betta smaragdina Ladiges, 1972. B. imbellis was found to locate in the southern tip of the Thailand peninsula. B. smaragdina was found in northeastern Thailand (called the E-sarn region), whereas B. splendens was widespread in the north, central plain, eastern, western, and upper southern region (Lertpanich and Aranyavalai, 2007, Monvises et al., 2009). Incidentally these more widespread Betta splendens fish have mainly been bred by selection into ornamental fish with exaggerated unpaired fins in terms of color, color pattern, and length, and fighter fish with regular-shaped unpaired fins which are small relative to body size. Another species of bubble-nesting betta from a Cambodian town bordering Thailand, B. stiktos, was identified by Tan and Ng (2005). An RAPD (randomly amplified polymorphic DNA) study of a limited number of betta fish in Thailand was reported in 2005 by Tanpitayacoop and Na-Nakorn (2005). In 2010, the phylogenetic relationship of three fighting fishes in Thailand was studied by Khongnomnan et al. (2010) using COI and 12S rRNA genes. In addition, Rüber et al. (2004a) used mitochondrial and nuclear genes to reconstruct a phylogenetic tree for a large number of mouth-brooding and bubble-nesting bettas; three of the latter were B. splendens, B. imbellis, and B. smaragdina.

We reported in 2012 two new species of bubble-nesting bettas, B. mahachaiensis and B. siamorientalis, in central and eastern Thailand based on morphological characters as well as short pieces of a mitochondrial COI (cytochrome c oxidase subunit I) gene and a nuclear ITS1 (internal transcribed spacer 1) DNA (Kowasupat et al., 2012a, Kowasupat et al., 2012b). Previously in 2010, in the preliminary work on the Betta sp. Mahachai using COI and 16S rRNA genes (Sriwattanarothai et al., 2010), we reported that, in the northeastern region of Thailand, there were two cryptic species of B. smaragdina distinct from the B. smaragdina belonging to the type-locality established earlier by Ladiges (1972).

Except one species, all wild bubble-nesting betta fish are still widespread despite human activities which are making their natural habitats less accessible in many places. The one major exception is Betta mahachaiensis whose already small and unique habitats have been dwindling fast also because of their proximity to the capital Bangkok and the accompanying human activities.

Despite barcoding based on COI having been proven successful in species identification in many cases (Ward et al., 2009, Zemlak et al., 2009, Pereira et al., 2013), its sole use has received criticisms for being inadequate in several ways (Dasmahapatra et al., 2010, Taylor and Harris, 2012). There are researchers who use both mitochondrial and nuclear DNA for species identification (Rüber et al., 2004a, Rüber et al., 2006), a practice that we followed in this study. To facilitate our discussion about the highly diversified smaragdina fish and other bettas in this study, we would like to define, for the present purpose, the terms B. smaragdina group, B. smaragdina type-locality complex, B. splendens group, B. splendens complex, and B. imbellis complex by adapting the indented classification recommended by Kizirian and Donnelly (2004) to the reconstructed phylogenetic tree in Fig. 1 (to be discussed) as shown in Table 1.

Fig. 1

The phylogenetic tree of the bubble-nesting bettas (belonging to the B. splendens group and the B. smaragdina group) reconstructed from the combined dataset of aligned COI and ITS1 sequences using partitioned Bayesian inference with reversible jump Metropolis-coupled Markov chain Monte Carlo (rjMC3): the posterior probabilities are indicated next to the nodes. Betta ocellata was used as the outgroup. The length of the scale bar is equivalent to 0.03 substitutions per nucleotide position. The horizontal lengths (from apex to base) of the isosceles triangles represent the most extreme divergences within the clades whereas the vertical lengths (length of base) represent the number of specimens. Each terminating branch without a triangle represents an individual specimen.

Table 1

Terms describing all the bettas in this study based on Fig. 1 and their relationships.

Betta splendens group

B. splendens complex

B. splendens

B. mahachaiensis

B. imbellis complex

B. imbellisa

B. siamorientalis

Betta smaragdina group

B. smaragdina type-locality complex

B. smaragdina (type-locality)a

B. stiktos

B. sp. (cf. smaragdina) 1

B. sp. (cf. smaragdina) 3

B. sp. (cf. smaragdina) 4

The Betta smaragdina group refers to all members of the group except B. stiktos.

Paraphyletic clade.

Here we report our new findings, based principally on COI and ITS1 DNA sequences, confirming the two cryptic species reported previously and one additional cryptic species of the B. smaragdina group in northeast Thailand. We also show new and well supported phylogenetic relationships among the members of the B. smaragdina group and also their relationships with respect to the members of the B. splendens group.

Materials and methods

Sample collection

The specimens of wild Thai B. smaragdina were collected from 19 out of 20 provinces in the northeast region, eight of which from the upper northern one third (in terms of area) (Sakon Nakhon basin) and eleven from the southern two thirds (Korat basin), whereas additional B. smaragdina fish were collected in Lao PDR from 2 provinces in the north and 2 provinces in the south. Wild B. imbellis were collected from 9 provinces in the southern tip of Thailand peninsula and from Malaysia (the type-locality one). Wild B. mahachaiensis were from central provinces of Samut Sakhon, Bangkok, and Samut Prakan. Wild B. splendens, the most widely distributed species, were caught from 15 provinces in the central plain, 4 provinces in the north, 4 provinces in the west, 4 provinces in the east, and 2 provinces in the upper peninsula of Thailand: Bangkok, the capital, is considered to locate at the center of the country. Wild B. siamorientalis were collected from eastern provinces of Chachoengsao, Sa Kaeo, Prachin Buri, Chon Buri, and a Cambodian province, Banteay Meanchey, bordering Thailand. B. stiktos were from Stung Treng, Cambodia. About 5 to 15 sites in each province were surveyed. Although betta fish are widespread in Thailand, we did not find the fish in a few provinces. The number of fish collected from each province was about 10 to 50. The rural ones were collected from sites far away from human settlements to better ensure that the fish were not those previously caught, reared, and later deemed undesirable and discarded haphazardly into natural waters.

The sacrificed fish were preserved in 95% ethanol for genetic analysis. Fish specimens were deposited at the National Science Museum, Thailand.

Sequence data were submitted to GenBank. Of the total 555 COI sequences, 249 were from Sriwattanarothai et al., 2010 (www.barcodinglife.org), plus 306 new COI sequences with accession numbers JQ818641-JQ818818, KF278818-KF278941 and KF381319-KF381322 (www.ncbi.nlm.nih.gov). All the 102 ITS1 sequences have accession numbers JQ818556-JQ818640 and KF381292-KF381318 (www.ncbi.nlm.nih.gov). The Kimura 2-parameter (K2P) distance metric was used for intraspecific and interspecific sequence comparisons (Kimura, 1980).

DNA extraction and sequencing

DNA extraction and sequencing were performed in two laboratories in Thailand: one at Prince of Songkla University (direct sequencing method), and the other at Mahidol University (DNA cloning method). The two methods slightly differed in terms of detailed conditions and primers used for ITS1. However, these two diverse methods gave congruent results which were integrated to recreate phylogenetic trees.

Direct sequencing method

DNA was extracted from muscle tissue by using the DNA extraction kit of Stratagene (Agilent Technologies). The barcoding region of COI was amplified as described by Ivanova et al. (2007) using VF2_t1, FishF2_t1 and FishR2_t1 (Ward et al., 2005), and FR1d_t1 (Ivanova et al., 2007) as primers (see Table 2).

Table 2

Primers used in PCR amplification and DNA sequencing of the COI and ITS1 regions.

Targeted gene	Primer name	Primer sequence (5′–3′)	References
COI	VF2_t1	TGTAAAACGACGGCCAGTCAACCAACCACAAAGACATTGGCAC	Ward et al. (2005)
COI	FishF2_t1	TGTAAAACGACGGCAGTCGACCTAATCATAAAGATATCGGCAC	Ward et al. (2005)
COI	FishR2_t1	CAGGAAACAGCTATGACACTTCAGGGTGACCGAAGAATCAGAA	Ward et al. (2005)
COI	FR1d_t1	CAGGAAACAGCTATGACACCTCAGGGTGTCCGAARAAYCARAA	Ivanova et al. (2007)
COI	M13F (-21)	TGTAAAACGACGGCCAGT	Messing (1983)
COI	M13R (-27)	CAGGAAACAGCTATGAC	Messing (1983)
ITS1	Betta_ITS1_F1	CACACCGCCCGTCGCTACTA	This study
ITS1	Betta_ITS1_R1	GTYCTTCMTCGACSCACGAG	This study
ITS1	Betta_ITS1_F2	ACTTGACTATCTAGAGGAAG	This study
ITS1	Betta_ITS1_R2	GTTCTTCATCGACGCACGAG	This study
ITS1	Betta_ITS1_F3	GTAGGTGAACCTGCGGAAGRATCATT	This study
ITS1	Betta_ITS1_R3	CACGAGYCGAGTGATCCACCGCTAA	This study
ITS1	Betta_ITS1_R4	TCCACCGCTAAGAGTTGTC	This study

Amplification was carried out with the total volume of 50 μl reaction mixture containing 500 ng genomic DNA, 1 × PCR buffer (200 mM Tris–HCl pH 8.4, 500 mM KCl), 2.5 mM MgCl2, 0.05 mM dNTPs, 1.2 U of Platinum® Taq DNA Polymerase (Invitrogen™), and 0.2 μM of each primer. PCR amplification was performed using a cycling program with initial denaturation at 95 °C for 2 min, followed by 40 cycles of 94 °C for 0.5 min, annealing at 52 °C for 0.5 min, extension at 72 °C for 1 min, and final extension step at 72 °C for 10 min. The PCR product was purified by running in 0.1% agarose gel. DNA sequencing was performed by 1st BASE DNA Sequencing Services (Malaysia) using sequencing primer M13F or M13R (Messing, 1983) (see Table 2).

The nuclear ITS1 DNA was amplified and sequenced with sequencing primers Betta_ITS1_F1 or Betta_ITS1_R1 (designed by Center for Genomics and Bioinformatics Research, Prince of Songkla University, Thailand) (see Table 2). The PCR reaction mixtures were carried out in a volume of 50 μl with 500 ng genomic DNA, 1 × PCR buffer (200 mM Tris–HCl pH 8.4, 500 mM KCl), 1.5 mM MgCl2, 0.2 mM dNTPs, 2 U Platinum® Taq DNA Polymerase (Invitrogen™) and 0.8 μM of each primer. The cycling program started with initial denaturation at 94 °C for 5 min, followed by 40 cycles of denaturing at 94 °C for 1 min, annealing at 54 °C for 1 min, extension at 72 °C for 1 min, and final extension step at 72 °C for 10 min. DNA sequencing was performed by 1st BASE DNA Sequencing Services (Malaysia) or Macrogen DNA Sequencing Services (Korea).

DNA cloning method

DNA was extracted from muscle tissue using approximately 5 mg of the fish specimen incubating with 180 μl of 50 mM NaOH at 95 °C for 10 min. The pH of the extract was then neutralized by addition of 20 μl of 1 M Tris–HCl, pH 8.0. After centrifugation, the supernatant having the DNA template was used for PCR reactions.

The barcoding region of COI was amplified as described by Ivanova et al. (2007) using VF2_t1, FishF2_t1 and FishR2_t1 (Ward et al., 2005), and FR1d_t1 (Ivanova et al., 2007) as primers (see Table 2) except that M13 tails were omitted from the primer sequences. The nuclear ITS1 DNA was amplified by primers Betta_ITS1_F2 and Betta_ITS1_R2, Betta_ITS1_F3 and Betta_ITS1_R3, or Betta_ITS_F2 and Betta_ITS_R4 designed by Mahidol University (see Table 2).

Amplification was carried out with the total volume of 25 μl containing 2 μl of the crude extract, 0.6 μl of 10 μM of each COI primer or 0.75 μl of 10 μM of each ITS1 primer, 0.5 μl of Terra™ PCR Direct Polymerase Mix (1.25 U/μl, Clontech), 10 μl of 2 × Terra™ PCR Direct Buffer (with Mg2 + and dNTP), and sterile distilled water. The PCR conditions were 98 °C for 2 min followed by 35 cycles of 98 °C for 10 s, annealing at 52 °C for 15 s and elongation at 68 °C for 1 min/kb.

After purification, the PCR products were ligated to pPrime cloning vector (5PRIME). Recombinant plasmids were verified by colony PCR (data not shown) prior to DNA sequencing by 1st BASE DNA Sequencing Services (Malaysia). All fragments were fully sequenced from both strands of the DNA.

Sequence alignment and phylogenetic analysis

All bidirectional sequences were assembled using Geneious version 5.1.7 (Biomatters Limited; www.geneious.com). The mitochondrial nucleotide sequences of the COI barcoding region were then visually aligned before phylogenetic analyses were performed.

The nuclear nucleotide sequences of the ITS1 region containing the end of 18S rRNA, ITS1, and the beginning of 5.8S rRNA were not as easily aligned, partly due to the quality of the sequences obtained from the direct sequencing method using for primers Betta_ITS1_F1 and Betta_ITS1_R1 (see Table 2) and also due to ITS1 having three relatively large indels and repeating regions. Since ITS1 was not commonly used for vertebrates, appropriate sequencing primers were not available. We thus designed two new forward primers (Betta_ITS_F2 and Betta_ITS_F3) and two new reverse primers (Betta_ITS_R3 and Betta_ITS_R4), corrected the Betta_ITS1_R1 primer to be Betta_ITS1_R2 (as shown in Table 2), and used the DNA cloning method to obtain better quality sequences. Among the three pairs of primers used in the DNA cloning method, the pair Betta_ITS_F2 and Betta_ITS_R4 yielded the most consistently complete sequences. As a result, a higher number of good quality ITS1 sequences were obtained so that the end of 18S rRNA and the beginning of 5.8S rRNA could be ascertained. Because these two regions were totally conserved, reliable ITS1 sequences could be extracted. All the acceptable ITS1 sequences were then aligned using MUSCLE (Edgar, 2004) and further adjusted visually. Phylogenetic trees were reconstructed from the aligned COI sequences, ITS1 sequences, and the combined sequences of the two using partitioned Bayesian inference with reversible jump Metropolis-coupled Markov chain Monte Carlo (rjMC3) (Huelsenbeck and Ronquist, 2001, Huelsenbeck et al., 2004) implemented in MrBayes version 3.2 (Ronquist et al., 2012). For COI and ITS, phylogenetic analyses were performed on both the full datasets and the partial datasets containing the taxa in the combined dataset (46 members).

In this study, the nuclear nucleotide sequences of a region of the RAG1 gene of some specimens above were also aligned. However, the differences among these sequences by themselves were not as discriminative as COI, ITS1 and combined COI + ITS1 sequences sufficiently to further differentiate the species and cryptic species. Naturally then, the inclusion of RAG1 as the combined RAG1 + COI + ITS1 sequence analyses did not change the topology of the inferred tree (data not shown).

Results

DNA sequences and their alignments

The COI sequences of all bettas here were 652 base pairs long starting at the third position in the same codon. There were 210 variant sites, 178 of which were informative for species identification. A majority of third positions (157 out of 218) varied, compared to only 37 first positions and 16 second positions out of 217 sites. One hundred sixty four sites contained exactly two possible bases, which yielded a rough estimate of the transition–transversion ratio at 9.25 (data not shown). The ITS1 sequences varied from 363 to 440 base pairs long and their alignment yielded a sequence of 470 base pairs long. Final alignment with the 291-base-pair ITS1 sequence of the Betta ocellata yielded a 480-base-pair sequence.

Phylogenetic trees

For reconstructing the phylogenetic tree from the combined (COI and ITS1) dataset, the latter was divided into four partitions: COI codons (each with three base positions) took up three partitions and the remaining one was for ITS1. Using the Bayesian information criterion (BIC) (Schwarz, 1978), + G was justifiable only for the first codon partition of COI in both the COI and combined datasets. The phylogenetic tree of the combined dataset, which was similar to that of the COI dataset (data not shown), is shown in Fig. 1. The number at each node represents the posterior probability from the Bayesian inference using reversible jump Metropolis-coupled Markov chain Monte Carlo. The posterior probabilities of all the (potential) speciation nodes in the tree reconstructed from the combined dataset were very high (≥ 0.97) indicating that all those bipartitions were extremely likely given the sequences. On the contrary, some of the corresponding posterior probabilities in the tree reconstructed from COI sequences (data not shown) were only fair, the most interesting being the one corresponding to the bipartition between the B. smaragdina group and the rest (1 in the combined tree versus 0.70 in the COI tree), indicating that the COI dataset alone might not be sufficient to reconstruct a highly reliable tree. In addition, the posterior probabilities of all the most recent speciation nodes in the tree reconstructed from the combined dataset were 1, supporting the species status of all the species and cryptic species and lending credence to the use of the combined dataset in the reconstruction of the phylogenetic tree.

The topologies of the trees reconstructed using the COI dataset and the one using the combined dataset were mostly identical (data not shown) while that based on the ITS1 dataset was noticeably different with much shorter branch lengths and much lower posterior probabilities (data not shown) indicating that COI pre-dominated over ITS1 in the combined dataset. Because base differences in COI result from substitutions while most of those in ITS1 result from indels, predominance of COI is to be expected. Even with the differences between the COI and ITS1 trees, the combined tree (see Fig. 1), upon which further discussions are based, had significantly higher posterior probabilities than the corresponding ones in the COI tree on average. However, this might be the effect of the longer sequences.

Species divergences

The intraspecific K2P divergences (Kimura, 1980) for the five known species of Thai bubble-nesting bettas (B. imbellis, B. mahachaiensis, B. siamorientalis, B. smaragdina, and B. splendens) and those of the B. spp. (cf. smaragdina) 1, 3, and 4 are shown in Table 3. B. mahachaiensis had the lowest intraspecific divergence of 0.03% (due to the one small and rather unique habitat). The intraspecific divergences among all the B. smaragdina group members ranged from a lower difference of 0.11% in B. sp. (cf. smaragdina) 1 to a much higher difference of 1.42% in type-locality B. smaragdina; the latter number was also the highest among the intraspecific divergences of all the taxa. The average intraspecific divergence among all the betta taxa discussed here was 0.41% which was lower than that among all the B. smaragdina group members (0.66%). In fact, the latter number, which should be very low especially after subdividing the group members into 4 clusters, was still much higher than the intraspecific divergences of all the non-smaragdina bettas, indicating that our subdivision might be, if anything, too conservative and suggesting the possibility of even more additional species.

Table 3

Intraspecific K2P divergences for the cytochrome oxidase c subunit 1 (COI) gene.

Species	K2P distances	N
Betta siamorientalis	0.0012	37
Betta imbellis	0.0022	36
Betta splendens	0.0023	189
Betta mahachaiensis	0.0003	79
Betta sp. (cf. smaragdina) 1	0.0011	44
Betta sp. (cf. smaragdina) 4	0.0081	6
Betta sp. (cf. smaragdina) 3	0.0032	9
Betta smaragdina (type-locality)	0.0142	144
Betta stiktos	0.0040	7
Betta ocellata	n/c	1
Total		552

The interspecific K2P divergences among these taxa are shown in Table 4. The type-locality B. smaragdina had percentage divergences that differed from B. spp. (cf. smaragdina) 1, 3, and 4 by 13.68%, 8.84%, and 8.21% respectively, all of which are greater than the cut-off points of 1%–2% suggested in the literature (see, for example, Schwarz, 1978 and Steinke et al., 2009), again indicating the possibility of additional new species. Among the four members of B. smaragdina group, B. sp. (cf. smaragdina) 1 (from northernmost part of the northeast region—Bueng Kan province) was the most distant, diverging by about 14% from both the type-locality B. smaragdina (from Korat basin and one half of Sakon Nakhon basin) and B. sp. (cf. smaragdina) 3 (from northern Lao PDR), and about 13% from B. sp. (cf. smaragdina) 4 (the other half of Sakon Nakhon basin). The average interspecific divergence among these B. smaragdina group members was 10.82%, which was about 16 times the corresponding average intraspecific divergence, well within the normal range of the same ratio for ornamental fishes (Steinke et al., 2009). Betta splendens had the percentage divergences different from B. mahachaiensis, B. siamorientalis, B. imbellis, B. smaragdina, and B. stiktos by 10.67%, 10.73%, 11.46%, 16.68%, and 17.03% respectively.

Table 4

Pairwise interspecific K2P distances for the cytochrome oxidase c subunit 1 (COI) gene.

	B. siamorientalis	B. imbellis	B. splendens	B. mahachaiensis	B. sp. (cf. smaragdina) 1	B. sp. (cf. smaragdina) 4	B. sp. (cf. smaragdina) 3	B. smaragdina (type–locality)	B. stiktos
B. siamorientalis	–
B. imbellis	0.0133	–
B. splendens	0.1073	0.1146	–
B. mahachaiensis	0.1099	0.1081	0.1067	–
B. sp. (cf. smaragdina) 1	0.1381	0.1399	0.1449	0.1524	–
B. sp. (cf. smaragdina) 4	0.1451	0.1389	0.1789	0.1670	0.1268	–
B. sp. (cf. smaragdina) 3	0.1642	0.1573	0.1694	0.1824	0.1429	0.0722	–
B. smaragdina (type–locality)	0.1506	0.1463	0.1668	0.1607	0.1368	0.0821	0.0884	–
B. stiktos	0.1453	0.1411	0.1703	0.1562	0.1421	0.0945	0.0933	0.0213	–
B. ocellata	0.2397	0.2422	0.2385	0.2350	0.2389	0.2300	0.2437	0.2340	0.2327

In this study, even though there were some possibly hybrid fish, none of them were from the wild. We have found discordant COI and ITS1 in the same individuals of B. imbellis, B. mahachaiensis, and B. siamorientalis bought from the Bangkok Chatuchak aquarium fish market, all of which the COI sequences were identified as B. splendens.

The results in this study (see Fig. 1) also show for the first time the positions of B. stiktos and B. siamorientalis relative to other species in the phylogenetic tree. This is congruent with the morphological minor differences between B. stiktos and type-locality B. smaragdina and between B. siamorientalis and B. imbellis. There are also only 2.13% and 1.33% sequence divergence differences between members of each pair.

Discussion

There was compelling evidence from COI sequences and subsequent phylogenetic analyses in our 2010 paper (Sriwattanarothai et al., 2010) that two B. spp. (cf. smaragdina) should enjoy a distinct species apart from the type-locality B. smaragdina. However, we decided to call both cryptic species because of the absence of the supporting evidence. In the present report we will maintain this position by referring to all newly identified B. smaragdina group members as cryptic species except, naturally, for the type-locality one.

Reconstruction of phylogenetic trees

In this study, we reconstructed phylogenetic trees using partitioned Bayesian interference with rjMC3 (Huelsenbeck and Ronquist, 2001, Huelsenbeck et al., 2004). In traditional likelihood-based inference, the best nucleotide substitution model has to be selected from a limited number of models (22 models in the relatively extensive jModelTest version 0.1.1 (Posada, 2008) before it is used in phylogenetic tree reconstruction. The rjMC3, on the other hand, allows all the 203 time-reversible substitution models to be parts of the posterior distribution. That is, in addition to the topology, branch lengths, and all other parameters, the substitution models and their substitution rates, allowed to vary across partitions, were also sampled from the posterior distribution. The mutation rates were also allowed to vary across partitions. In order to determine whether there was rate variation across sites (+ G) or a proportion of invariable sites (+ I) for each partition, jModelTest was used to gauge the relative likelihoods of the maximum-likelihood (ML) trees using substitution models with and without + G and + I. It turned out that + I was not required in all the partitions while + G was probably needed. Thus, Bayesian inference was carried out for each of the three datasets with rate variation across sites for all the partitions and the sampled shape parameters were inspected using Tracer version 1.5 (Rambaut and Drummond, 2003). Only those partitions whose shape parameters concentrated at low positive values (high values of shape parameters correspond to low variation) would be the candidates to retain + G. The inference was repeated with and without + G for each candidate partition in each dataset and the likelihoods were compared using the BIC to decide whether the extra shape parameter was justifiable. The result was used to run rjMC3 for each of the three datasets with two runs, each with one cold chain and three heated chains with the incrementing heating temperature of 0.1. One million generations were performed for each run with the sample frequency of 500 and 250,000 generations were discarded as the burn-in period. All the indicators in MrBayes were inspected to ensure the convergence and the samples from both runs were combined to reconstruct the phylogenetic trees.

Biodiversity of the B. smaragdina group members in northeast Thailand

In 2010, based on COI and 16S RNA results, we reported that B. smaragdina comprised a greater diversity than that known up to then (Sriwattanarothai et al., 2010). There were at least two more cryptic species in addition to the one known as the type-locality from the province of Nong Khai (changed from Korat as reported by Sriwattanarothai et al., 2010, according to reliable information from Horst Linke). This change in the provincial origin of the type-locality does not affect our results in that the Korat and Nong Khai fish are the same. The present study represents a more comprehensive coverage of B. smaragdina species over the length and width of northeast Thailand, i.e., from the northernmost provinces to the southernmost provinces. Our tools were mainly comparing short pieces of DNA from the mitochondria (COI) and that from the nucleus (ITS1) to distinguish among the apparently B. smaragdina-like bubble-nesting fighting fish. As a result of our present work, the relative phylogenetic positions of B. splendens, B. imbellis, and type-locality B. smaragdina remain the same as in our 2010 publication (Sriwattanarothai et al., 2010, Rüber et al., 2004a). Regarding B. sp. Mahachai, the present findings also support the fish's position relative to the three mentioned above as reported previously (Sriwattanarothai et al., 2010).

Phylogenetic analyses based on COI and ITS1 data separately (data not shown) gave 3 and 2 cryptic species of the B. smaragdina group respectively, albeit their positions on the phylogenetic trees being slightly different. Knowing that the ITS1 is generally more conservative than COI in percent base change for a typical organism, combining COI and ITS1 data yielded a phylogenetic tree showing B. sp. (cf. smaragdina) 1 well separated in one branch, whereas the other three members of the B. smaragdina group, including the hitherto type-locality, were more clustered together (see Fig. 1). One interesting point emerging from the work was that the type-locality B. smaragdina genetically resembles the B. stiktos (from Stung Treng, Cambodia) the most, making it worthwhile to explore their differences or similarities further. In other words, is B. stiktos different enough from the best-known B. smaragdina (type-locality) to be another species? Because the interspecific K2P divergence between them and the branch length leading to the B. stiktos clade were the lowest among all such values within the B. smaragdina group (see Table 4 and Fig. 1), if these values of divergence and branch length were good enough to differentiate B. stiktos from the type-locality B. smaragdina, the other B. smaragdina group members should also enjoy the species status as well. Obviously, other additional characters, DNA and anatomical ones, should also enhance their differentiation.

The most significant agreement between the ITS1 and the combined trees was the monophyly of the B. smaragdina group (including B. stiktos) which was uncertain as of Sriwattanarothai et al. (2010) and only moderately supported in the COI tree (posterior probability = 0.70, data not shown). However, the amino acid sequences translated from all COI sequences and the ITS1 alignment clearly supported monophyly of the B. smaragdina group (results not shown). The contribution of ITS1 in this regard was rather significant since it was consistent with the morphological evidence. The combined tree also supported the species statuses of B. mahachaiensis (Kowasupat et al., 2012a), B. siamorientalis (Kowasupat et al., 2012b) and B. stiktos (Tan and Ng, 2005), with the posterior probabilities of 1 for all three clades. Given these species statuses, it is likely that there are three new cryptic species (all also with posterior probabilities of 1) within the B. smaragdina group, whose statuses should nevertheless be further ascertained. Our alignment also revealed that the only accepted B. smaragdina up to then in Sriwattanarothai et al. (2010) is not the type-locality (neither from Nong Khai nor Korat) as first believed. The latter fish is found by the present study to be our B. sp. (cf. smaragdina) 4. From their DNA sequences identity, the actual type-locality reported here is the same fish as that referred to as B. sp. (cf. smaragdina) 2 by Sriwattanarothai et al. (2010). Another clade with a rather high uncertainty by Sriwattanarothai et al. (2010) was the one containing B. splendens and B. mahachaiensis with the posterior probability of only 0.63 (B. imbellis was very close to this clade). Now this clade was highly supported in our analysis with the posterior probability of 0.99.

The position of B. siamorientalis in the phylogenetic tree in Fig. 1 indicates that B. imbellis is paraphyletic with respect to B. siamorientalis even though the latter is monophyletic. This fact should not diminish the species status of both species even though they could not satisfy the criterion of reciprocal monophyly (CRM). To indicate this subtlety in the tree, Kizirian and Donnelly (2004) suggested that the two taxa together could be referred to as B. imbellis complex which is a monophyletic taxon and the binomial B. imbellis should be eliminated. According to this view, the species status of B. siamorientalis is stronger than that of B. imbellis. However, we see no real need in eliminating the binomial B. imbellis, especially given that its intraspecific K2P divergence was only 0.22%. The same can be said about the position of B. stiktos with respect to the type-locality B. smaragdina. However, because the intraspecific K2P divergence of the latter was very high (1.42%) and the subtree of the B. smaragdina type-locality clade almost totally covered that of the B. stiktos clade, the name B. smaragdina complex should be better at portraying this relationship and suggesting that there might be other cryptic species apart from the three in this study.

In this work we have constantly borne in mind the controversies surrounding the use of just one short piece of mitochondrial DNA for species identification (Ward et al., 2009, Zemlak et al., 2009, Dasmahapatra et al., 2010, Taylor and Harris, 2012) and tried to make our work as integrative as possible (Pires and Marinoni, 2010). We have thus used both a mitochondrial and a nuclear DNA for species identification. However, one important aspect that needs addressing is the differentiating morphological data which so far have not been forthcoming. The most we can say at this stage about morphological difference is that most (68 out of 73) specimens of the B. sp. (cf. smaragdina) 1 have a spider web pattern across the caudal rays versus spots between rays of other B. smaragdinas; thus we cannot cite this trait as truly belonging to this particular cryptic species. Also the B. sp. (cf. smaragdina) 1 displays a unique, among B. smaragdinas, pelvic fin flickering when displaying aggression (Sriwattanarothai et al., 2010). Incidentally, we wish to show by the following statements that short DNA sequences and knowledge of local situations may be of value in deciding species status. Earlier we mentioned that most bettas sold for ornamental or fighting purposes and some of the caught ones, despite their various differences in external appearance, have been found by COI barcoding to have extremely low divergence (0.23%, see Table 3) showing all of them to be B. splendens. Some of these undesired fighter fish, regularly released by their owners into “natural waters”, which are caught by us near populated areas also have a variety of shapes and colors on the fins and body. Without the knowledge of how they have found their way into natural waters and breed there, and their identification by DNA sequences as described above, one would think that these fish belong to species other than B. splendens. Thus external appearance of a mature fish may sometimes be misleading. These released B. splendens fighters generally, like the traded fighters, have stouter bodies and relatively smaller unpaired fins than the wild ones caught in remote areas.

The most consistent external feature of all the B. smaragdina group members is the dark blue to blue iridescent scales covering most of the area of the operculum, all with B. smaragdina body and fins. We have tried to distinguish these four B. smaragdina group members by other obvious external features and behaviors without too much success so far. It is thus worth pursuing further other traits such as internal anatomical structures, external patterns and colors, and behavioral displays. At this juncture we would like to remark that the B. smaragdina group and the B. splendens group may be distinguished by the iridescent yellowish-green to bluish-green scales covering almost all areas of the operculum in the former and the double parallel vertical bars of red or green to bluish-green color in the latter. The primary node in the phylogenetic tree forking the two groups into two distinct branches is also worth noticing. The nucleic acid base transition-transversion patterns and the amino acid substitution characteristics of the two groups are also very different (data not shown).

It should be interesting and informative to find explanations in terms of biogeographical past, recent, and present for the disparate biodiversity among the B. smaragdina group members in different parts of the northeast of Thailand (Sinsakul et al., 2002, Dudgeon, 2008). Also the very low divergence among B. splendens should also be explained, especially given great expanses from the north to central to upper south of Thailand occupied by the fish. B. imbellis is located in a smaller area of the lower southern peninsula which is near the equator. B. mahachaiensis and B. siamorientalis occupy habitats of smaller size still in central and eastern Thailand.

Using ITS1 to identify species of bettas

The alignment of ITS1 sequences revealed that there were very few differences among the ITS1 of all the bettas in this study. As a result, the phylogenetic tree reconstructed from the ITS1 dataset contained very short branch lengths with very low posterior probabilities. However, in addition to the main contribution of ITS1 in supporting the monophyly of the B. smaragdina group mentioned earlier, the DNA sequence can also be used to identify most species of bubble-nesting bettas in Thailand. Table 5 presents the genetic markers that can be used to identify these species. The doubly underlined bases are unique to that particular species. We can see that the genetic marker for B. mahachaiensis contains five unique bases while other markers contain at most one. The markers for B. imbellis and B. sp. (cf. smaragdina) 3 contain no unique bases but the combinations of singly underlined bases are enough to identify these species. Notice that there is no unique base marker for B. smaragdina (type-locality) and B. sp. (cf. smaragdina) 4 which means that if an ITS1 sequence of a bubble-nesting betta in this study contains none of the markers in Table 5, the ITS1 sequence belongs to either B. smaragdina (type-locality) or B. sp. (cf. smaragdina) 4.

Table 5

Genetic markers for the ITS1 sequences of the bettas in this study.

The doubly underlined bases and the combinations of singly underlined bases are unique to the species.

Limitations of this work

There have been criticisms of species identification based on COI barcoding despite its success in doing so in many instances (Rüber et al., 2004b, Rambaut and Drummond, 2003). However, the combined use of mitochondrial and nuclear short pieces, e.g., COI and ITS1 in this study, should make identification stand on a firmer ground. Our results should be enhanced further by the use of RAD-seq tags methodology which yields multiple-loci DNA sequences of the genome and thus is better at distinguishing closely related species (Emerson et al., 2010). Naturally, if and when morphologically differentiating criteria are available, identification will then be more generally acceptable.