Phylogenetic analysis
We analysed only populations with at least three out of four molecular markers sequenced; thus some of the literature data and some sequences available at GenBank were not included. The available Milnesiumsequences, together with those from outgroup taxa (see Table S1–S2), were aligned independently for each marker, using MAFFT version 7 (Katoh, Misawa, Kuma, & Miyata, 2002; Katoh & Toh, 2008). For 18S rRNA and 28S rRNA, the G-INS-I strategy was applied to consider the secondary structure of RNA, whereas for ITS-2 and COI, the default setting was used. The obtained alignments were checked manually in BioEdit and then trimmed to 1159 (18S rRNA), 887 (28S rRNA), 841 (ITS-2) and 569 bp (COI). Next, the four obtained alignments were concatenated in SequenceMatrix (Vaidya, Lohman, & Meier, 2011). We used the Bayesian Information Criterion in PartitionFinder version 2.1.1 (Lanfear, Hua, & Warren, 2016) to find the most suitable substitution model for posterior phylogenetic analysis. As COI is a protein coding fragment, the alignment was previously divided into three data blocks corresponding to the three codon positions. The best fit-model was GTR+I+G for all the partitions, with the exception of the first and third codon positions in COI, for which SYM+I+G and K81UF+G had a better fit respectively.
Phylogenetic inference was carried out using two different approaches. The first one utilized the same exact methodology as Morek & Michalczyk (2020) in MrBayes v3.2 (Ronquist & Huelsenbeck 2003). Additional Bayesian Inference analyses were carried out in BEAST 1.10.4 (Drummond, Suchard, Xie, & Rambaut, 2012) using different clock (strict, lognormal relaxed, and exponential relaxed) and tree (coalescent with constant population size versus Yule) models. Different nucleotide substitution models and clock rates were set for each partition. In order to calibrate the molecular clock, and in the absence of any other reliable calibration, we used two calibration points: (i ) we constrained the divergence between Hetero- and Eutardigrada using a normal distribution according to the estimates by Guidetti et al. (2017), and (ii ) we set a lognormal distribution for the origin ofMilnesium using the estimated age of the only know fossil assigned to this genus, Milnesium swolenskyi Bertolani & Grimaldi, 2000 as lower bound.
The analyses ran for 100 million generations in the CIPRES Science Gateway (Miller, Pfeiffer & Schwartz, 2010), sampling every 10,000 steps. The best combination of clock and tree priors, was relaxed exponential clock model with Yule process as tree prior, which were selected according to the Bayes factors calculated in TRACER 1.6 (available from http://beast.bio.ed.ac.uk). The consensus tree for this best combination was built using TREEANNOTATOR (distributed with BEAST). Additionally, mutation rates for the four utilised markers, estimated in TRACER, are given in Table 2.
In all the trees, clades recovered with posterior probability (PP) between 0.95 and 1.00 were considered well supported, those with PP between 0.90 and 0.94 were considered moderately supported and those with lower PP were considered poorly supported or unsupported. The consensus tree was visualised in FigTree v.1.4.3 (available from http://tree.bio.ed.ac.uk/software/figtree).