Molecular systematics and morphological variation in the Mesopotamian spiny eel Mastacembelus mastacembelus (Teleostei: Mastacembelidae)

Abstract

The Mesopotamian spiny eel, Mastacembelus mastacembelus, an endemic freshwater fish restricted to the Euphrates and Tigris river drainages of Mesopotamian, and its related taxa in Afrotropica and South Asia present suitable elements for studies of the evolutionary history of mastacembelids across these heterogeneous areas. To understand and clarify the geographic variation of the Mesopotamian spiny eel from the Persian Gulf basin, we investigated the molecular and morphological characteristics of individuals collected from its whole distribution range in the Tigris (Karkheh, Karoon and Zab), Zohreh and Persis (Mond and Helleh) main river drainages in Iran. The morphological analysis revealed a separation of the studied populations; however, the genetic distances were relatively small and the molecular tree did not present a distinct branching pattern for the main river drainages based on the COI mitochondrial marker. Six haplotypes were identified: Helleh, Mond, Karkheh, Karoon, Zab and Zohreh tributaries each showed a separate haplotype. The COI marker indicates that the Mesopotamian eel has greater affinity with the African species than with the Southeast Asian species. Estimation of divergence times indicates the Mesopotamian spiny eel apparently diverged about 1 million years ago (Mya) (95% highest posterior density [HPD]: from 390,000 years ago to 2.1 Mya) in different river drainages of the Persian Gulf basin.

Keywords:

• Freshwater spiny eels
• Mastacembelids
• Phylogeny
• Mesopotamia

Introduction

Freshwater spiny eels/mastacembelids (Teleostei: Synbranchiformes: Mastacembelidae) are a predominantly riverine family of freshwater fish distributed in the Old World from the tropical parts of Africa to Southeast Asia, containing the three genera Macrognathus, Sinobdella and Mastacembelus and 93 species (Coad 2015; Hastings et al. 2015; Fricke et al. 2021).

To explore the diversity, monophyly and colonization history of Mastacembelus eels, the first molecular phylogenetic analysis of this group was performed by Brown et al. (2010), using a multigene dataset of mitochondrial and nuclear genes applying several methods. The study highlighted the Africa–Asia split of mastacembelid eels at ~19 Mya, which is considerably younger than the split between the two continents, suggesting a dispersal scenario for their current distribution. Also based on the ancestral range estimation, it has been demonstrated that mastacembelids and their constituent genera originated and diverged preliminary in Asia, with Mastacembelus suggested to have reached Africa (the Afrotropical Mastacembelus lineage) shortly before its divergence from the Mesopotamian/Middle East spiny eel, Mastacembelus mastacembelus, between 27.5 and 12 Mya (Lavoué 2019). However, Day et al. (2017) concluded the monophyly of an African Mastacembelus radiation plus the Mesopotamian spiny eel, M. mastacembelus (with uncertain position), and highlighted the origination of the spiny eel in Asia and its colonization of Africa c. 15.4 Mya (95% highest posterior density [HPD]: 23.9–8.8 Mya), presumably through the Middle East.

The Mesopotamian spiny eel is interesting in term of phylogeography, and it has not been well studied in the Middle East. The morphological status of the Mesopotamian spiny eel populations from Karakay Reservoir, Tohma Stream and Tigris River were investigated by Çakmak and Alp (2010) using morphometric and meristic traits. Tutar (2015) compared three populations of the species from the Turkish part of the Tigris River drainage based on the combined mitochondrial DNA (12S rRNA, 16S rRNA, tRNAPhe and tRNAVal genes). However, morpho-molecular variations of the species have not been investigated among the Iranian populations so far, although Iran includes a large part of its distribution. Iranian populations of this species inhabit different river tributaries and habitats with different ecological conditions that have prompted us to initiate a more detailed investigation, both morphologically and molecularly. Hence, this study aims to understand and clarify the geographic variations of the Mesopotamian spiny eel from the Persian Gulf basin, using integrated molecular and morphological characteristics of the individuals collected from its whole distribution range in the Iranian part of the Persian Gulf basin.

Materials and methods

Sampling

A total of 151 specimens of Mesopotamian spiny eel were collected from its whole distribution range in the Tigris, Zohreh and Persis (Mond and Helleh) river drainages (Persian Gulf basin) in Iran using an electrofishing device (ELEMAX, SHX 2000) during fieldwork in 2016 and 2017 (Figure 1; Table I). An attempt to obtain samples in the Kor River basin was not successful (see Discussion). After anesthesia with 1% clove oil solution, the pectoral fin (or, from some specimens, tissue from below the dorsal fin) on the right side of each specimen was removed and preserved in 96% ethanol for the molecular analyses. A total of 126 adult individuals were fixed in 10% formaldehyde for morphological studies. Voucher specimens were deposited in the Zoological Museum, Collection of Biology Department, Shiraz University (ZM-CBSU).

Table I. Sampling sites of Mesopotamian spiny eel and details of specimens

Sampling site	N	River drainage	Accession no.
Choman River* (n = 4, M1607, M1608, M1609, M1610)	9	Tigris	ON159366, ON159367, ON159368, ON159369
Gomasiab River* (n = 1, M1641)	7	Tigris	ON159382
Kashkan River* (n = 2, M1616, M1617)	23	Tigris	ON159374, ON159375
Saimareh River* (n = 2, M1614, M1615)	11	Tigris	ON159372, ON159373
Chardavol* (n = 1, M1638)	0	Tigris	ON159381
Sazbon River (Ilam)	14	Tigris
Sabzab River (Andimeshk)* (n = 4, M1581, M1582, M1611, M1612)	7	Tigris	ON159362, ON159363, ON159370, ON159371
Zohreh River (Sardasht)	6	Zohreh
Kheirabad River* (n = 2, M1541, M1542)	12	Zohreh	ON159360, ON159361
Fahlian River* (n = 2, M1583, M1584)	12	Zohreh	ON159364, ON159365
Dadin spring* (n = 2, M1635, M1637)	9	Persis	ON159379, ON159380
Dalaki River* (n = 2, M1539, M1540)	7	Persis	ON159358, ON159359
Shahpur River	1	Persis
Simakan spring* (n = 3, M1631, M1632, M1633)	8	Persis	ON159376, ON159377, ON159378
Total n = 25	N = 126

N = number of specimens used for morphological analyses (collected in this study).

n = number of specimens used for molecular analyses.

*: Populations used for molecular analyses.

Figure 1. Sampling localities of the Mesopotamian spiny eel from Iran. Red circles indicate sampling locations defined in Table I.

Molecular analyses

DNA extraction, polymerase chain reaction and sequencing

A total of seven specimens of the Mesopotamian spiny eel from the Persis (Mond and Helleh) river drainages, four specimens from the Zohreh River drainage and 14 specimens from the Tigris River drainage were selected for molecular analyses (Table I). Total genomic DNA was extracted from the preserved tissues using the standard salt extraction method (Bruford et al. 1992). We amplified the DNA barcode region of the mitochondrial COI marker using a polymerase chain reaction (PCR) procedure with the universal primers FishF1 and FishR1 (Sharina & Kartavtsev 2010) and the following amplification protocol: initial denaturation at 94°C for 5 min; 35 amplification reaction cycles consisting of denaturation at 94°C for 60s, annealing at 55°C for 60s, and extension at 72°C for 60s per cycle; and a final extension cycle at 72°C for 5 min. PCR products were evaluated with a 100 bp ladder on 1% agarose gel electrophoresis. Sequencing was carried out using a forward primer by Macrogen Service Centre (Seoul, South Korea). The obtained sequences were deposited in the national center for biotechnology information (NCBI) GenBank under specific accession numbers (see Table I).

Phylogenetic analyses

In addition to the samples sequenced here, sequences of one Mastacembelus species from southeastern Asia (LT577001.1 M. armatus) and another from Africa (LT576985.1 M. taiaensis) were retrieved from NCBI with M. armatus as the outgroup. We also included the only COI sequence of Mesopotamian spiny eel from the Middle East in the NCBI (with the accession number LT577004, from Turkey). The sequences were edited and aligned using BioEdit v. 7.0.0 (Hall 1999) and checked for the presence of unexpected stop codons using Mega 7 software (Kumar et al. 2016). The model suggested for the COI was determined using jModelTest 2.1.10 and GTR+G + I selected according to the Akaike information criterion. Bayesian inference (BI) was done using MrBayes 3.1.2 (Ronquist & Huelsenbeck 2003) with 10 million generations with four Markov chain Monte Carlo sampling and with a sampling frequency of 100; 25% of the trees were discarded as burn-in. A maximum likelihood (ML) (Felsenstein 1985) gene tree was inferred using RAxML 7.2.5 with 10,000 bootstrap replicates (Stamatakis 2006) and the GTR+G + I model to examine the robustness of the Bayesian results. Mean genetic distances among populations were calculated using the Kimura two-parameter model (K2P) implemented in Mega v. 7.0, and this software was also used to extract variable nucleotide sites from the sequence alignment.

Estimation of divergence times

Divergence times of Iranian populations were estimated using the time tree provided by Day et al. (2017). The sequences examined here (belonging to other species) were retrieved from the NCBI, mostly from Day et al.’s (2017) study. We included African and Asian species as monophyletic, and then three nodes from the previously published time tree were considered to calibrate the tree: mean age of the family at 29 Mya, mean age for divergence of Mastacembelus at 21 Mya and its divergence in Africa at 15.5 Mya. Rates among branches were assumed to follow a lognormal distribution. We specified a normal prior distribution for the mean age of the nodes with a standard deviation (SD) of 3 for the family, 2 for the genus and 1.5 for the African clade using BEAST v. 1.7.2 software. The analysis was run for 100 million generations under the GTR+G + I model, with empirical base frequencies, Yule process and a random starting tree with parameters sampled every 10,000 steps, and Tracer v. 1.4 was used to assess convergence and effective sample size (ESS) following a preburn-in of 10%. Finally, we generated maximum clade credibility trees with mean node heights using TreeAnnotator v. 1.7.4, and the time tree was visualized in FigTree v. 1.4.4. We considered Macrognathus zebrinus the outgroup and included one sequence per haplotype from our 25 new sequences.

Morphological analyses

The morphological analyses were carried out using 126 specimens covering the whole distribution range of the Mesopotamian spiny eel in the Tigris River drainage (n = 71), the Mond and Helleh river drainages (n = 25) and the Zohreh River drainage (n = 30), which all drain into the Persian Gulf. Morphometric characters were measured with digital calipers and recorded to 0.01 mm. Methods for counts and measurements follow Plamoottil and Abraham (2013) and include the following 30 characters: total length (TL), standard length (SL), head length (HL), head depth (HD), head width (HW), mouth width (MW), eye diameter (ED), interorbital distance (InOD), preorbital distance (Prod), postorbital distance (POD), predorsal distance (Prdd), preanal distance (Prad), length of dorsal fin (Ldf), length of anal fin (Laf), length of pectoral fin (Lpcf), length of caudal fin (LCF), distance between pectoral and anal fin (Dpcaf), maximum body depth (MaxBD), minimum body depth (MinBD), length of base of pectoral fin (BLPF), length of spinous anal fin ray (LSAF), length of base of soft anal fin ray (LBSAF), length of base of soft dorsal fin ray (LBSODF), length of soft dorsal fin ray (LSDF), height of soft dorsal fin (Hsdf), height of soft anal fin (HSAF), height of anal fin spine (HAS), depth of dorsal fin (Ddf), depth of anal fin (DAF) and length of base of spinous dorsal fin (LBSpF). In addition, five meristic characters were counted under a stereomicroscope as follows: the numbers of dorsal fin spines (DFS), anal fin spines (AFS), pectoral fin rays (PFR), anal fin rays (AFR) and dorsal fin rays (DFR).

To eliminate any size effects, the morphometric data were standardized according to the equation (Mmp/SL)*100, where Mmp: measured morphometric parameter, and SL: standard length, except for head characters (HD, HW, ED, Prod, POD, InOD) that divide into HL (head length). Meristic characters are independent of fish size and were used as raw data. Canonical discriminant analysis (CDA) was carried out to determine the success of classification based on the morphometric measurements and meristic counts separately, using IBM SPSS v. 26. Principal component analysis (PCA) was also performed, to assess individual relationships within and between populations, based on the morphometric and meristic characters separately. To examine the suitability of the data for PCA, the Kaiser–Meyer–Olkin (KMO) measure of sampling adequacy and Bartlett’s test of sphericity were performed (KMO: 0.72; P < 0.001) using IBM SPSS v. 26. Then, PCA analysis was performed to visualize individual relationships using PAST 3.

Results

Phylogenetic analyses

The COI sequences including 536 nucleotides were obtained from 25 specimens of the Mesopotamian spiny eel, of which 526 nucleotides were constant and 10 were parsimony informative. Two different phylogenetic analyses (BI, ML) for Mastacembelus spp. produced nearly identical topologies; thus, only the topology from the ML analysis is presented, along with the statistical values from the ML and BI algorithms (Figure 2).

Figure 2. Phylogenetic estimates for Mastacembelus mastacembelus based on the COI mitochondrial marker. The topology is taken from the ML, and numbers above branches are bootstrap values/posterior probabilities. ML: maximum likelihood.

Some molecular variations were recovered in the Iranian Mesopotamian spiny eel. Table II lists the 10 variable nucleotide sites found in the COI barcode region. Specimens of Zohreh and Persis are more similar to each other than to the Tigris, and Tigris specimens including Karkheh, Karoon and Zab tributaries showed more molecular variations (Table II). Specimens of Karkheh tributary (Gomasiab, Seimareh, Kashkan, Chardavol) showed a separation from the rest of the samples (Zohreh, Mond, Zab, Helleh, Karoon tributaries) (Figure 2). Specimens of the three main river systems (Tigris, Zohreh, Persis) were not nested in three separate groups; however, six haplotypes were identified in total and Helleh, Mond, Karkheh, Karoon, Zab and Zohreh tributaries each showed a separate haplotype.

Table II. Variable sites of Mesopotamian spiny eel relative to KY609156.1 Mastacembelus armatus complete mitochondrial genome

We found genetic distances (K2P) of between 0.2% and 1.3% for the COI among the tributaries (Table III). The maximum genetic distance was observed between specimens of Kharkeh tributary and Choman. The genetic distances among Tigris, Zohreh and Persis populations varied from 0.2% to 0.7%.

Table III. Average genetic distance between pairs of tributaries based on Kimura two-parameter (K2P). Bold values report maximum and minimum distance values between tributaries

Tributary	1	2	3	4	5	6
Dalaki, Dadin (Helleh)	0
Simakan (Mond)	0.004	0
Kheirabad, Fahlian (Zohreh)	0.002	0.002	0
Kharheh	0.010	0.010	0.008	0
Andimeshk (Karoon)	0.006	0.006	0.004	0.011	0
Choman (Zab)	0.008	0.004	0.006	0.013	0.010	0

Divergence time

Inclusion of Macrognathus as outgroup and Asian and African species as monophyletic clades suggested M. mastacembelus is the sister of the African mastacembelids (Figure 3). The diversification of Mastacembelus in the Middle East is relatively recent compared with diversification in Africa, and the Mesopotamian spiny eel diverged apparently about 1 Mya (95% HPD: from 390,000 years ago to 2.1 Mya) in different river drainages of the Persian Gulf basin.

Figure 3. Divergence time estimates of the major cladogenetic events for species of Mastacembelus based on the COI gene using BEAST analysis. The blue bars illustrate the extent of the 95% HPD credibility intervals for each divergence time. HPD: highest posterior density.

Morphological analyses

All traits were subjected to the CDA, which showed the random assignment of individuals to their basin was 85.7% based on the morphometric characters (Figure 4) and 57.9% for meristic characters (Figure 5). The classifications were 85.9% for Tigris, 92% for Persis and 80% for Zohreh populations based on the morphometric characters and 54.9% for Tigris, 70% for Zohreh and 68% for Persis based on the meristic characters. Discriminant functions (DF) 1 and 2 together explained 100% of the variance, with DF1 accounting for 68.7% and DF2 31.3% of all variance for the morphometric characters; and with DF1 explaining 82.2% and DF2 17.8% of all variance for the meristic characters. The test of functions for DF 1 and 2 was significant for both morphometric and meristic data (Wilks’s lambda, P < 0.001).

Figure 4. Scatter plot of the canonical discriminant analysis based on morphometric variables of the Mesopotamian spiny eel from three river drainages of Iran.

Figure 5. Scatter plot of the canonical discriminant analysis based on meristic variables of the Mesopotamian spiny eel from three river drainages of Iran.

The PCA carried out on the morphometric measurements extracted 13 factors with eigenvalues > 1, explaining 94.39% of the variance. The first principal component (PC) described 37.47% of the variation and the second 14.96%. Therefore, the first PC was the most important component contributing to separation among the basins. The most significant loadings on PC 1 were HW and MW and on PC 2 they were HD and HW, all from head characters. With regard to factor loading, values above 0.30 are considered significant (Lombarte et al. 2012); in the present study significant factors on PC1 were HW and MW and significant factors on PC2 were HD, HW and LBSpF. A visual inspection of the plotted PC1 and PC2 scores showed a high degree of overlap among the three basins (Figure 6).

Figure 6. Principal components analysis (PCA) of morphometric characters of Mesopotamian spiny eel from three river drainages of Iran. Plus symbol: Tigris; multiplication symbol: Zohreh; diamond: Persis.

PCA on the meristic measurements extracted three factors with eigenvalues > 1, explaining 97.85% of the variance. The first PC described 88.06% of the variation and the second 5.97%. The most significant loadings on PC 1 were AFR and DFR, and the most significant loading on PC 2 was DFS (Figure 7).

Figure 7. Principal components analysis (PCA) of meristic characters of Mesopotamian spiny eel from three River drainages of Iran. Plus symbol: Tigris; multiplication symbol: Zohreh; diamond: Persis.

Discussion

The Middle Eastern species M. mastacembelus is restricted to the inland waters of the Middle East, occurring in a wide range of isolated water bodies in the Euphrates and Tigris river drainages/Mesopotamia, in Iran, Turkey, Syria and Iraq (Çakmak & Alp 2010; Coad 2015). The Mesopotamian spiny eel and its relatives in Afrotropica and south Asia present suitable elements for studies of the evolutionary history of mastacembelids across these heterogeneous areas. Day et al. (2017) included the monophyly of an African Mastacembelus radiation plus the Middle Eastern species, M. mastacembelus, based on five molecular markers and found that the position of this taxon is uncertain (as sister to the African radiation or nested within it, although support for the former scenario was stronger).

Mastacembelus is suggested to have diverged sometime during the Early Miocene, coinciding with the closure of the Tethys Sea c. 18–20 Mya (Day et al. 2017). During the Miocene and Pliocene, the major Levantine river systems and Euphrates were connected to each other (Alwan et al. 2016). Probably via headwater capture, the genus entered from the Tigris–Euphrates river system to rivers of southern Turkey and the Levant and subsequently disappeared from the Levant and North Africa. Mastacembelus is not found in eastern Iran. According to the above scenario and distribution pattern of Mastacembelus, this genus may have reached the Tigris–Euphrate basin across the Iranian land mass followed by range loss in eastern Iran, as many rivers in eastern Iran are saline and regularly dry up. However, past and present species distribution modeling is suggested.

Presently available data on the distribution of Mesopotamian spiny eel reveals that in Iran the species is restricted to Tigris, Zohreh and Persis drainages. We did not record M. armatus from the Kor River drainage in our extensive survey, in contrast to findings by Mokhayer (1981). Coad (2015) noted that the record of M. armatus by Mokhayer (1981) from the Kor River basin is probably a misidentified M. mastacembelus. There were also no observations of the Mesopotamian spiny eel from the Kor River in our fieldwork or in the literature in recent decades.

Although significant differences in morphometric traits were detected among some populations of the Mesopotamian spiny eel collected from Karakaya Reservoir, Tohma Stream in Euphrates and Tigris river in Diyarbakır by Çakmak and Alp (2010), molecular differences among populations of the species have not previously been investigated, except for a comparison by Tutar (2015) of the three populations from the Turkish part of the Tigris basin based on the combined mitochondrial DNA (16S rRNA, 12S rRNA and tRNA genes), in which no differences were found. According to our CDA results, there was a significant difference in morphometric characters of eels from different basins which could be used with great accuracy to distinguish the Tigris, Zohreh and Persis populations. According to the PCA results, although visual inspection of plotted PC1 and PC2 scores showed a high degree of overlap among the three basins, some characters of the head (HW, HM) showed significant differences among basins (sampled eels from the Zohreh basin have smaller measurements for these characters than those from the two other basins) that may be interpretable from an adaptive perspective as being caused by different ecological conditions.

The genetic distance analysis showed relatively low genetic distances among Iranian tributaries. In concordance with the genetic distance analysis, the phylogenetic trees did not show any congruent pattern including separate clusters of the Tigris, Zohreh and Persis (although all tributaries formed separate groups), revealing that the morphometric differences may be independent of COI genetic distances, highlighting a need for further molecular marker studies in the future. In any case, age estimation of the splitting of Mesopotamian spiny eel from African species, observation of different haplotypes and molecular variations of this species in the inland waters of Iran, along with the morphological differences among major river systems, may indicate a long presence of this genus in Iranian waters and may confirm the possible scenario of dispersal from Southeast Asia to Africa through Iran. More accurate conclusions will be reached when more studies are conducted on the populations of this species in the Middle East. Another reason for the morphological and molecular variations of this species in Iran may be related to the high climatic and habitat diversity from southern to western Iran. To reconstruct the distribution route of the genus Mastecembelus, population genetic studies with a sufficient number of samples and past climatic modeling are suggested.

In conclusion, from the results of the molecular divergence time of masacemblids based on COI, it can be concluded that the Mesopotamian spiny eel diverged apparently about 1 million years ago (Mya) in different river drainages of the Persian Gulf basin.

Separation among the populations of the Mesopotamian spiny eel from the different basins of Iran was shown based on morphometry; however, the genetic distances were relatively low based on COI for the main river drainages. Access to more data on this species in Mesopotamia and more comprehensive analyses would be helpful to indicate the route of Mestacembelus dispersal from Southeast Asia to Africa.

Acknowledgements

The project was financially supported by Shiraz University and was approved by the Ethics Committee of the Biology Department of the university (SU 9331009). We thank R. Sadeghi. H. Mehraban, M. Masoudi, H. Darvishnia, G. Sayyadzadeh and H. Zareian for their help in sampling. We also thank F. Zarei for help in preparing Figure 1.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the Shiraz University.