A molecular assessment of species boundaries and relationships in the Australian brine shrimp Parartemia (Anostraca: Parartemiidae)

Specimen collection and preservation

We collected Parartemia specimens from 84 salt lakes in Australia (mainly Western Australia) between September 2017 and August 2022 (Fig. 1 and Supplementary Table S1). Specimens were collected using a dip net and then euthanised by freezing and preserved in 100% ethanol. Specimens from a site in Tasmania (provided by the Tasmanian Museum and Art Gallery) and eight sites in Western Australia (from the Department of Biodiversity, Conservation and Attractions, the Western Australian Museum and the Stantec Australia Pty Ltd) were also used (details in Supplementary Table S1).

Fig. 1.

Approximate locations of Parartemia collection sites. Sites in close proximity may not always be distinguishable in the figure (see detailed information in Supplementary Table S1). The hydrological information, shown in blue, is sourced from the national surface water database of Geoscience Australia (www.ga.gov.au). The full Australian map includes the state and territory boundaries, WA, Western Australia; SA, South Australia; NT, Northern Territory; Qld, Queensland; NSW, New South Wales; ACT, Australian Capital Territory; Vic., Victoria; and Tas., Tasmania. No specimens of P. minuta were collected in this study but GenBank sequences for this species from Lake Buchanan in Qld and an unspecified site in the Paroo area in NSW were included (details in Supplementary Table S1).

Identification of morphospecies

All Parartemia specimens were identified to morphospecies using the morphological criteria of Timms (2012). All known morphospecies were represented in the samples except for P. auriciforma, P. yarleensis and P. minuta. Parartemia auriciforma is only known from one site in central Australia, P. yarleensis is only known from inland South Australia and P. minuta occurs in inland areas in South Australia, Victoria, New South Wales and Queensland (Timms et al. 2009). Our attempts to obtain either fresh or ethanol-preserved specimens of these three species were unsuccessful. However, for P. minuta, we included 28S and 16S sequences from GenBank (see below) and morphological data from Timms (2012). We collected two groups of specimens whose morphology did not match those of any described species (see ‘Results’) – herein these groups are referred to as Parartemia sp. ‘y’ and Parartemia sp. ‘z’.

DNA extractions, PCR amplification and sequencing

Genomic DNA was extracted from thoracic segments V through VII or thoracopods III through VII using a Masterpure Complete DNA and RNA Purification Kit (Epicentre) following the manufacturer’s instructions. Taxonomically significant characters (head, genital segments and abdomen) were left intact for later cross-checking with molecular findings if required. Negative controls, i.e. assays with reagents but no added tissue or DNA, were included in every extraction to check for contamination.

One nuclear genetic region (28S) and two mitochondrial genetic regions (16S and COI) were used. Details of the PCR primers and amplicon lengths are given in Table 1. Except for the universal primers LCO1490 and HCO2198 (Folmer et al. 1994), all primers were designed specifically for Parartemia because it was not possible to amplify the target region from all species or all populations of a species with preexisting primers (see Table 1, Supplementary Tables S2 and S3).

PCR reaction volumes were 25 μL containing 5 µL of GoTaq Reaction Buffer (Promega), 0.5 μL of dNTPs (10 mM per nucleotide), 0.25 μL of each of forward and reverse primers (10 µM), 2 μL of MgCl₂ (25 mM), 0.35 μL of bovine serum albumin (10 µg µL^–1), 0.125 μL of GoTaq G2 Hot Start Taq Polymerase, 1 μL of DNA and adjusted to the final volume using PCR grade water. PCR reactions for all three markers were (i) 95°C initial denaturation temperature for 5 min; (ii) 40 cycles of denaturation at 94°C for 1 min, annealing at 48°C for 1 min and extension at 72°C for 45 s; and (iii) a final extension at 72°C for 7 min. PCR products were purified using Exo-SAP purification (Dugan et al. 2002) and sequenced in the forward and reverse directions in an automatic ABI 3700 sequencer (Applied Biosystems) by Macrogen Inc. (South Korea).

Sequence data

The sequencing chromatograms were visualised in Chromas (ver. 2.6.5, Technelysium Pty Ltd, Australia). Forward and reverse sequences were compared, and any ambiguities were corrected. The consensus 16S sequences were aligned using MAFFT online version (see http://mafft.cbrc.jp/alignment/server/) with the G-INS-i strategy, whereas the consensus COI sequences were aligned using MUSCLE (ver. 5, see http://www.drive5.com/muscle/; Edgar 2004) in MEGA X (ver. 10.2.6, see https://www.megasoftware.net; Kumar et al. 2018). The 28S sequences were aligned as described for 16S except that we used the Q-INS-i strategy and made manual adjustments by eye in regions containing indels. We found no evidence that nuclear copies of COI had inadvertently been included in our Parartemia dataset as, for example, there were few amino acid substitutions, and no indels or stop codons, in the translated sequences (see Raupach and Radulovici 2015). Haplotypes in the 28S, 16S and COI datasets were identified using DnaSP (ver. 6.12.03, see http://www.ub.edu/dnasp/; Rozas et al. 2017). All new haplotypes have been deposited in GenBank with accession numbers listed in Supplementary Table S1. This study also used GenBank sequences for Parartemia (including one 28S and two 16S haplotypes of P. minuta) and a range of outgroup taxa (see Supplementary Tables S1 and S4). A multigene (28S, 16S and COI) concatenated dataset was assembled in MEGA X (Kumar et al. 2018), using individuals for which data from at least two of the three targeted genetic regions were available.

Phylogenetic analysis

Phylogenetic analyses were conducted using Bayesian inference (BI) and maximum likelihood (ML) frameworks on the concatenated dataset and on single locus 28S, 16S and COI datasets. The former is widely recognised as the most reliable approach for phylogenetic analysis (Wiens and Moen 2008). The concatenated dataset was used to investigate species relationships within Parartemia. The 16S and COI markers were also used to investigate the validity of Parartemia morphospecies, as they are effective in identifying crustacean species (Remigio et al. 2001; Costa et al. 2007; Raupach and Radulovici 2015).

The best nucleotide substitution models (GTR + I + G for 28S and TrN + G + I for 16S and COI) were selected using jModelTest 2 (ver. 2.1.9, see https://github.com/ddarriba/jmodeltest2; Darriba et al. 2012) based on the Bayesian information criterion (BIC). Maximum likelihood molecular clock tests conducted in MEGA X (Kumar et al. 2018) revealed that the partitioned concatenated dataset and the single locus datasets did not adhere to the assumptions of a strict clock. BI analysis for each dataset was separately conducted in BEAST (ver. 1.10.4, see https://beast.community/beast; Suchard et al. 2018), employing the identified substitution models (see above) with the uncorrelated relaxed clock and coalescent constant population size as tree priors. The analysis was run for 50 million generations, and the estimated sample size (ESS) was checked by looking at the log-output file in Tracer (ver. 1.7.2, see https://github.com/beast-dev/tracer/releases/tag/v1.7.2; Rambaut et al. 2018). A burn in of 25% of the initial trees was discarded, and the final tree was produced in TreeAnnotator (a BEAST-distributed program) and visualised in FigTree (ver. 1.4.4, A. Rambaut, see http://tree.bio.ed.ac.uk/software/figtree/). The analysis was also performed for each dataset using the uncorrelated relaxed clock and Yule process as tree priors, but the ESS was low (<200), even after increasing the MCMC chain length to 100 million generations, and thus the results were discarded.

The ML phylogenetic analysis was performed separately for the concatenated, 28S, 16S and COI datasets on the IQ-TREE web server (see http://iqtree.cibiv.univie.ac.at; Trifinopoulos et al. 2016) using 5000 ultrafast bootstrap replicates and the above-mentioned substitution models.

Validity of Parartemia morphospecies and pairwise distances

Three species delimitation analyses were employed to assess the validity of Parartemia morphospecies using the 16S and COI datasets (without outgroups). Firstly, the Assemble Species by Automatic Partitioning (ASAP), which relies on pairwise genetic distances (Puillandre et al. 2021), was conducted on the ASAP website (see https://bioinfo.mnhn.fr/abi/public/asap) using the Kimura (K80) substitution model (Kimura 1980) and default settings. Of the 10 best ASAP partition schemes included in the results, the one with the lowest ASAP score and another one with the smallest number of partitions were chosen. Secondly, the maximum likelihood implementation of the multi-rate Poisson tree processes (mPTP) (Kapli et al. 2017) was performed on the mPTP webserver (see https://mptp.h-its.org) with default settings using a BI phylogenetic tree generated by BEAST (Suchard et al. 2018) based on the same parameters as outlined in the ‘Phylogenetic analysis’ section (see above). Lastly, the General Mixed Yule Coalescent (GMYC) analysis (Fujisawa and Barraclough 2013) was performed in R (ver. 4.3.3, R Foundation for Statistical Computing, Vienna, Austria, see https://www.r-project.org/) using the same BI phylogenetic tree as used for mPTP (see Michonneau 2016 for further details).

MEGA X (Kumar et al. 2018) was used to compute uncorrected p-distances and Kimura two-parameter (K2P) (Kimura 1980) distances between Parartemia haplotypes in the 28S, 16S and COI datasets.

General information

After aligning and trimming, a total of 28 28S (674 bp excluding gaps), 100 16S (474 bp excluding gaps) and 161 COI (658 bp) haplotypes of Parartemia were identified in sequences from 50, 113 and 232 individuals respectively.

Substantial genetic divergence was present in the 16S (e.g. p-distance of up to 20.34%) and especially the COI (e.g. p-distance of up to 25.23%) genetic regions in Parartemia. As expected, divergence in the 28S region was more limited (e.g. p-distance of up to 2.8%), with few differences between most species (Supplementary Table S5). Accordingly, the relationships among Parartemia species were poorly resolved in 28S phylogeny (Supplementary Fig. S1), although the 28S data did provide strong support for the monophyly of Parartemia (see below).

Monophyly of Parartemia and species relationships

All Parartemia haplotypes formed a well-supported monophyletic group in both BI and ML concatenated phylogenetic trees (BPP = 1 and 97% bootstrap value; Fig. 2 and Supplementary Fig. S2), as well as in the BI and ML 28S (Supplementary Fig. S1) and 16S trees (see Fig. 3 and Supplementary Fig. S3) and the ML COI tree (see Supplementary Fig. S4). In the BI COI tree, all Parartemia haplotypes formed a monophyletic group, but the node support was low (see Fig. 4).

Fig. 2.

Bayesian inference (BI) phylogenetic tree for the Parartemia concatenated dataset (COI, 16S and 28S). Maximum likelihood (ML) phylogenetic tree is available in Supplementary Fig. S2. Bayesian Posterior Probability (BPP, when ≥0.80) and bootstrap values from the ML tree (when ≥80%) are indicated at nodes (BPP/bootstrap). For nodes where one value was above the threshold and the other was below, the latter is indicated by hyphen (-). BPP values of 1 and bootstraps of 100% are indicated by asterisks (*). Node with bootstrap value (x) indicates that the node-level species composition is not supported by the ML phylogenetic tree. Species names are given based on Parartemia morphotaxonomy (morphospecies). Coloured bars to the right of the tree indicate species groups identified within Parartemia in this study, named alphabetically (A–E).

Fig. 3.

Bayesian inference (BI) phylogenetic tree of Parartemia 16S haplotypes. Maximum likelihood (ML) phylogenetic tree is available in Supplementary Fig. S3. Bayesian Posterior Probability (BPP, when ≥0.80) and bootstrap values from the ML tree (when ≥80%) are indicated at nodes (BPP/bootstrap). For nodes where one value was above the threshold and the other was below, the latter is indicated by hyphen (-). BPP values of 1 and bootstraps of 100% are indicated by asterisks (*). Each haplotype is denoted by its GenBank accession number followed by its morphospecies name. The green, blue, magenta and orange vertical rectangles indicate recovered species partitions based on Assemble Species by Automatic Partitioning (ASAP) lowest score scheme (ASAP1), ASAP least partition scheme (ASAP2), multi-rate Poisson Tree Processes (mPTP) and Generalised Mixed Yule Coalescent (GMYC) analyses respectively. Species hypotheses (SH) determined by this study are indicated by black bars under the SH heading. Grey-coloured bars to the right of the coloured rectangles indicate species groups (SG) identified within Parartemia in this study, named alphabetically (A–E).

Fig. 4.

Bayesian inference (BI) phylogenetic tree of Parartemia COI haplotypes. Maximum likelihood (ML) phylogenetic tree is available in Supplementary Fig. S4. Bayesian Posterior Probability (BPP, when ≥0.80) and bootstrap values from the ML tree (when ≥80%) are indicated at nodes (BPP/bootstrap). For nodes where one value was above the threshold and the other was below, the latter is indicated by a hyphen (-). BPP values of 1 and bootstraps of 100% are indicated by asterisks (*). Each haplotype is denoted by its GenBank accession number followed by its morphospecies name. The green, blue, magenta and orange vertical rectangles indicate recovered species partitions based on Assemble Species by Automatic Partitioning (ASAP) lowest score scheme (ASAP1), ASAP least partition scheme (ASAP2), multi-rate Poisson Tree Processes (mPTP) and Generalised Mixed Yule Coalescent (GMYC) analyses respectively. Species hypotheses (SH) determined by this study are indicated by black bars under the SH heading. Grey-coloured bars to the right of the coloured rectangles indicate species groups (SG) identified within Parartemia in this study, named alphabetically (A–E; species group E is not shown due to the absence of COI data for P. minuta).

Although some aspects of the relationships among Parartemia species were not well resolved, particularly for deeper divergences, some species groups and divergent species were evident. The haplotypes of five species (P. bicorna, P. contracta, P. informis, P. laticaudata and P. veronicae) were invariably grouped together in a single usually well-supported clade in the BI and ML concatenated, 16S and COI trees (species group A in Fig. 2–4). A sixth species, P. triquetra, for which there were no 28S or COI data, was also included in this group in the 16S trees (Fig. 3 and Supplementary Fig. S3). Similarly, Parartemia sp. ‘z’ and the P. boomeranga and P. longicaudata morphotypes always formed a single group (species group B in Fig. 2–4). Parartemia cylindrifera and Parartemia sp. ‘y’ also consistently grouped together, forming a well-supported clade in the BI trees; this clade was present but less strongly supported in the ML trees (species group C in Fig. 2–4). Parartemia zietziana and the different lines of P. purpurea (see below) also usually grouped together but the clade was not always well supported (species group D in Fig. 2–4). The concatenated and 16S BI trees also indicated that, although P. minuta and P. mouritzi showed considerable divergence from each other, they formed a base group in the phylogeny (species group E in Fig. 2 and 3; no COI data are available for the former species). Three species (P. serventyi, P. extracta and P. acidiphila) were each distinctive (see Fig. 2–4).

The morphology of the medial process (MP) or medial space (MS) in the heads of males in species that occurred in the same species groups (A–E) in the molecular phylogenies were usually similar (Table 2, Fig. 5). For example, the six species in group A all had an undivided medial process (Fig. 5). Apart from P. yarleensis, which was not included in the molecular phylogeny, P. serventyi was the only other species to show this feature (Table 2) and, although the P. serventyi haplotype clade was distinctive, it was closest to species group A clade in the concatenated trees (Fig. 2 and Supplementary Fig. S2). Another example is group C, which was formed by P. cylindrifera and Parartemia sp. ‘y’, the only Parartemia species that have a small bifid medial process (Table 2). Parartemia minuta and P. mouritzi in species group E had different MS structures but are the only known Parartemia species possessing basolateral spines (see Fig. 5).

Fig. 5.

Photographs showing the morphology of the medial process (MP) or medial space (MS) in the heads of male Parartemia from 19 species confirmed in this study. The male head of P. minuta has been redrawn from Timms (2014). Coloured bars indicate five species groups (A–E) in Parartemia identified in the molecular phylogenies. The grey-coloured bar indicates the three molecularly distinctive species. The tree topology is derived from the collapsed 16S phylogenetic tree (Fig. 3).

Species delimitation

The three species delimitation methods identified different numbers of partitions among the Parartemia morphospecies. For the 16S dataset, the ASAP scheme with the lowest score (ASAP1; threshold distance 0.01) gave 56 partitions and the ASAP least partition scheme (ASAP2; threshold distance 0.06) and the mPTP and GMYC methods yielded fairly similar results of 27, 30 and 24 partitions respectively (Fig. 3 and Table 3). For the COI dataset, the GMYC (31) and mPTP (38) methods and especially ASAP1 (62; threshold distance 0.01) gave more partitions than ASAP2 (15; threshold distance 0.16) (Fig. 4 and Table 3). As is described below, the ASAP2 scheme showed the closest match with the Parartemia morphospecies. The other methods indicate the presence of one or more partitions within most morphospecies (Table 3). Our species hypotheses mainly conform with the ASAP2 scheme (see Table 3 and below).

For both the 16S and COI genetic regions, the maximum amount of intraspecific genetic distance varied among species and tended to be larger in species that were sampled from a greater number of sites (see Table 4 and Supplementary Table S1). For any one species, the maximum intraspecific genetic distance was always less than its minimum genetic distance from another species, but overall there was a limited amount of overlap between minimum interspecific and maximum intraspecific distances for both the COI and 16S datasets (Table 4, Supplementary Tables S6 and S7).

In total, 8 of the 18 Parartemia morphospecies (16 described and 2 newly identified undescribed species) were each represented by a single well-supported clade (or a distinct haplotype when the species was represented by a single haplotype) that exactly corresponded to a single ASAP2 partition in both 16S and COI datasets (Fig. 3 and 4; also see Table 3). Another two species, P. minuta and P. triquetra (respectively represented by two haplotypes and a single haplotype), each corresponded to a single ASAP2 partition in the 16S dataset but were not included in the COI dataset (Table 3). For this group, genetic distances were highest in P. extracta, which was collected from seven sites plus two 16S sequences from GenBank and had maximum p-distances of 7.53% for 16S and 3.65% for COI. They were also high in P. informis, which was collected from 12 sites and had maximum p-distances of 6.28% for 16S and 12.16% for COI (Table 4 and Supplementary Table S1). The geographic distributions of the divergent 16S lineages within each of these species did not overlap (Supplementary Fig. S5 and S6). This species group also included the putative new morphospecies Parartemia sp. ‘y’, which was collected from two sites and represented by two 16S and four COI haplotypes (Table 3), with maximum p-distances of 0.21 and 1.37% respectively (Table 4). The minimum p-distance between the haplotypes of Parartemia sp. ‘y’ and those of any other species (P. cylindrifera) was 12.16% for 16S and 15.65% for COI (Supplementary Table S6).

The remaining eight morphospecies fell into four categories: (i) morphospecies that formed single ASAP2 partitions in the COI but not the 16S dataset, (ii) morphospecies that formed single ASAP2 partitions in the 16S but not the COI dataset, (iii) multiple morphospecies that were combined into a single ASAP2 partition in both the 16S and COI datasets, and (iv) a single morphospecies that was spilt into multiple ASAP2 partitions in both the 16S and COI datasets.

Category (i) comprised P. cylindrifera and P. serventyi, each of which formed multiple (5) partitions in 16S ASAP2 analysis (Table 3), despite each being represented by a single well-supported clade in the phylogenetic trees (see Fig. 2–4) and only a single partition in the COI ASAP2 analysis. These two species, which were sampled from respectively 12 (plus one COI and one 16S sequence from GenBank) and 13 sites (Supplementary Table S1), had the highest 16S and COI genetic distances of all sampled species (Table 4). Within each species, the distributions of the divergent 16S lineages did not overlap, although those of lineages C and E in P. serventyi bordered each other (see Supplementary Fig. S7 and S8).

Category (ii) included P. laticaudata and P. veronicae, which were combined into a single partition in the COI ASAP2 analysis (Table 3). This was surprising because these two species did not form a monophyletic group in any of the phylogenetic analyses (see Fig. 2–4) and were partitioned via the more conserved 16S region. Although the minimum genetic distances between the haplotypes of these two species, e.g. p-distances of 9.24% for 16S and 14.44% for COI, were among the lowest found in the Parartemia species (see Supplementary Table S6), that between the 16S haplotypes of P. verconicae and P. informis was only 8.39% (Supplementary Table S6). Category (ii) also included the putative new species Parartemia sp. ‘z’, which was combined into a single partition with P. boomeranga and P. longicaudata in the COI analysis; it also grouped with these species in the ML 16S tree (Supplementary Fig. S3). Parartemia sp. ‘z’ was represented by single 16S and COI haplotypes. The minimum p-distances between these haplotypes and those of the P. longicaudata–P. boomeranga morphotypes, at 9.39% for 16S and 13.98% for COI, were low for interspecific comparisons (Supplementary Table S6). Nevertheless, the haplotypes of Parartemia sp. ‘z’ were separated from those of P. longicaudata–P.boomeranga in most of the phylogenetic analyses (see Fig. 2–4) and formed their own partition in all species delimitation analyses except COI ASAP2 (Table 3).

Category (iii) comprised P. boomeranga and P. longicaudata, which in combination corresponded to a single ASAP2 partition in both the 16S (a single partition was also suggested by the 16S GMYC analysis; Fig. 3 and Table 3) and (along with Parartemia sp. ‘z’) in the COI datasets (see above; also see Fig. 4 and Table 3). Together, these species formed a single clade in the phylogenetic trees (which also included Parartemia sp. ‘z’ in the ML 16S tree; Supplementary Fig. S3), but neither morphotype formed an exclusive monophyletic group within this broader clade (see Fig. 2–4). Similarly, although the other species delimitation methods identified multiple partitions within the broader P. boomeranga–P. longicaudata clade, none of these partitions exactly corresponded to haplotypes representing one or the other morphotype (Fig. 3 and 4). The amount of divergence (p-distance) within the P. longicaudata–P. boomeranga complex ranged up to 7.52% for 16S and 11.7% for COI, which was less than that in some other species (Table 4). The P. longicaudata morphotype was common and widely distributed whereas we only found the P. boomeranga morphotype at five sites in the northern Wheatbelt within 125 km of each other (Supplementary Table S1 and Supplementary Fig. S9). The P. boomeranga–P. longicaudata complex included a widespread lineage (E) whose distribution overlapped or bordered that of the other four lineages (A–D) whose distributions did not overlap (see Supplementary Fig. S9).

Parartemia purpurea was the only morphospecies in category (iv). Representatives of this taxon fell into three well-supported clades (a, b and c) in the phylogenetic trees (see Fig. 2–4). Some details of the relationship among these clades varied between datasets and tree-building methods (Fig. 2–4). Regardless, each of the three clades corresponded to separate ASAP2 partitions (also for the mPTP and GMYC methods) in both the 16S and COI datasets (Fig. 3 and 4) and differed from each other by p-distances of at least 13.6 and 18.24% for the 16S and COI data respectively (see details in Supplementary Table S6), which are greater than the maximum intraspecific distance recorded in any other morphospecies (Table 4). On this basis, we propose that ‘P. purpurea’ comprises three cryptic species that herein are called P. purpurea (a), P. purpurea (b) and P. purpurea (c) (see Fig. 2–4). The morphology of some specimens of each species were checked but no characteristic differences were found. All three species were found within 85 km of each other in the Esperance hinterland, but their geographical distributions did not overlap (see Fig. 6).

Fig. 6.

Distribution of three cryptic species of P. purpurea. Site details are in Supplementary Table S1.

Morphology of new species

The overall morphology of the putative new species Parartemia sp. ‘y’ was like that of P. cylindrifera (with which it also had the most molecular similarities). However, the ventral processes in the males of Parartemia sp. ‘y’ featured nearly rectangular distolateral outgrowths and obtuse distomesial corners, whereas those of other species including P. cylindrifera had round distolateral corners protruding ventrally and round distomesial corners (Fig. 7).

Fig. 7.

Photographs of putative new species of Parartemia showing distinctive morphological features relative to morphologically similar species. Comparison of male head of Parartemia sp. ‘y’ (a) with that of P. cylindrifera (b). Comparison of male head and posterior thoracic and anterior abdominal segments of Parartemia sp. ‘z’ (c, d) with those of the P. boomeranga (e, f) and P. longicaudata (g, h) morphotypes. Morphological features mentioned in the text are indicated.

The overall morphology of Parartemia sp. ‘z’ was similar to the P. boomeranga and P. longicaudata morphotypes (with which it also had the most molecular similarities) but the thickened medial edges of the ventral processes in the males were strongly concave (boomerang-shape) contrasting with their weakly concave (boomerang-shape) or cylindrical shape in P. boomeranga–P. longicaudata (Fig. 7). Parartemia sp. ‘z’ was also easily discernible from P. boomeranga–P. longicaudata by the substantially larger lateral lobes on the eleventh thoracic segment in males, which were more than double the size of lateral bulges on the first genital segment (Fig. 7).

Parartemia species POP4 and POP5

To determine the species identity of POP4 and POP5 from Remigio et al. (2001), we trimmed our 16S sequences and compared them to two short 16S sequences (335 and 337 bp) from POP4 and POP5 on GenBank (data not presented). The sequence for POP4 (accession number AY014794) was identical to the corresponding region in one of three P. laticaudata haplotypes. The sequence for POP5 (accession number AY014795) was similar (minimum and maximum p-distances of 2.2 and 4.3% respectively) to the corresponding region in one of five haplotypes found in P. purpurea (c), but highly divergent from the next closest taxa based on 16S, which were P. purpurea (b) and P. purpurea (a) (minimum p-distances of 16.1 and 17.1% respectively). These findings indicate that POP4 corresponds to P. laticaudata and POP5 to P. purpurea (c).

Molecular divergence

The amount of molecular divergence found in Parartemia was large, as was also noted by Remigio et al. (2001). Using COI p-distances as an example, the maximum divergence in Parartemia (25.23%) was higher than that reported for a range of other branchiopod genera, such as Triops (15.6%; Meusel and Schwentner 2017), Limnadopsis (18.6%; Schwentner et al. 2011), Eocyzicus (19.0%; Schwentner et al. 2014), Ozestheria (21.2%; Schwentner et al. 2015) and Artemia (21.8%; Muñoz et al. 2008). It was also greater than that detected in other invertebrate genera that have a long history in Australian salt lakes or their precursors, for example, Coxiella gastropods (18.9%; Lawrie et al. 2023) or Australocypris giant ostracods (18.9%; Rahman 2024). The extreme conditions in salt lakes may accelerate the rate of molecular evolution in halophilic cladocerans and anostracans, like Parartemia (Hebert et al. 2002). Such rate heterogeneity would help to explain why divergence levels in Parartemia are high relative to branchiopods from fresh or low salinity water but not necessarily in comparison with Artemia or other genera of halophilic crustaceans. It may also be that some Parartemia lineages are particularly old or that habitat specialisation has driven divergence among some lineages (see below).

Phylogenetic relationships

Some aspects of the phylogeny of Parartemia, particularly deeper divergences, were not fully resolved in this study. This may be because the selected genetic markers were not suitable for higher-level taxonomic resolution in Parartemia or because these deeper divergences may be associated with rapid cladogenesis (e.g. see Kang et al. 2008; Whitfield and Kjer 2008; Pinceel et al. 2013). However, some groups of related species and some distinctive species were evident. The fact that species thought to be closely related based on the molecular data usually also showed similarities in the structures of their medial processes or medial space in the males adds evidence that these groups are real. Like Remigio et al. (2001), we found that P. minuta (which occurs in central and eastern Australia) was usually positioned at the base of the phylogeny and, in our study, usually grouped with P. mouritzi (which occurs in Western Australia and was not sampled by Remigio et al. 2001). These are the only two Parartemia species known to possess basolateral spines (see Timms 2012), providing further evidence that they are sister species. Like Remigio et al. (2001), we found that P. extracta was a distinctive species. Remigio et al. (2001) found two new species of Parartemia, which they called POP4 and POP5 and that we have identified as P. laticaudata and P. purpurea (c) respectively. The species groups that Remigio et al. (2001) called clades A and B were not apparent in our samples. For example, that study placed P. longicaudata (which included their POP1, POP3 and POP6 samples), P. cylindrifera, P. purpurea (c) (then POP5) and P. zietziana in clade A whereas our study usually placed P. purpurea (c) and P. zietziana in the same clade but P. longicaudata in a different clade and P. cylindrifera in yet another one. The discrepancies between the results of this study and those of Remigio et al. (2001) can be attributed to the results of the latter study being based on only a short 16S fragment and only a small subset of species (see Introduction). Increased taxonomic representation has been shown to alter perceptions of species relationships in a range of taxa, including Branchinella fairy shrimps in Australia (Pinceel et al. 2013).

On the basis that new anostracans species are expected to arise from widespread species in peripheral habitats, Rogers (2015) predicted that anostracan genera will comprise a series of small clades, each containing a basal widespread species and a derived species with a narrow distribution at the periphery of the widespread one (see also Rogers and Aguilar 2020). Our phylogenies provide several examples of widespread and narrowly distributed sister species including (with the widespread species listed first) – P. cylindrifera and Parartemia sp. ‘y’, P. longicaudata–P. boomeranga and Parartemia sp. ‘z’, P. informis and P. triquerta and P. minuta and P. mouritzi (see Timms et al. 2009 and below for distributional information). However, other aspects of the relationships between these sister species do not fit with the above prediction. For example, relative to their widespread counterparts, the narrow range species were usually ancestral in the phylogenies or had either overlapping (e.g. P. cylindrifera and Parartemia sp. ‘y’) or disjunct distributions (e.g. P. informis and P. triquerta). Intraspecific divergence within some widespread Parartemia species may more closely conform with the above prediction and will be examined in detail in a future study on phylogeographic patterns in Parartemia.

Species delimitation

The most recent previous estimate of the number of Parartemia morphospecies was 18 (Timms 2010). Our study suggests that Parartemia consists of at least 21 species. These species comprise: (i) 13 described morphospecies that have been confirmed using the molecular data in this study; (ii) 2 putative new morphospecies; (iii) 3 cryptic species of P. purpurea; (iv) 1 species encompassing the P. longicaudata and P. boomeranga morphotypes; and (v) 2 described morphospecies that were not included in this study. These data suggest that Parartemia is the most speciose genus of halophilic invertebrates found in Australia, notwithstanding that the taxonomy of some of the other invertebrate groups is rudimentary (see Lawrie et al. 2021). Based on the currently available data, the next most diverse genus is Coxiella (a gastropod) with at least 15 species, although recent molecular evidence indicates that these species are spread over several unrecognised genera (four clades in Lawrie et al. 2023). Australocypris with 10 species is the most species-rich genus in Mytilocypridinae giant ostracods (Rahman 2024). Parartemia is much more species-rich than Artemia, which only has nine recognised species even though Artemia is essentially globally distributed (Rogers 2013; Asem et al. 2023b). Compared to Artemia, speciation in Parartemia is likely facilitated by their heavy, sinking resting eggs, which tend to retard dispersal (McMaster et al. 2007; Timms et al. 2009). Australia also has a diverse range of Branchinella fairy shrimps in fresh or low-salinity water (Pinceel et al. 2013; Rogers and Timms 2014). This diversity and that in Parartemia is probably linked to a long history of aridity and persistent ephemeral water bodies in the Australian landscape (Rogers and Timms 2014). Rogers (2015) and references therein have argued that species diversity in anostracans will be enhanced in an older landscape like Australia where more occupied habitats will favour speciation by habitat specialisation over colonisation of vacant habitats. This fits with the high levels of ecological specialisation apparent among Parartemia species, e.g. with regard to pH or substrate geochemistry (see Timms et al. 2009; Timms 2012).

Our species hypotheses for Parartemia are based on the most conservative results from the species delimitation analyses (i.e. the ASAP scheme with the least number of partitions), which tended to match morphological species boundaries. This fits with the suggestion that morphological diagnosis of anostracan species is generally straightforward (Rogers and Aguilar 2020), although this is not always the case (e.g. see Ketmaier et al. 2008; Pinceel et al. 2013; Rogers and Aguilar 2020). Many Parartemia morphospecies comprised multiple divergent lineages that were partitioned in some of the species delimitation analyses. Other than for the three lineages of P. purpurea (discussed below), we have interpreted the divergent lineages within morphospecies as examples of intraspecific variation rather than as cryptic species for the following reasons. (1) The presence of morphologically similar but genetically divergent lineages may be linked to an accelerated rate of molecular evolution in Parartemia (see above). (2) With the exception of P. longicaudata, the geographic distributions of the divergent lineages of morphospecies were essentially nonoverlapping. Also, conspecific divergent lineages never co-occurred in the same water body (see Supplementary Fig. S5–S9). (3) Species delimitation methods are sensitive to the population structure of a species and sometimes delimit genetically divergent populations within a species (Luo et al. 2018; Gaytán et al. 2020). (4) Defining species based on molecular data alone is challenging (Jörger and Schrödl 2013; Fišer et al. 2018). (5) We wanted to avoid the pitfalls of taxonomic inflation (see Padial and De la Riva 2006).

We found two putative new morphospecies of Parartemia – Parartemia sp. ‘y’ and Parartemia sp. ‘z’ – in our samples. These taxa exhibited both distinctive morphological and genetic characteristics relative to other described species. The evidence is stronger for Parartemia sp. ‘y’, which was found at two sites, than for Parartemia sp. ‘z’, which was only found at a single site. The results of this study provide the first evidence of cryptic species in Parartemia, namely P. purpurea (a), P. purpurea (b) and P. purpurea (c), although cryptic species have been reported in a range of anostracans (e.g. Pinceel et al. 2013; Rogers 2014) and other branchiopods (e.g. Schwentner et al. 2013; Meusel and Schwentner 2017). The three P. purpurea lineages were consistently separated from each other in the species delimitation analyses and did not always form a single well-supported clade in the phylogenetic trees. No morphological differences were detected among individuals from these lineages, including in relation to characters important in species diagnosis in Parartemia (see ‘Introduction’), although it is possible that further scrutiny could reveal some subtle differences. Ecological or physiological specialisation may have been crucial in driving cryptic speciation in P. purpurea (e.g. Rogers 2014).

The two new Parartemia morphospecies as well as the three cryptic species of P. purpurea were found in only the Esperance hinterland region of Western Australia. This brings the total number of Parartemia species found in this region to 9 (out of 21), 5 of which have only been found in this region (this study and Timms et al. 2009). In general, the region is well known for the diverse range of invertebrates found in its salt lakes (see Timms 2009). Coxiella gastropods (Lawrie et al. 2023) and Australocypris ostracods (Rahman et al. 2023) also have high species richness and multiple endemic species in this region. The region hosts a large number and variety of salt lakes (see Timms 2009 for details), which are probably important factors driving divergence and diversity. It also seems likely that this region has served as a refugium in evolutionary time, limiting extinctions, possibly because of favourable hydrological or climatic conditions (e.g. see Jansson 2003; Davis et al. 2013).

Our results suggest that the P. boomeranga and P. longicaudata morphotypes comprise a single monophyletic lineage. These results are based on a comprehensive sampling of these morphotypes. Parartemia boomeranga was sampled across its documented range, including three lakes near the type locality (unnamed lake east of Gunyidi, −30.12, 116.24; see Timms 2010). Parartemia longicaudata was also sampled across its entire known geographical range, from Esperance to the Houtman Abrolhos Islands, including the neotype locality (Pink Lake in Esperance; see Supplementary Fig. S9 and Timms 2010). Morphological differences between P. boomeranga and P. longicaudata are minor and confined to the thickened medial edges of the ventral processes in males (boomerang-shaped in the former and cylindrical in the latter) and the shape of the medial space (concave in the former and convex in the latter) (see Fig. 7 and Timms 2010). Regardless, the P. boomeranga and P. longicaudata morphotypes did not form reciprocally monophyletic groups in the molecular phylogeny, so each morphotype must have evolved more than once or may be environmentally induced. Following the guidelines set by the International Code of Zoological Nomenclature (International Trust for Zoological Nomenclature 1999), P. longicaudata being named first has nomenclatural priority over P. boomeranga.

Conservation implications

Timms et al. (2009) provided a detailed assessment of the conservation status of Parartemia species. This assessment can be updated using the improved taxonomic and associated distributional information on Parartemia generated by this study.

Timms et al. (2009) proposed that P. boomeranga (then Parartemia sp. ‘c’) should be assessed as potentially vulnerable. Timms (2012) later advised that this species was extremely rare and at risk of extinction. However, we encountered this species at five sites within the reported range of this species, including Lake Moore (in the northern wheatbelt region in Western Australia) from which this taxon has previously been recorded (Timms et al. 2009) and three sites that are near the type locality (see above). Furthermore, the molecular results suggest that the P. boomeranga morphotype is not a valid species but is synonymous with the common and widely distributed P. longicaudata, which is not regarded as threatened (see Timms et al. 2009).

The two new species and the three cryptic species of P. purpurea discovered in this study all have narrow geographic distributions. Except for P. purpurea (c), these species are known from between only one and three sites in the Esperance hinterland region. All these species have sites in the Kau Rock Nature Reserve or the nearby Beaumont Nature Reserve and are therefore offered some protection but are nonetheless somewhat vulnerable given their rarity and restricted distributions. Parartemia purpurea (c) is known from five or six sites (depending on whether we resampled Remigio et al.’s (2001) POP5 collection site, the exact location of which was not reported) but all are in the Esperance hinterland and none are in nature reserves. The above distributional information does not consider 20 other sites in the Esperance hinterland that were not sampled in this study but are reported to contain the P. purpurea morphotype (see Timms et al. 2009; Timms 2010; Rogers and Timms 2014). Regardless, as an area of evolutionary importance with high species diversity (see above and Timms et al. 2009), the Esperance hinterland should be a priority site for Parartemia conservation (see Moritz 2002; Davis et al. 2013).

Limitations and future work

Future phylogenetic studies on Parartemia should aim to include P. auriciforma and P. yarleensis, which were missing from this study, as well as more samples from salt lakes in remote areas, which have typically been poorly studied and may harbour a high proportion of undiscovered species (Timms 2010). It is also important to formally describe and name the new species discovered in this study so they can be taken into account in conservation planning and legislation (e.g. see Mace 2004; Padial and De la Riva 2006). A better understanding of the evolutionary significance of the divergent lineages within morphospecies is also needed. This could be facilitated via assessments of phylogeographic structures of these species (e.g. see Seidel et al. 2009).

The phylogenetic results of this study are based on a total of three different genetic markers, which were not sufficient to resolve the relationships among some lineages, particularly those reflecting deeper divergences. Future studies should include more markers and loci to improve the resolution of phylogenetic relationships in Parartemia, although there is no simple answer as to the minimum number of markers or loci needed to produce a robust phylogeny or whether markers should be selected at random or systematically (Gatesy et al. 2007).