Introduction
The Persian Gulf and Gulf of Oman, in the north-west of the Indian Ocean, are connected by the Strait of Hormuz. In the Persian Gulf, the thermohaline circulation system and other physical factors (e.g. winter surface isotherms) play a crucial role in shaping the current species distribution (Sheppard et al., Reference Sheppard, Al-Husiani, Al-Jamali, Al-Yamani, Baldwin, Bishop, Benzoni, Dutrieux, Nicholas, Durvasula, Jones, Loughland, Medio, Nithyanandan, Pillingm, Polikarpov, Price, Purkis, Bd, Saburova, Samimi-Namin, Taylor, Wilson and Zainal2010). The surface water of the Indian Ocean with lower-salinity enters the Persian Gulf through the northern part of the Hormuz Strait and is deflected along the Iranian coast, forming a basin-wide cyclonic circulation in the southern Persian Gulf (Yao, Reference Yao2008). Relatively, the Persian Gulf is geologically young, as its basin was completely dried during the Pleistocene glacial periods about 17,000 years ago (Barth & Khan, Reference Barth, Khan, Abuzinada, Barth, Krupp, Böer and Al Abdessalaam2008). Around 12,500 years ago, the Strait of Hormuz opened up, and marine water started to enter the central basin. The present shoreline was established about 6000 years ago (Lambeck, Reference Lambeck1996); therefore, the corresponding biota is also young. Consequently, the Persian Gulf has an impoverished fauna and flora, not only because of its young age, but mostly because of its harsh environment (Sheppard et al., Reference Sheppard, Andrew and Roberts1992).
Studying population genetic structures enhances our knowledge of how populations evolve under evolutionary forces such as gene flow. Small passively dispersed marine organisms can maintain gene flow between populations over long distances (Gómez et al., Reference Gómez, Carvalho and Lunt2000). On the other hand, many taxa show high levels of genetic differentiation and form cryptic species complexes, even across small spatial scales within the range of their dispersal (e.g. Zeller et al., Reference Zeller, Reusch and Lampert2006; Fontaneto et al., Reference Fontaneto, Boschetti and Ricci2008; Xiang et al., Reference Xiang, Xi, Wen, Zhang, Wang and Hu2011).
Genetic diversity of holoplanktonic organisms such as copepods has been largely underestimated in the past owing to their high dispersal capacities, large population sizes, wide geographic ranges and high tolerance to variable environmental factors (Goetze, Reference Goetze2003). Pelagic species may become geographically segregated (allopatric speciation) by oceanographic barriers, isolation by distance, selection due to environmental adaptations or recent historical events (Palumbi, Reference Palumbi1994; Blanco-Bercial et al., Reference Blanco-Bercial, Álvarez-Marqués and Bucklin2011). Historical events, such as glaciation events or bottleneck effects followed by a sudden population expansion may result in establishment of unique haplotypes (Nei et al., Reference Nei, Maruyama and Chakraborty1975; Nuwer et al., Reference Nuwer, Frost and Armbrust2008). Moreover, ecology and habitat affinity of populations are probably the main drivers for the speciation of the diverse marine copepod fauna apart from environmental barriers such as continental landmasses and large-scale ocean circulation (Goetze, Reference Goetze2005). Assortative mating can lead to sympatric and micro-allopatric speciation (Palumbi, Reference Palumbi1994; Doebeli et al., Reference Doebeli, Dieckmann and Macdonald2000) and genetic variability might occur, even though morphological differences are absent, leading to the formation of cryptic and pseudo-cryptic species complexes (e.g. Goetze, Reference Goetze2003; Goetze & Ohman, Reference Goetze and Ohman2010; Cornils et al., Reference Cornils, Wend-Heckmann and Held2017). Furthermore, evolutionary processes are not only driven by selective pressures, but also influenced by genetic drift, random changes of gene frequencies in populations (Lande, Reference Lande1976).
Symbiotic copepods have one or more planktonic phase for dispersal, infection, host switching, mating, and presumably, predator-avoidance (Kearn, Reference Kearn2004; Huys, Reference Huys, Martin, Olesen and Høeg2014). Usually, the primarily planktonic naupliar stages play a key role in dispersal. The first copepodid stage is infective, and the subsequent stages, including adults, then typically establish a symbiotic association with a host organism.
Copepods of the genus Clausidium Kossmann, Reference Kossmann1874 are external associates of burrowing decapod families Callianassidae and Upogebiidae (Kihara & Rocha, Reference Kihara and Rocha2013). Males of these copepods are often found attached to their larger female partner. The dispersal capacity of adult Clausidium individuals is not known, due to the cryptic lifestyle of their hosts. Members of this copepod genus are certainly not proficient swimmers, as they resemble a larval stage compared with free-living marine planktonic copepods. Little is known about the lifestyle, the dispersal capacities and the genetic variability of this parasitic copepod genus. The life cycle is likely similar to species belonging to Hemicyclops as this is a common member of the same family (Clausidiidae). Descriptions of some parts of life cycles of Hemicyclops have been published but these have mostly focused on the copepodid stages (Itoh & Nishida, Reference Itoh and Nishida2007), but it is doubtful if naupliar stages have been described. Based on Hemicyclops species' life cycle, it would be expected that members of genus Clausidium should have six naupliar stages followed by six copepodid stages (the sixth counts as an adult) (Geoff Boxshall, personal communication, June 2018). This is the typical pattern for many families of parasitic or associated copepods. The nauplius stages are generally planktonic and the primary dispersal phase. The first copepodid stage is the infective stage that locates the host, and the remaining copepodid stages are found associated with the host. If the copepod is only loosely associated with the host, copepodid stages of the Clausidiidae may occasionally occur in plankton samples (these were erroneously given the generic name Saphirella; Geoff Boxshall, personal communication, June 2018).
Two species of Clausidium occur in the study region, Clausidium persiaensis in the Persian Gulf and Clausidium iranensis in the Gulf of Oman (Sepahvand et al., Reference Sepahvand, Rastegar-Pouyani, Kihara and Momtazi2017, Reference Sepahvand, Kihara and Boxshall2019). In the present study, we used the mitochondrial marker gene, cytochrome oxidase c subunit 1 (COI) to determine the population structure of these two species in order to establish the importance of environmental factors and life history traits on intraspecific diversity of either species. We hypothesize that the populations of both species undergo similar evolutionary forces and that this is reflected in a similar genetic structure.
Materials and methods
Study areas, sampling and sample processing
In total, 130 specimens were collected from bodies of the burrowing shrimps Callianidea typa (C. persiaensis) and Neocallichirus jousseaumei (Clausidium iranensis) from five localities along the coastal region of the Persian Gulf and the Gulf of Oman (Figures 1 and 2, Table 1). Single individuals, from each locality, were preserved in 99.5% ethanol. All specimens were morphologically identified to species level using the original species descriptions and a identification key to Clausidium copepod species (Sepahvand et al., Reference Sepahvand, Kihara and Boxshall2019).
Map of sampling stations in the Persian Gulf and the Gulf of Oman. The arrows show the circulation of water within the Persian Gulf.
(A) Callianidea typea (female), habitus, dorsal view, and (B) its parasite Clausidium persiaensis. Confocal laser scanning microscopy maximum projections, female habitus; (C) Neocallichirus jousseaumei. Habitus, lateral view; (D) Clausidium iranensis. Confocal laser scanning microscopy maximum projections, Female habitus. Scale bar, A and C 1 cm, B and D 100 μm.
Detail of collection localities for two species of Clausidium
DNA isolation, PCR and sequencing
From 79 morphologically identified specimens, DNA was extracted using the whole specimen in 40 μl of Chelex (InstaGene Matrix, Bio-Rad) or 100 μl DNeasy Tissue Extraction Kit (QIAGEN, Hilden, Germany) following Gollner et al. (Reference Gollner, Stuckas, Kihara, Laurent, Kodami and Martinez Arbizu2016).
A 658 base pair fragment of the mitochondrial gene COI, was amplified using universal primers LCOI (5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′) and HCOI (5′-TAA ACT TCA GGG TGA CCA AAA AAT CA-3′) (Folmer et al., Reference Folmer, Black, Hoeh, Lutz and Vrijenhoek1994). PCR conditions for the above-mentioned primers were described by Gollner et al. (Reference Gollner, Stuckas, Kihara, Laurent, Kodami and Martinez Arbizu2016). Sequences were manually edited using Bioedit v.7.1.30 (Hall, Reference Hall1999) and Chromas v2.23 (available at www.technelysium.com.au). To check for stop codons, COI sequences were translated into amino acids using Mega 6 (Tamura et al., Reference Tamura, Stecher, Peterson, Filipski and Kumar2013).
Population genetic and demographic analysis
To depict the evolutionary relationships among the mtDNA haplotypes, a haplotype network was constructed using TCS version 1.18 (Clement et al., Reference Clement, Posada and Crandall2000), and then rearranged and drawn using the tool tcsBU (Múrias dos Santos et al., Reference Múrias dos Santos, Cabezas, Tavares, Xavier and Branco2015). The number of mtDNA haplotypes, the levels of haplotype diversity (h) and nucleotide diversity (π) were calculated for each population, using DnaSP 5.00 (Librado & Rozas, Reference Librado and Rozas2009). The pairwise fixation index (FST) was evaluated in Arlequin 3.5 (Excoffier & Lischer, Reference Excoffier and Lischer2010), with 10,000 permutations to measure the genetic differentiation among the populations.
To infer the historical demography, separately for each species, Tajima's D statistics of neutrality (Tajima, Reference Tajima1989) and the frequency distribution of pairwise differences between mtDNA haplotypes (i.e. mismatch distribution) were employed. A significantly negative deviation from zero can be interpreted as being the result of past population expansion and/or purifying selection, whereas a significantly positive value can result from balancing selection and/or a decrease in population size. Tajima's D values was evaluated using the coalescent algorithm implemented in DnaSP 5.0, in which the observed value is compared to a null distribution generated by 10,000 replicates, given an empirical population sample size and the observed number of segregating sites. To determine whether population has gone through an expansion, the mismatch distribution, with Harpending's raggedness index (Harpending, Reference Harpending1994), was calculated using Arlequin 3.5 (Excoffier & Lischer, Reference Excoffier and Lischer2010) under the sudden expansion model with 1000 bootstrap replicates. If the raggedness index value is insignificant (P > 0.05), demographic expansion cannot be rejected (Schneider & Excoffier, Reference Schneider and Excoffier1999).
Results
A COI fragment of 658 bp length of 74 specimens was sequenced, including 49 and 25 sequences belong to C. persiaensis and C. iranensis respectively. For C. iranensis, 10 segregating sites, substantial haplotype diversity (0.76), low nucleotide diversity (0.0016) and only one parsimony informative site were identified (Table 2). The haplotype network consists of 10 haplotypes (accession number pending), of which eight are rare haplotypes with a single individual and two major haplotypes sharing the remaining 17 individuals. All orphan haplotypes, two from Djod and six from Chabahar populations, were separated only by one mutation step from the main haplotype, giving a star-shaped topology (Figure 3). There is a significant genetic differentiation between the two populations belonging to the Gulf of Oman, Djod and Chabahar (Φst = 0.51, P < 0.05). The population's demography of C. iranensis supports a recent expansion model due to significant negative value of the neutrality test (Tajima's D = −1.93; P < 0.05). The skewed unimodal shape of mismatch distribution and insignificant value of the raggedness index (r = 0.15, P < 0.01), in line with results from neutrality test and network topology, reflects a scenario of recent population expansion (Figure 4, Table 3).
COI haplotype network of Clausidium iranensis. Coloured circles represent individual haplotypes and the sizes of the circles are proportional to the frequency of the haplotypes. Each point mutation is symbolized by a mark on the lines connecting haplotypes.
Mismatch distributions for the COI locus for C. iranensis.
Genetic diversity indices for C. iranensis
N, number of individuals sampled; S, number of segregating sites; H n, number of haplotypes; H d, haplotype diversity; π, nucleotide diversity.
Mismatch distribution and significance values (P) of Harpending's raggedness index (HRI) and sums of squared deviation (SSD) for C. iranensis and C. persiaensis
For C. persiaensis the data set included 24 parsimony informative sites, 44 segregating sites, substantial haplotype diversity (0.96) and nucleotide diversity (0.0078) (Table 4). The haplotype network includes 29 haplotypes, with one main haplotype including 13 individuals sharing by all populations. Twenty-five haplotypes represent private haplotypes (Figure 5), of which four haplotypes include more than one specimen. The haplotype network indicates close linkage among haplotypes as the Hormozgan population was shared among the samples collected from Qeshm and Bushehr and its private haplotypes are partially close to other two populations (Figure 5).
Median-joining haplotype network for COI sequences from C. persiaensis taxa. Coloured circles represent individual haplotypes and the sizes of the circles are proportional to the frequency of the haplotypes. Each point mutation is symbolized by a mark on the lines connecting haplotypes.
Genetic diversity index of populations of C. persiaensis
N, number of individuals sampled; S, number of segregating sites; H n, number of haplotypes; H d, haplotype diversity; π, nucleotide diversity.
F-statistics showed insignificant genetic differentiation between Hormozgan and Busheher (Φst = −0.09, P < 0.05) and between Bushehr and Qeshm (Φst = −0.00032, P < 0.05), otherwise it was significant between Hormozgan and Qeshm (Φst = 0.067, P < 0.05). AMOVA test, shows no detection of any genetic structure among populations; F ST = 0.07. Tajima's D values were −1.67 (P < 0.05), indicating a recent event of population expansion (Figure 6, Table 3). The skewed bimodal shape of mismatch distribution and significant value of the raggedness index (P = 0.009), in contrast to results from the neutrality test, cannot demonstrate a scenario of recent population expansion (Figure 6, Table 3). Therefore, it can be excluded that populations are constant.
Mismatch distributions for the COI locus for C. persiaensis.
Discussion
This study is the first to assess the population structure of copepods associated with ghost shrimp between the Persian Gulf and Gulf of Oman. The two Clausidium species were studied alive in close association with their ghost shrimps (Sepahvand et al., Reference Sepahvand, Rastegar-Pouyani, Kihara and Momtazi2017, Reference Sepahvand, Kihara and Boxshall2019). To our knowledge, there is no information on the genetics, behaviour and ecology of Clausidium copepods, because of the cryptic lifestyle of their hosts. Our analyses were based on a comparatively large number of specimens (total N = 130) of two different Clausidium species, from locations spanning across more than 2000 km of coastline. Contrary to previous findings of significant population genetic structure across the region in different taxa (e.g. Afkhami et al., Reference Afkhami, Schubart and Naderloo2016; Ghanbarifardi et al., Reference Ghanbarifardi, Aliabadian and Esmaeili2018), the smaller values of Φst observed in the present study may derive from greater gene flow between subpopulations of Clausidium copepods, in particular in C. persiaensis. The real distribution and connectivity of intertidal animals in the region and closely related areas (i.e. West Indian Ocean) are determined by oceanographic regime, environmental conditions and historical events (Tsang et al., Reference Tsang, Achituv, Chu and Chan2012; Afkhami et al., Reference Afkhami, Schubart and Naderloo2016; Rahimi et al., Reference Rahimi, Rezvani Gilkolaie, Ghavam Mostafavi, Jamili and Rahnema2016; Ghanbarifardi et al., Reference Ghanbarifardi, Aliabadian and Esmaeili2018).
Our study suggests that each species of Clausidium copepod shows a pattern of generally high haplotype diversity and high connectivity, regardless of their host species. Clausidium persiaensis shows no spatial genetic differentiation, while C. iranensis, even though populations are in proximity, slight spatial genetic differentiation was observed. A plausible scenario for such a haplotype diversity of Clausidium copepods is likely to be the difference in occupancy levels observed between the two copepod species as a consequence of the behaviour or physiology of their hosts, since their attachment mechanisms are the same (Sepahvand et al., Reference Sepahvand, Brown and Gholamifard2020). The hosts are different in burrowing patterns, grooming and nutrition behaviours (Griffis & Suchanek, Reference Griffis and Suchanek1991; Sepahvand et al., Reference Sepahvand, Sari, Tudge and Bolouki2014) and these parameters may influence the number of copepods occupying each host and consequently the haplotype diversity. Sepahvand et al. (Reference Sepahvand, Brown and Gholamifard2020), claimed that, although the number of Clausidium copepods increases with the host size, the two host species vary in the degree of symbiont invasion, with large C. typa (host for C. persiaensis) hosting ~7 times as many symbionts as the similarly sized N. jousseaumei (host for C. iranensis). Factors such as host gender and host species also affect the density of copepods (Corsetti & Strasser, Reference Corsetti and Strasser2003).
Dispersal capacity and oceanographic barriers may play a central role for genetic differentiation in copepods (Goetze & Ohman, Reference Goetze and Ohman2010; Norton & Goetze, Reference Norton and Goetze2013). Hence, population genetic structure of Clausidium copepods could be driven by the oceanography and biological traits along the southern coast of Iran, since the pattern of genetic connectivity corresponds to the ocean current in the Persian Gulf and Gulf of Oman. It is likely that there are different magnitudes of gene flow taking into account the dissimilar dispersal abilities of the two species. Whereas both adult burrowing shrimps and Clausidium copepods are able to move, their movement is limited compared with free-living species.
In summary, it seems that ecological, morphological, life-history differences of the host and Clausidium copepods as well as abiotic (historical events, ocean current) parameters construct the population structure of the copepods. However, based on Sepahvand et al. (Reference Sepahvand, Brown and Gholamifard2020), host specificity and microhabitat selection of Clausidium copepods as the most important parameters of population size may have the most effective role on haplotype diversity of this group of copepods. Future studies should focus on getting a better understanding of the genetic population structure of C. iranensis and C. persiaensis as well as connectivity among populations and detailed analysis of life cycle strategies of their hosts.
Acknowledgements
We are very grateful to Dr Sahar Khodami from Senckenberg am Meer, German Center for Marine Biodiversity Research, Wilhelmshaven for help and encouragement during this study. We appreciate Dr Astrid Corlins from Alfred-Wegener-Institut, Helmholtz-Zentrum für Polar- und Meeresforschung for her helpful comments on the manuscript. The authors express their gratitude to Dr Alireza Sari from University of Tehran for his help in this study.
