A total of 254 individuals belonging to 29 species, including both our new samples and the sequences taken from GenBank and the Barcode of Life Data systems (BOLD systems) were analysed (Suppl. material 1). A total of 44 new sequences of three molecular markers were used. The specimens were photographed and then a small piece of right pectoral fin of each individual was excised and preserved in absolute alcohol. The whole samples were stored in formaldehyde for long-term storage. All collected fish specimens were identified as described by Randall (1995). Voucher specimens are deposited in the Zoological Museum, Ferdowsi University of Mashhad (ZMFUM), Iran content-type="institution" (Suppl. material 1).

Genomic DNA was extracted from the muscle tissue following the BioGene kit protocol. The reactions volume of polymerase chain reaction (PCR) was 25 μl, including 2 μl template DNA, 9.5 μl ddH₂O, 12.5 μl Master Mix and 0.5 μl of each forward and reversed primer with 10 μM concentration. Two mitochondrial and one nuclear gene were amplified as follows: 617 bp of the CO1 was amplified using universal primers CO1LBC_F and CO1HBC_R (Ward et al. 2005), 1100 bp of the Cyt b was amplified using UniF and UniR primers (Orrell et al. 2002) and 1440 bp of the nuclear gene RAG1 was amplified using nested PCRs product (1^st PCR: 2510F and RAG1R1; 2^nd PCR: RAG1F1 and RAG1R2 primers) (López et al. 2004). PCR conditions for the COI, Cyt b, and Rag1 genes were as reported by Sanciangco et al. (2011). The PCR products were evaluated with 1.5% agarose gel for band characterisation, purified and then sent for sequencing at Microsynth Company, Switzerland. Ocyurus chrysurus (Bloch, 1791) and Rhomboplites aurorubens (Cuvier, 1829) from the family Lutjanidae Gill, 1861 were used as outgroups, since Lutjanid species has a close phylogenetic relationship with Haemulidae (Betancur-Rghton et al. 2013; Near et al. 2013).

The quality of new sequences was probed using the CLC Genomics Workbench v.3.6.5 software (QIAGEN Bioinformatics). The sequences were aligned in MAFFT v.7.463 and edited with BIOEDIT v.7.2 software (Hall 1999). Genetic interspecific distances of the mitochondrial COI gene were calculated by MEGA X (Kumar et al. 2018) using the Kimura two-parameter (K2P) model. The best-fit nucleotide substitution model was identified by the Akaike Information Criterion (AIC) in the jModelTest v.2 (Posada 2008). AIC was TIM2+I+G for COI, GTR+I+G for Cyt b, GTR+I+G for the two mtDNA regions and TIM2+I+G for the mitochondrial and nuclear markers. Maximum Likelihood (ML) analysis was performed using RAxMLv.7.2.6 (Stamatakis 2014) with 1000 replicate bootstrap. Maximum Likelihood (ML) analysis was performed using RAxMLv.7.2.6 (Stamatakis 2014) with 1000 replicate bootstrap. Bayesian Inference (BI) was run in MrBayes v.3.2.7 (Ronquist et al. 2012) with two independent Markov Chain Monte Carlo for 50,000,000 generations. The first 25% of trees were discarded as part of the burn-in phase. The remaining trees were combined to generate a 50% majority rule consensus tree. The posterior probability and bootstrap values were calculated from the BI and ML for each clade and values above 0.95 and 75 were deemed as significant supports, respectively. A haplotype network was constructed using the median-joining method in the PopART v.1.7 package (Leigh and Bryant 2015). The parsimony informative and variable sites, as well as haplotype diversity were calculated using DnaSP v.5.10 to test patterns of population structure in D. pictum and D. labiosum (Librado and Rozas 2009). Initially, genetic divergence between the putative species of Diagramma was estimated and the analysis of molecular variance (AMOVA) as Φ_ST values was performed in GenAlEx v.6.5 (Peakall and Smouse 2006) based on pairwise comparisons. Demographic changes in the effective population size were investigated for both species using two neutrality tests (Fu’s Fs and Tajima’s D) and the Ramos-Onsins and Rozas test R2 in DnaSP v.5.10 (Librado and Rozas 2009). The time elapsed since most recent population expansion was calculated herein for Diagramma pictum and D. labiosum (Rogers and Harpending 1992) using the τ = 2μt value, where τ is the time elapsed since the start of expansion in years, μ is the mutation rate per generation for the sequence (the number of base pairs per divergence rate (1% or 2% Mya). generation time (in years) and t is the time elapsed since the beginning of expansion in years. As reported by Grandcourt et al. (2011), the generation time is 2.9 years for Diagramma pictum. The mutation rate reported for the mtDNA (COI) in previous reef fish studies ranged from 1% to 2% per million years (Bowen et al. 2001).

A calibrated time-tree of Plectorhinchinae was estimated from the combined data set (COI, Cyt b and RAG1: 3157 total bp) in BEAST v.2.6.2 (Bouckaert et al. 2019). Three independent MCMC chains were run for 20 million generations each with sampling every 1,000 generations. Runs were implemented using an uncorrelated lognormal relaxed clock model to calculate rate heterogeneity and the Yule prior method for the tree prior. Convergence of runs and the effective sample size (> 200) were checked in Tracer v.1.6 (Rambaut 2016). The three runs were combined by LogCombiner v.2.6.2 (Bouckaert et al. 2019). The first 25% of trees were discarded in the TreeAnnotator v.2.6.2 (Bouckaert et al. 2019). The consensus tree with divergence time estimates was visualised using FigTree v.1.4.4 (Rambaut, 2016). Due to a lack of fossil record for Haemulidae, two calibration points were used following Near et al. (2013) and Tavera et al. (2018): (1) the Haemulidae-Lutjanidae split in early Eocene (mean 52.5, 95% HPD interval 57.2–46.8) and (2) the Parapristipoma-Plectorhinchus split in the late Eocene (mean 34.5 Ma, 95% HPD interval 42.6–27 Ma).

The ancestral area in Plectorhinchinae was inferred using the BioGeoBEARS analysis in RASP v.4 (Yu et al. 2020). Geographic distribution of each species was derived from literature (Fricke et al. 2022) and sampling of our ranges (Suppl. material 1). Four biogeographic regions were defined, based on present-day distribution and genetic diversification reported in previous studies (Ma et al. 2016; Stiller et al. 2022): Western Atlantic (WA); Eastern Atlantic (EA); Western Indian Ocean (WIO); Central Indo-Pacific (CIP). The input trees of BEAST analysis were used to infer the ancestral area probability. Six different biogeographic models were evaluated to select the best-fit using the AIC cumulative weight (AICc_wt).

Molecular species delimitation was used from four different methods, including; (A): the Automatic Barcode Gap Discovery (ABGD) (Puillandre et al. 2012) uses the pairwise genetic distance to detect a barcode gap while assuring that intra and interspecific genetic distances do not overlap and are executed in the website (http://wwwabi.snv.jussieu.fr/public/abgd) with prior intraspecific divergence (0.001 to 0.1), steps = 10 and a relative gap width of 1 using the Kimura 2-parameter model (Kimura (K80)); (B): Assemble Species by Automatic Partitioning (ASAP) (Puillandre et al. 2021), based on the distance threshold and a hierarchical clustering algorithm to create a list of best partitions. The partition score combines two metrics: the barcode gap width and the probability of panmixia (p-value), in which a lower score indicates a better partition. The method that was done using the website https://bioinfo.mnhn.fr/abi/public/asap with the Kimura (K80) model; (C): Bayesian Poisson Tree Process (bPTP) (Zhang et al. 2013), a method that estimates the number of mutations between two branching events using the branch length information of the input phylogeny. The bPTP approach was adopted on its website (http://species.hits.org/ptp/) using a non-ultrametric tree of ML inference as the input, 500,000 MCMC generations and 10% burn-in; (D): Generalised Mixed Yule Coalescent (GMYC) (Pons et al. 2006) method was run on its website (http://species.h-its.org/gmyc), using the ultrametric tree produced by BEAST v.2.6.2 (Bouckaert et al. 2019).

We assembled a total of 217 COI sequences (34 newly developed plus 183 archived; 617 bp) from 29 nominal Plectorhinchinae species. D. punctatum and Plectorhinchus pictus were sequenced for the first time. Novel COI barcodes were provided for Plectorhinchus flavomaculatus, P. gaterinus and P. gibbosus for the first time from the Northern Indian Ocean (NIO). Of 29 nominal species studied, the coalescent-based delimitations (bPTP, GMYC) yielded 27 MOTUs, whereas the distance-based methods (ABGD and ASAP) yielded 26 MOTUs (Fig. 1). ABGD and ASAP identified Plectorhinchus lessonii, and P. diagrammus and P. albovittatus and P. caeruleonothus as a MOTU. bPTP and GMYC methods considered Plectorhinchus lessonii and P. diagrammus as one MOTU. Additionally, haplotype networks, based on the COI gene, identified two main distinct haplogroups (Fig. 2A). The haplogroup I was composed of Diagramma pictum from the NIO to WP, which was separated by 15 mutation steps from haplogroups II, represented by D. punctatum from the Red Sea. In the first haplogroups, samples of D. pictum were mixed with haplotypes of D. labiosum from Australia. As suggested by the genetic distance method that uses the K2P model, the highest interspecific genetic distance was found between Plectorhinchus macrolepis and P. schotaf (21.3%). In comparison, the lowest distance was observed between D. pictum and D. labiosum (0.53%) and P. lessonii and P. diagrammus (1.2%) (Table 1). The highest intraspecific genetic variation was observed in P. gibbosus (0.69%), while values close to 0 to 0.09% were detected in most of the Diagramma and Plectorhinchus species (Table 1). Based on all delimitation analyses, D. labiosum specimens from Queensland, Australia (type locality) are genetically analogous to D. pictum and there is no evidence of their genetic differentiation. Cyt b sequence of a single D. labiosum specimen (DQ784750), deposited at the GenBank of Australia with D. pictum specimens in the phylogenetic tree, was identified as a monophyletic lineage (genetic divergence: 0.51%) (data are not shown). Seventeen haplotypes were detected based on COI sequences between D. labiosum and D. pictum (Table 2). Seventeen haplotypes were detected, based on COI sequences between D. labiosum and D. pictum (Table 2). The highest haplotype and nucleotide diversity indices were found in the samples of D. labiosum (Australia) (Hd = 0.600, π = 0.0021), while specimens of D. pictum (Hd = 0.560: NIO, π = 0.0018: WP) demonstrated the lowest values (Hd = 0.579, π = 0.001) (Table 2), where the diversity indicators were not influenced by sample size reduction. ΦST value revealed a non-significant level between Australia (D. labiosum) and IWP, except Australia (D. pictum) (Φ_ST = –1.783) (Table 2). Test of neutrality found negative and non-significant Tajima’s D values, while Fu’s Fs (-8.801) and R2 (0.038) results were significant (Table 2). These significant values (Fu’s Fs, R2), in combination with the estimates of τ for both species, indicate a shallow population history with recent expansion (Table 2). The unimodal mismatch distribution for the pooled sample of D. pictum supported a recent historical expansion (Fig. 3).

Figure 1. 

Bayesian Inference tree, based on the COI showing the results of species delimitation algorithms. The vertical bars indicate species delimitation analyses. The black circles indicate nodes with supports (ML bootstrap BP ≥ 75% and BI probability PP ≥ 0.95), orange circles (BP ≤ 75% and BI PP ≥ 0.95) and blue circles (PP ≤ 0.95 and ML BP ≤ 75%).

Figure 2. 

Median‐joining network obtained of 617 bp of mitochondrial COI for Diagramma species (A). Map depicting the geographical distribution of 17 haplotypes (B). Circle sizes are relative to haplotype frequency. Numbers between haplotypes represent mutational steps between haplotypes. Colour circles indicate species of D. pictum (red: Persian Gulf to Bangladesh, green: China and Indonesia) and D. labiosum (orange: Australia), D. punctatum (black: Red Sea).

Figure 3. 

Mismatch distribution inferred from partial mtDNA COI for the pooled dataset of D. pictum.

We analysed 254 COI DNA sequences (617 bp) which show 325 variable sites; the Cyt b gene consisting of 1117 bp showed 295 variable positions; and Rag1 consisting of 1440 bp showed 184 variable sites. Neither gaps nor stop codons were observed in the final alignment. ML and BI analyses with the combined dataset (mtDNA and nuclear) revealed that the phylogenetic relationships of sweetlips species had the highest resolution (Fig. 4). The ML and BI phylogenetic trees, based on two mtDNA (Cyt b + COI) and one nuclear (Rag1) gene, had a congruent topology, with a slight difference in their support values. The phylogenetic relationships observed in the analysis of the concatenated data (Cyt b + COI + Rag1) revealed that Parapristipoma, as a basal position, is the sister clade of the Plectorhinchus species with high bootstrap and posterior probabilities (Fig. 4). Within the Plectorhinchus clade, two well-supported subclades were formed: (1) The first clusters the most WIO species (subclade I): Plectorhinchus schotaf, P. makranensis, P. playfairi, P. sordidus, P. chubbi, P. unicolor, P. plagiodesmus, P. gibbosus, P. macrolepis, P. chrysotaenia, P. cinctus and P. pictus; (2) the second clusters the most CIP species: P. gaterinus, P. chaetodonoides, P. picus, P. albovittatus, P. lessonii, P. lineata, P. macrospilus, P. diagrammus, P. polytaenia, P. vittatus, P. orientalis, P. caeruleonothus along with D. pictum and D. punctatum indicating that the Diagramma species is nested within Plectorhinchus and should be recognised as Plectorhinchus pictum and Plectorhinchus punctatum (subclade II) (PP = 1; BS = 100%).

Figure 4. 

Phylogeny recovered by the Bayesian Inference (BI) and Maximum Likelihood (ML) analyses for the Cyt b, COI and RAG1 dataset. The black circles indicate nodes with supports (ML bootstrap BP ≥ 75% and BI probability PP ≥ 0.95), orange circles (BP ≤ 75% and BI PP ≥ 0.95) and blue circles (PP ≤ 0.95 and ML BP ≤ 75%).

The Haemulidae diverged from Lutjanidae ~ 49 Ma (early to middle Eocene, 95% HPD: 43.6–55.2 Ma) and then Plectorhinchinae divided into two main clades (A–B) and two subclades (I–II) (Fig. 5). The first split created Clade A containing Parapristipoma at ~ 35 Ma (late Eocene, 95% HPD: 32.4–37.1 Ma). Clade B diverged at ~ 21 Ma (early Miocene, 95% HPD: 18.4–23.5 Ma) into two subclades. Within each subclade, there were two periods of maximum diversification from the Miocene to Pliocene epochs. The best-fit biogeographic models recovered for the ancestral range reconstruction was the BAYAREALIKE (AICc_wt: 0.71). As suggested by the ancestral area reconstruction using the combined dataset (Cyt b + COI + Rag1), 17 dispersal and eight vicariance events occurred during the evolution of the studied Plectorhinchinae. The deepest node of early Plectorhinchinae constructed probably has its origin in the IWP with the dispersal and vicariance events occurring from the late Eocene to early Oligocene epoch (Fig. 6). However, more comprehensive taxonomic sampling of the Parapristipoma species can help find a definitive point for the origin of the subfamily Plectorhinchinae. The ancestors of the earliest Plectorhinchus lineages were supposed to be of WIO origin with vicariance and dispersal events taking place in the Miocene. Within subclade I, there seems to be a burst of species diversification in the WIO at 20 Mya. Within subclade II, the vicariance event affected speciation in the CIP as the ancestral area at 15 Mya.

Figure 5. 

Time-calibrated phylogeny, based on combined mitochondrial (COI, Cyt b) and nuclear (RAG1) genes, inferred from the BEAST analysis of the Plectorhinchinae. Numbers in the horizontal purple bars represent the 95% Highest Posterior Density intervals for each estimated event. The black circles (pp ≥ 0.95) and orange circles (pp < 0.95). A: Parapristipoma, B: Plectorhinchus.

Figure 6. 

Ancestral area reconstruction obtained using BioGeoBEARS for Plectorhinchinae. Biogeographic areas are shown in the circles and the colours in the squares indicate the current species distribution.

According to the molecular data, Plectorhinchus is divided into two subclades (Fig. 4), which display high concordance with the division proposed by Tavera et al. (2018). It can reinforce the phylogenetic value of the geographic distribution and body-colour pattern in distinguishing Plectorhinchus lineages. The first group (subclade I) mostly consists of species with uniformly dark colours (Plectorhinchus schotaf, P. makranensis, P. playfairi, P. sordidus, P. chubbi, P. unicolor, P. plagiodesmus, P. gibbosus, P. macrolepis and P. chrysotaenia) or black pots on the upper sides of the body (P. cinctus and P. pictus), along with more profound body depth (Randall 1995, 1997; McKay 2001), which are those typically distributed in the WIO (Fig. 4) (Fricke et al. 2022). The second group (subclade II) was comprised of species with spots on the whole body (P. gaterinus, P. chaetodonoides and P. picus) or horizontal stripes alongside (P. albovittatus, P. lessonii, P. lineata, P. macrospilus; P. diagrammus, P. polytaenia, P. vittatus and P. orientalis), except for P. caeruleonothus and Diagramma species (Fig. 4) (Randall 1995, 1997; Johnson et al. 2001; McKay 2001). The species includes representatives with shorter bodies from the CIP (Fricke et al. 2022). The colour variation observed in subclade II compared to subclade I can be attributed to the substrate type (Randall et al. 1997; McKay 2001; Johnson et al. 2001). The colour variation of the members of subclade II inhabiting the coral reef may be linked to social communication, mimicry and camouflage defence mechanisms to avoid predation adjusted to their habitat (Randall 1998; McMillan et al. 1999).

In our analyses, specimens from 29 nominal species were delimited between 26 and 27 consensus MOTUs. In general, low interspecific genetic divergences were observed in some of the MOTUs within Plectorhinchus, which focuses on the barcoding analysis in closely-related species. This situation was formerly reported in haemulid within Anisotremus (Bernardi et al. 2008) and Haemulon (Rocha et al. 2008). In the meanwhile, the fact that P. lessonii and P. diagrammus, along with P. albovittatus and P. caeruleonothus, grouped together in all the delimitation analyses, exhibits their very recent divergence (< 1 Ma) compared to other Plectorhinchus members. The morphological variability, including minor differences in colouration between D. pictum and D. labiosum determines taxonomic classification into subspecies or species status (Allen 1997; Johnson et al. 2001; Bogorodsky et al. 2014; Bray 2017; Fricke et al. 2022). Our molecular analyses did not yield any genetic evidence to confirm D. labiosum as a distinct species (Figs 1, 2, Table 1). A relatively low level of genetic divergence between D. pictum and D. labiosum found in both markers (COI and Cyt b) may indicate a recent speciation event (Table 1). These distances are considerably shorter than those reported in haemulid species (Johnson and Wilmer 2015; Damadi et al. 2020) and other related species of marine fishes, which is at least 2% between closely-related species (Ward et al. 2005; Hsu et al. 2007). Demographic tests (Fu’s, R2 and mismatch distribution) indicate that D. labiosum and D. pictum share a common ancestor and expansion events over the last few thousand years coincided with the Holocene glaciations that restricted larval dispersal and ocean circulation deflection (Yokoyama et al. 2001; Wilson 2013; Rosser et al. 2020). These features might have reduced coral reef habitats and colonisation from northwest Australia, especially in the reef fish (Ovenden et al. 2002; Salini et al. 2006; Sulaiman and Ovenden 2010). The incongruence between colour variation and the level of genetic divergence is not rare in fish and could be attributed to some other reasons: (1) phenotypic plasticity rather than reproductive isolation, (2) incomplete lineage sorting and (3) introgression through hybridisation (Bowen et al. 2006; Schultz et al. 2007; DiBattista et al. 2012; DiBattista et al. 2015; Bernal et al. 2017). There are several lines of evidence that favour the incomplete lineage sorting over other scenarios; however, there is still a paucity of evidence and further analyses using genomic data are warranted. Firstly, D. labiosum and D. pictum showed extensive mtDNA sharing and coalescence of haplotypes, which could be interpreted as evidence of a very recent divergence (Figs 1, 2). Secondly, their distribution is allopatric (D. labiosum: Australia and D. pictum: IWP, except Australia) (Fig. 2B) (Fricke et al. 2022). Hybridisation events and lineage sharing usually occur in sympatric and parapatric areas rather than in allopatric areas due to the greater opportunities for hybridisation in the secondary contact zones (Hobbs et al. 2009; DiBattista et al. 2016). In contrast, incomplete lineage may take place in both sympatric and allopatric areas (DiBattista et al. 2012; DiBattista et al. 2015). Thirdly, regional colour morphs usually represent minor molecular divergences due to recent genetic divergence, which may not be related to ecological and behavioural adaptations (Taylor and Hellberg 2003; Drew et al. 2008; Drew et al. 2010). Fourthly, the Diagramma species are distinguished only based on colour patterns (Johnson et al. 2001; Heemstra and Heemstra 2004; Bogorodsky et al. 2014) which evolve more quickly than the other morphological traits and genetic characters (Bowen et al. 2006; Schultz et al. 2007). Colouration could be an important criterion to distinguish closely-related species in haemulid fishes (Johnson and Wilmer 2015; Damadi et al. 2021). Fifthly, the most recent divergence occurred between D. pictum and D. labiosum around 28–56 million years ago (Table 1), which is more recent than the complete lineage sorting time. A divergence time of 2.9 Mya seems sufficient for these lineages for complete sorting (1× generation time) (Grandcourt et al. 2011).

Our biogeographic findings exhibit a complex history highlighting the role of both biological and geographical factors in the diversification and evolution within the Haemulids. The first diversification, i.e. the split of Haemulids from Lutjanids as a sister lineage in the Eocene, could be related to the ice sheet formed in the Antarctic and reduced atmospheric CO₂ level (Zachos et al. 2001; Coxall et al. 2005; Lear et al. 2008). These climatic and oceanographic changes resulted in a sea-level drop, upwelling and extinction of warmer climate fauna (Zachos et al. 2001; Livermore et al. 2007; Ma et al. 2016). This extinction formed vacant niches and the emerging Indo-Australian Archipelago probably gave rise to new reef habitats and high species-richness in the cooler oceans (Cowman and Bellwood 2011; Ma et al. 2016). This divergence time corresponds to the one shown by Near et al. (2012) for the same historical event in Haemulids. The historical biogeography suggests that the common ancestor of Plectorhinchinae might have its origin in the IWP and subsequently spread and diversified in the western Indian Ocean and the Central Indo-Pacific. This scenario could be firmly demonstrated by complete biogeographic sampling of some taxa and extra genetic/genomic analyses. For Plectorhinchinae, divergence was first observed in East Atlantic (EA) species (Parapristipoma) with the fewest species from the Plectorhinchus lineage, approximately 33.8 Mya, which coincided with the onset of the terminal Tethyan event (TTE) that disconnected West Tethyan and Arabian marine biodiversity hotspots (Adams et al.1983; Harzhauser et al. 2007). The diminishing connectivity between two regions might have set the scene for vicariant splits, which appear to be linked to the early diversification of this clade. This divergence time is also congruent with the EA and IWP divergence in other species of reef fish, particularly the Groupers and Boxfishes (Santini et al. 2013; Ma et al. 2016). Within the Plectorhinchus IWP Clade, two factors, the Indo-Australian Archipelago (IAA) uplift and major global cooling, were associated with the diversification of sweetlips subclades. These factors were inferred from other reef fish in the Miocene period (Zachos et al. 2001; Bellwood and Wainwright 2002; Shevenell et al. 2004; Cowman and Bellwood 2011; Ma et al. 2016). An estimate of divergence times of each Plectorhinchus subclade suggests further diversification occurred during the Miocene that slow indicated by the long branches, followed by slower diversification in Pliocene to Pleistocene epochs (Fig. 5). The IAA uplift during the early to mid-Miocene in subclade II (CIP Plectorhinchus) created new islands, coinciding with the lowering of sea level, leading to novel shallow reef habitat discontinuity. During Plio-Pleistocene, Diagramma species diverged from each other around 4.1 to 2.3 Mya (Fig. 5) and D. punctatum, an endemic species, formed in the Red Sea, coinciding probably with high temperature and salinity, sea-level changes during glacial cycles and the narrow Strait of Bab al Mandab. These factors led to a decrease in unsuitable habitat and the water flow through this strait (Rihm and Henke 1998; Riegl and Purkis 2012; Stevens et al. 2014; Torquato et al. 2019) that has been documented in marine fish (Froukh and Kochzius 2007; DiBattista et al. 2013; Bogorodsky et al. 2017). The major global cooling could lay the ground for emergence of new environmental niches likely driven by adaptive radiations in subclade I (WIO) during the Miocene (Zachos et al. 2001; Shevenell et al. 2004).