The molecular dataset was complemented with sequences retrieved from GenBank corresponding to a selection of 134 Niphargus species/populations covering the entire geographical distribution and morphological disparity of the genus, plus four of the nine additional recognized niphargid genera (v. Horton and Lowry, 2013), viz. Carinurella Sket, 1971, Pontoniphargus Dancău, 1968, Niphargobates Sket, 1981 and Niphargellus Schellenberg, 1938. We completed the dataset with a series of more distant outgroups in the families Gammaridae (6 genera), Gammaracanthidae (1), Eulimnogammaridae (1), Acanthogammaridae (1), Pallaseidae (1), Pontogammaridae (1), Anisogammaridae (1) and Crangonyctidae (3). We rooted the tree with a phoxocephalid and an urothoid, in accord with the most recent hypothesis on amphipod supra-family-level relationships (Lowry and Myers, 2013). Data on all sequences and taxa used are shown in Table S1.
Multiple sequence alignment was performed using MAFFT 7 online version (http://mafft.cbrc.jp/alignment/server/; Katoh and Standley, 2013), using the default FFT-NS-1 algorithm for the cox1 and H3 datasets, and the Q-INS-i option for the LSU sequences since it considers the secondary structure of the RNA sequences. In order to evaluate the effect on the topology of the indels in the LSU alignment, an additional LSU matrix was generated by removing ambiguously aligned regions with Gblocks v. 0.91b (Talavera and Castresana, 2007), and specifying a minimum length of two positions for a block and selecting only positions with a gap in less than 50% of the sequences if they are within an appropriate block. Nucleotide substitution saturation was assessed with the Xia’s test (Xia et al., 2003) implemented in DAMBE v. 5.2.64 (Xia, 2013). Since evidence of substitution saturation was detected on third coding positions of the cox1 alignment (see Results), we also performed the phylogenetic analyses by excluding these characters from the dataset. Genetic divergence among Haploginglymus taxa was estimated through the calculation of the uncorrected pairwise cox1 distances (p-distance) in MEGA7 (Kumar et al., 2016). In order to test the phylogenetic congruence among the three molecular markers, single-gene trees were inferred and the resulting topologies were compared. Optimal partitioning strategy and evolutionary models for each alignment were assessed with PartitionFinder (Lanfear et al., 2012) under the Bayesian Information Criterion (BIC), whereas the phylogenetic inference was conducted using IQTREE multicore v. 1.3.12 (Nguyen et al., 2015) performing 1000 ultrafast bootstrap approximation replicates. Finally, a phylogenetic analysis of the three molecular markers combined was independently carried out under both a maximum likelihood and a Bayesian framework using IQTREE and MrBayes 3.2 (Ronquist et al., 2012), respectively. For the latter, two independent analyses consisting of four chains each were run for 5·107 generations specifying a sampling frequency every 1000 generations, and setting up a burn-in fraction of 35%. MCMC convergence and effective sample size (ESS) estimates were checked with TRACER v. 1.6 (Rambaut et al., 2013).