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).