Morphological and genetic variability in cosmopolitan tardigrade species—Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf, 2010

Paramacrobiotus fairbanksi was described from Alaska (USA) based on integrative taxonomy and later reported from various geographical localities making it a true cosmopolitan species. The ‘Everything is Everywhere’ (EiE) hypothesis assumes that the geographic distribution of microscopic organisms is not limited by dispersal but by local environmental conditions, making them potentially cosmopolitan. In the present work we report four new populations of P. fairbanksi from the Northern Hemisphere which suggests that the ‘EiE’ hypothesis is true, at least for some tardigrade species. We also compared all known populations of P. fairbanksi at the genetic and morphological levels. The p-distances between COI haplotypes of all sequenced P. fairbanksi populations from Albania, Antarctica, Canada, Italy, Madeira, Mongolia, Spain, USA and Poland ranged from 0.002 to 0.005%. In total, twelve haplotypes (H1-H12) of COI gene fragments were identified. We also report statistically significant morphometrical differences of species even though they were cultured and bred in the same laboratory conditions. Furthermore, we also discuss differences in the potential distribution of two Paramacrobiotus species.

Paramacrobiotus fairbanksi was described from Alaska (USA) based on integrative taxonomy and later reported from various geographical localities making it a true cosmopolitan species.The 'Everything is Everywhere' (EiE) hypothesis assumes that the geographic distribution of microscopic organisms is not limited by dispersal but by local environmental conditions, making them potentially cosmopolitan.In the present work we report four new populations of P. fairbanksi from the Northern Hemisphere which suggests that the 'EiE' hypothesis is true, at least for some tardigrade species.We also compared all known populations of P. fairbanksi at the genetic and morphological levels.The p-distances between COI haplotypes of all sequenced P. fairbanksi populations from Albania, Antarctica, Canada, Italy, Madeira, Mongolia, Spain, USA and Poland ranged from 0.002 to 0.005%.In total, twelve haplotypes (H1-H12) of COI gene fragments were identified.We also report statistically significant morphometrical differences of species even though they were cultured and bred in the same laboratory conditions.Furthermore, we also discuss differences in the potential distribution of two Paramacrobiotus species.
The Phylum Tardigrada currently consists of ca.1500 species 1 that inhabit terrestrial and aquatic environments throughout the world 2 .Currently there are 33 families, 159 genera, 1464 species and 21 additional subspecies within this phylum 1 .Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf, 2010 3 was described from Alaska (USA) and reported from the Antarctic, Italy, Poland and Spain 4 (reported as Macrobiotus richtersi Murray, 1911 5 ) [6][7][8][9] .It is a large-size (up to 800 μm) parthenogenetic Paramacrobiotus found mostly in mosses and can be shortly characterized by white or transparent cuticle without pores, three bands of teeth in the oral cavity, three macroplacoids and a microplacoid in pharynx (richtersi group), smooth lunules under all claws, granulation on all legs, and eggs with reticulated conical processes without caps or spines.Paramacrobiotus fairbanksi is a triploid species 8 inhabiting various locations throughout the globe.The species is an omnivore, i.e., it feeds on algae, cyanobacteria, fungi, nematodes and rotifer 10 .However, dietary preferences have been observed to differ between juveniles and adults (juveniles prefer green alga and adults favour rotifers and nematodes 10 ).
The 'Everything is Everywhere' hypothesis, which was proposed at the beginning of the twentieth century 11,12 suggests that microorganisms and small invertebrate are not dispersal-limited on large geographical scales and have potentially cosmopolitan distribution.Microscopic organisms are often considered cosmopolitan species, as, the presence of specific adaptations allows them to being disperse easily.These adaptations include (a) the possibility of easy passive dispersion (by wind, rivers, sea currents, other animals, etc.), (b) the presence of very

Culture procedure
Specimens of the populations from Albania, Canada, Madeira and Mongolia were cultured in the Department of Animal Taxonomy and Ecology (Faculty of Biology, Adam Mickiewicz University in Poznań) according to the protocol described by Roszkowska et al. 41 .In summary, tardigrades were cultured in small Petri dishes in spring water mixed with distilled water (1:3) with the rotifers and nematodes added as food ad libitum.All cultures were kept in the environmental chamber at a temperature of 18 °C and in darkness.

Microscopy, morphometrics and morphological nomenclature
Specimens were extracted from cultures and prepared for light microscopy analysis.They were mounted on microscope slides in a small drop of Hoyer's medium and secured with a cover slip 42,43 .Slides were then placed in an incubator and dried for two days at ca. 60 °C.Dried slides were sealed with transparent nail polish and examined under an Olympus BX41.

Molecular data analysis
The amplified nuclear and mitochondrial barcode sequences were edited using the BioEdit software 53 .Comparison of obtained molecular markers with those deposited in GenBank and homology search were performed using Blast application (Basic Local Alignment Search Tool 54 ).The COI haplotypes were generated using the DnaSP v5.10.01 program 55 and were translated into amino acid sequences using the EMBOSS-TRANSEQ application 56 to check for internal stop codons and indels.Then all sequences obtained in our study, and the sequences downloaded from the GenBank database as originating from P. fairbanksi, were aligned with ClustalW using default settings.Alignment sequences were trimmed to 689, 572 and 574 bp for 28S rRNA, 18S rRNA and COI barcodes, respectively.The calculation for the uncorrected pairwise distances (p-distances) was performed for COI sequences using the MEGA X 57 .
All obtained sequences have been deposited in GenBank (for the accession numbers please see Table 3).The slides prepared from exoskeleton/voucher after DNA extraction of P. fairbanksi were deposited at the Department of Animal Taxonomy and Ecology, Institute of Environmental Biology, Adam Mickiewicz University, Poznań, Uniwersytetu Poznańskiego 6, 61-614 Poznań, Poland.www.nature.com/scientificreports/Reconstruction of genetic relationships among COI haplotypes and genealogical connections was carried out using the Network 4.6.1.3software (www.fluxu xengi neeri ng.com).The median-joining algorithm (MJ) 58 and substitution rates with the weight of 3 for transitions and 1 for transversions (transition: transversion ratio (ti:tv)) were applied.The star contraction pre-processing was generated to delete all superfluous median vectors and links.Additionally, the maximum parsimony post-processing was calculated.In turn, signatures of population expansion, equilibrium or decline in P. fairbanksi were inferred from the neutrality tests calculation (Tajima D 59 and Fu F S 60 , respectively) computed in the DnaSP v5.10.01 program and Arlequin v.3.5.software 61 .Analyses were performed with 1000 replicates.

Statistical analysis
We used the Analysis of Variance (ANOVA) test with post hoc comparison of pairs of measurements, applying Benjamini-Hochberg correction to statistically analyze the differences in morphometrics between different populations of P. fairbanksi.Measurements of the body and buccal tube length (BL and BTL, respectively) were used as the dependent and the populations as grouping variables.Normal distribution in residuals was checked using the Shapiro test.Other morphometric traits, i.e., stylet support insertion points (SSIP), external width of buccal tube (BTEW) and placoids (M1-macroplacoid 1, M2-macroplacoid 2, M3-macroplacoid 3, Mi-microplacoid, MR-macroplacoid row, PR-placoid row) were also analysed.All the analyses were performed in R 4.1.3 62.The level of statistical significance was considered at p < 0.05.Principal Component Analysis (PCA) was performed using the R script from Stec et al. 63 .The analysis was performed for data from eggs and animals.For animals, both absolute values (raw measurements in μm) (BLm, BTLm, SSIPm, BTEWm, M1m, M2m, M3m, Mim, MRm and PRm) and relative pt values (BLpt, SSIPpt, BTEWpt, M1pt, M2pt, M3pt, Mipt, MRpt and PRpt) were used.For eggs, absolute values (raw measurements in μm) were used.All analyses were carried out using the R software program 48 .The "imputePCA" function of the R package "missMDA ver.1.17" was used to impute missing data in the animal data set using the PCA imputation technique 64 .Cross-validation (function "estim_ncpPCA") was used to determine the number of components utilized to impute the missing data.The PCA function of the software "FactoMineR ver.2.3" 65 was used to perform PCAs on the scaled data.The software "ggplot2 ver.3.3.2","plyr ver.1.8.6", and "gridExtra ver.2.3" were used to depict PCAs 66,67 .The presence of a structure in the PCA data was tested using a randomization approach on the eigenvalues and statistics according to Björklund 68 and an in-house R script developed by MV in Stec et al. 63 .PERMANOVA was done on the PCs with the R packages "vegan ver.2.5.6" and "pairwiseAdonis ver.0.3" 69 , with the species hypothesis generated by phylogenetic techniques as the independent variable.Using the Benjamini-Hochberg correction, the level for multiple post hoc comparisons was adjusted independently for adults and eggs 70 .In total, 106 tardigrade specimens (16 Albanian, 16 Antarctic, 17 Canadian, 15 Madeiran, 14 Mongolian, 15 Polish, 4 Italian and 9 Alaskan) were measured and later used in the analyses for animals.Furthermore, differences in egg morphology between populations were studied and tested using ANOVA.Egg bare diameter (EBD), full diameter (EFD) and processes height (PH) were characters for the populations used as the dependent variable to determine compared groups and Benjamini-Hochberg corrections.In total, 100 tardigrade eggs (15 Albanian, 16 Antarctic, 15 Canadian, 15 Madeiran, 6 Mongolian, 15 Polish and 18 Alaskan) were measured and used in the analyses.All the analyses were performed in R 4.1.0.The level of statistical significance was considered at p < 0.05.An ecological niche modelling (ENM) approach was used to predict the current potential distribution of P. fairbanksi and P. gadabouti.The ENM was performed with the use of the MaxEnt software, ver.3.2.0.( https:// biodi versi tyinf ormat ics.amnh.org/ open_ source/ maxent/).MaxEnt performs the model with the fewest possible occurrence data and takes presence-only (PO) data.The model generates models of habitat appropriateness by handling continuous and categorical variables using regularization parameters 71,72 .The raster package in R was used to extract climatic raster values, and for ENM evaluation, version 0.3.1 of ENMevaluate in R was used.The bioclimatic variables available in MERRAclim Dataset 19 were used as environmental variables for MaxEnt modelling.We used MERRAclim Dataset because it provides a global set of satellite-based bioclimatic variables that includes Antarctica, which is one of the localities for P. fairbanksi.The 19 global bioclimatic datasets from the 2000s at 5 arcminutes resolution (mean value) 73 consist of temperature layers (BIO1-BIO11) and humidity layers (BIO12-BIO19).Using ENMTools 74 correlations among environmental layers were tested for all 19 variables and correlation coefficients >|0.7| were removed.Based on correlation results, six bioclimatic variables with correlation coefficients <|0.7| were selected.The temperature layers are in degrees Celsius multiplied by 10 and the humidity layers are in kg of water/kg of air multiplied by 100,000 73 .The receiver operating characteristic (ROC) plot's area under curve (AUC) was used to assess the model's accuracy 71 .AUC describes the relationship between www.nature.com/scientificreports/ the proportion of correctly anticipated presences and the proportion of absences of mistakenly-projected species in the model 75 .The AUC gauges the effectiveness of the model with a value between 0 and 1.Furthermore, AUC values > 0.9 indicate excellent accuracy, 0.7 to 0.9 indicate good accuracy, and values below 0.7 indicate low accuracy 71,76,77 .The jackknife test was used to estimate the model's variable relevance.Later, the MaxEnt results were masked to the landmasses using shapefile from https:// www.natur alear thdata.com/ downl oads/ 10m-physi cal-vecto rs/ 10m-land/.The localities for P. fairbanksi are from Table 1 and P. gadabouti from Kayastha et al. 22 .The coordinate list is provided in SM.01 and the R script for ENM in SM.02.The html MaxEnt output provided as SM.03a and SM.03b.

Ethics declarations
All procedures were conducted in accordance with the guidelines.Also, none of the moss samples were collected from the region which requires permission.
The randomization test in PCA demonstrated that only the first two PCs explained greater variation than the null model (no data structure) for both animal and egg datasets (SM.11).As a result, only the initial two PCs were maintained and used for additional investigation and interpretation.Furthermore, the ψ and ϕ statistics of the PCA were significantly distinct from what they anticipated under the null assumption (animals: ψ = 60.72 www.nature.com/scientificreports/p < 0.001, ϕ = 0.82 p < 0.001; animals pt: ψ = 13.14 p < 0.001, ϕ = 0.43 p < 0.001; eggs: ψ = 17.62 p < 0.001, ϕ = 0.50 p < 0.001).The first two components of the PCA of animals' absolute measured value (Fig. 5) explained 90% of the overall variation (83.7% for PC1 and 6.7% for PC2) and for animals' pt indices (Fig. 6) explained 65% of the overall variation (46.3% for PC1 and 18.7% for PC2).PCA of egg measurements (Fig. 7) described 68% of the total variance with the first two components (52.5% for PC1 and 15.5% for PC2).PERMANOVA revealed that species identity has a substantial overall effect on PCs (p < 0.001, Tables 12, 13, 14).Raw morphometric data for all the populations in the present study are given in the Supplementary Materials (SM.04-SM.07).R script for single characters as well as measurement files for both adults and eggs, are provided in the Supplementary Materials (SM.08,SM.09a and SM.09b).All the test results from R are provided in Supplementary Materials (SM.10).Results of PCA randomization tests in the Supplementary Materials (SM.11).www.nature.com/scientificreports/

Genetic comparisons and phylogeographical analyses of different populations of the P. fairbanksi
The COI sequences of P. fairbanksi from Albania, Canada, Madeira and Mongolia were 623-689 bp-long, and represented three haplotypes: haplotype H11 was observed in the population from Albania, haplotype H1 was identified in P. fairbanksi from Mongolia, and haplotype H4 was found in populations from Canada and Madeira (for details see Table 3 and Fig. 8A,B).No stop codons, insertions or deletions were observed.The translation was successfully carried out with the -2nd reading frame and the invertebrate mitochondrial codon table.The p-distances between COI haplotypes of all sequenced P. fairbanksi populations deposited in GenBank, i.e., from Antarctica, Italy, Spain, the USA and Poland ranged from 0.002 to 0.005% (an average distance of 0.003%) (Fig. 8B).In total, twelve haplotypes (H1-H12) of COI gene fragments were identified after comparing all available COI sequences of P. fairbanksi.Overall, the median joining COI haplotype network showed a star-like radiation.Interestingly, the most frequent haplotype H4 was present in populations from Italy, Madeira and Canada.This central haplotype H4 was surrounded by ten haplotypes (H1, H3, H5-H12) that differed from it by one mutational step.One haplotype (haplotype H2 from Spain) differed from central haplotype H4 by two mutational steps.In several geographical regions, i.e., the USA, Albania, Italy, Poland and Spain there were regional endemic haplotypes.Surprisingly, the second haplotype that occurred in different localities was haplotype H1 and this haplotype was common for three populations, from Mongolia, Poland and Antarctica.
In turn, the 18S rRNA sequences of P. fairbanksi from Albania, Canada, Madeira and Mongolia were 917-1547 bp-long (Table 3) and no nucleotide substitution was found (except for a single undetermined "N" base).Compared with the data available in GenBank sequences of P. fairbanksi (sequences were aligned and trimmed to 572 bp), they showed only one nucleotide substitution.A comparison was performed with the sequences from the following geographical localities: Antarctica (GenBank: MN960302 9 ), Poland (GenBank: MH664941-42 78 ), USA (GenBank: EU038078 79 ) and Italy (GenBank: MK041027-29 8 ).The 28S rRNA molecular marker (694-805 bp-long) was very conservative.No nucleotide substitution was found for all obtained sequences even after comparing (and trimmed to 689 bp) with GenBank sequences from Antarctica (GenBank: www.nature.com/scientificreports/MN960306-MN960307 9 ) and Poland (GenBank: MH664950 78 ).Nevertheless, one unidentified base was found in the sequence originating from the Polish population.Demographic expansion was preliminarily tested based on the value of neutrality tests that confirmed a neutral model of observed polymorphism.Negative significant values for Tajima's D were found, indicating a high number of low-frequency polymorphisms in the COI sequences dataset and potential population size expansion (Fig. 8C).In turn, values of Fu's Fs test statistic for COI data were negative, but non-significant: − 7.794, 0.10 > P > 0.05 (graphical results not shown).

Predictions of the distribution of the two parthenogenetic Paramacrobiotus species
Ecological niche modelling of potential distribution based on available location data was performed for two parthenogenetic species with verified records from various realms, i.e., P. fairbanksi and P. gadabouti.The study is limited to bioclimatic variables.The stimulated model predicted good accuracy for the overall model with an AUC for P. fairbanksi of 0.778 and the overall model with an AUC for P. gadabouti of 0.843.The suitability for P. fairbanksi seems moderate (peach areas on the map in Fig. 9A) to good (yellow areas on the map in Fig. 9A) with the most suitable habitats in the northern hemisphere (green areas on the map in Fig. 9A).Paramacrobiotus gadabouti shows maximal suitability around areas with a Mediterranean climate, although it also has a wider distribution than just the Mediterranean biomes (tropics, subtropics, and temperate regions) (Fig. 9B).

Morphometric comparison of different populations of the P. fairbanksi
Based on morphometric analyses, there is clear variation in measurements of morphological features between populations of P. fairbanksi from different regions of the world.However, the identification of this species is still possible based on morphological characters, in fact the morphometric data are shown to not affect the species identification because of the overlap in measurements of all measured structures.Therefore, it is valid to suggest the correct classification of all the specimens collected from different regions based on their morphology only.Even though the egg processes of Polish and Albanian populations are similar, the EBD of the Polish population www.nature.com/scientificreports/are the and those from the Albania are largest.The EFD as well as egg processes of the Madeiran population are largest while those of the Polish population are the smallest.Additionally, body length values of the Polish and USA populations of P. fairbanksi are smaller compared to the other populations studied.Kaczmarek et al. 9 suggested that the differences in measurements between different populations of this species are caused by conditions, i.e., specimens from cultures and specimens from wild populations.However, in the present study all measurements were based on specimens from cultured populations, i.e., Albanian, Canadian, Madeiran and Mongolian.Thus, we can suggest that the phenomenon described by Kaczmarek et al. 9 (that dwarfing is caused by suboptimal conditions, high culture densities and inbreeding and that it might be due to the result of ongoing speciation) is unlikely to be true.Similarly, the suggestions that harsh conditions in Antarctica may favour laying larger eggs while in cultures the eggs are smaller because of the lack of such selective pressure 9 seems untrue as the egg size of specimens from Antarctica overlaps with egg sizes of specimens from Albania, Canada and Mongolia, which were sampled from cultured populations in the present study.

Genetic comparison of different populations of the P. fairbanksi
Cytochrome oxidase subunit I gene (COI) sequences is one of the most reliable barcodes to investigate genetic variation with phenotypic plasticity since COI is a genetic marker with a high genetic variation compared to multiple other DNA barcodes 80 .Various studies combining COI variation and phenotypic plasticity were conducted throughout different invertebrates' phyla, including tardigrades 9,21,[81][82][83][84] , proving the marker's accuracy in this group of organisms.The result showed high genetic homogeneity between organisms with wide geographical distribution together with clearly visible morphological differences known as phenotypic plasticity 82 .This phenomenon could be possible explanation for the subtle morphological variation observed in this study between P. fairbanksi populations, accompanied by very low mitochondrial barcode variation.
Furthermore, several studies uncovered data incongruence between mitochondrial and nuclear markers, e.g., for earthworms 85 or corals 86 , suggesting that occasionally COI may fail as a barcode marker due to hybridization events.Many studies have already shown that Wolbachia (presence shown by Mioduchowska et al 87 in P. www.nature.com/scientificreports/fairbanksi) can increase the speciation rate and can affect COI haplotypes 88 .However, the nuclear markers tested for P. fairbanksi have been consistent for the studied The exact causes and mechanisms of the phenotypic plasticity in the morphology of adults and eggs of P. fairbanksi remains unknown, although, it has unsurprisingly been shown, that some physical traits differ in chosen cultured tardigrades depending on the temperature and food abundance 7 .If the morphological variation in P. fairbanksi is an effect of phenotypic plasticity, it is unclear which factors could cause various morphotype expressions.The specimens from Mongolia, Albania, Canada and Madeira that were measured in our study come from populations cultured in similar laboratory conditions but were started with different counts of founders of various ages, kept in variable densities and with different numbers of generations that had passed culture, so no answer can be proposed at this moment.No molecular markers that correspond with morphological features in tardigrades have been suggested so far.Future studies with higher-resolution markers designed for intrapopulation variation should be performed to determine if any pattern of genetic diversity concordant with morphological variation can be observed.

Parthenogenesis and wide distribution
The phenomenon where parthenogenetic (asexual) lineages occupy a wider geographical range, but sexual populations are restricted to a limited area, is termed 'geographical parthenogenesis' 98 .Guidetti et al. 8 concluded that the difference in the dispersal potential of tardigrades is associated with the two types of reproduction, www.nature.com/scientificreports/i.e., parthenogenetic species show a very wide distribution, inhabiting more continents, while the amphimictic species show a very limited or punctiform distribution.A similar pattern was shown for arthropods where parthenogenesis has been linked with higher dispersal abilities 99 (for example, the freshwater ostracod Eucypris virens 100 (Jurine, 1820) the scorpion species Liocheles australasiae 101 (Fabricus 1775) are parthenogenetic for multiple generations in captivity 102,103 and are widely distributed 104,105 .Similar cases are found in many animals and plants (ref [106][107][108][109] )).However, Baker et al. 99 also suggested that parthenogenesis indicates morphological variation as a result of epigenetic mechanisms.Furthermore, Mioduchowska et al. 87 provided molecular evidence of the presence of the bacterial endosymbiont Wolbachia based on next generation sequencing in tardigrades.
Wolbachia have an effect on the evolution as well as the ecology of their hosts and have been found to cause effects including cytoplasmic incompatibility, feminization, male killing, and induced parthenogenesis 110 It has been noted that at the intraspecific level, even individuals from the same population can undergo morphological changes in their characters to diversify within niches available to the species 111 .Similarly, Kihm et al. 104 proposed epigenetic factors as a main cause for variability in tardigrade Dactylobiotus ovimutans 112 Kihm, Kim, McInnes, Zawierucha, Rho, Kang & Park, 2020 egg morphology, although the population was cultured under controlled laboratory conditions.Despite being rare, it is known that intraspecific variation is caused by external environmental conditions, epigenetics and seasonality 112 .

"Two faces" of cosmopolitism in the Paramacrobiotus
Ecological niche modelling is an important and useful tool that has been used to address issues in many fields of basic and applied ecology 113 .It effectively predicts habitat suitability for rare and poorly studied taxa 114,115 .P. fairbanksi presence is linked to the presence of suitable microhabitats, like moss patches, and their life strategy can make them less likely to be affected by general climatic conditions.However, bioclimatic variables used in the study may be a good predictor of the possibility of the occurrence of suitable microhabitats.We investigated the possible distribution of P. fairbanksi and compared it with other widely distributed species of the genus Paramacrobiotus, i.e., P. gadabouti.Paramacrobiotus fairbanksi already reported from various continents exhibit a cosmopolitan distribution covering different types of environments, whereas P. gadabouti, although also potentially cosmopolitan, has a clear affinity to areas with a Mediterranean climate.Its distribution is poorly known due www.nature.com/scientificreports/ to lack of sampling in many habitats.Such differences clearly show us that even when we consider some of the species to be cosmopolitan, specific patterns of distribution can be completely different.However, we must also stress that the number of known localities for both species relatively low and, in the future, when the number of records of these species will be higher, a distribution pattern may look different.

Conclusions
Paramacrobiotus fairbanksi described originally from Alaska, USA, is now known from almost all zoogeographic realms.The identification of this species is still possible based on morphometric characters alone because the morphometric data are shown to not affect the species identification because of the overlap in measurements of all measured structures.Moreover, the analysis shows low genetic variability among P. fairbanksi populations from various geographical locations, which may in general suggest that interspecies genetic variability in tardigrades is very low too or could be' the effects of Wolbachia infection.The species fits the 'Everything is Everywhere' hypothesis and is an example of a parthenogenetic species with wide distribution.Despite very low genetic variation, some indiscrete morphological variations were observed.Since all the studied populations were cultured and bred in the same laboratory conditions, such variation may have been caused by epigenetic effects, and were not the result of different temperatures, food sources and seasonality.

Figure 1 .
Figure 1.A World map with indicated sample number from Table 1 along with haplotypes of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 3 found in different localities (see also Fig. 8).The world map is from https:// www.wpmap.org/ blank-world-map-with-antar ctica/ blank-world-map-jpg/ and the figure was prepared in Corel Photo-Paint 2021.

Figure 2 .
Figure 2. (A) Differences in the body length (BLm); (B) differences in the buccal tube length (BTLm); (C) differences in the stylet support insertion point (SSIPm); (D) differences in the Microplacoid length (Mim).The studied populations of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 are AL Albania; AQ Antarctic; CA Canada; IT Italy; MD Madeira; MN Mongolia; PL Poland; US USA.Minimum, maximum, median, first quartile and third quartile for each population are presented.All measurements are in micrometres [μm].White boxplots represent cultured population and light grey boxplots represent wild population.

Figure 3 .
Figure 3. (A) Differences in the pt of body length (BLpt); (B) differences in the pt of external width of buccal tube (BTEWpt); (C) differences in the pt of Macroplacoid 1 (M1pt); (D) differences in the pt of Macroplacoid 3 (M3pt).The studied populations of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 are AL Albania, AQ Antarctic, CA Canada, IT Italy, MD Madeira, MN Mongolia, PL Poland, US USA.Minimum, maximum, median, first quartile and third quartile for each population are presented.White boxplots represent cultured population and light grey boxplots represent wild population.

Figure 4 .
Figure 4. (A) Differences in the egg full diameter (EFD); (B) differences in the egg bare diameter (EBD); (C) differences in the egg processes height (PH).The studied populations of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 are AL Albania; AQ Antarctic; CA Canada; MD Madeira; MN Mongolia; PL Poland; US USA.Minimum, maximum, median, first quartile and third quartile for each population are presented.All measurements are in micrometres [μm].White boxplots represent cultured population and light grey boxplots represent wild population.

Figure 5 .
Figure 5. Results of PCA for animal measurements, 1st and 2nd Principal Components.Score scatterplots presented in top-left quadrants; boxplots of single component scores presented in top-right and bottom-left quadrants and loading plot presented in bottom-right.

Figure 6 .
Figure 6.Results of PCA for animal pt indices, 1st and 2nd Principal Components.Score scatterplots presented in top-left quadrants; boxplots of single component scores presented in top-right and bottom-left quadrants and loading plot presented in bottom-right.

Figure 7 .
Figure 7. Results of PCA for egg measurements, 1st and 2nd Principal Components.Score scatterplots presented in top-left quadrants; boxplots of single component scores presented in top-right and bottom-left quadrants and loading plot presented in bottom-right.

Figure 8 .
Figure 8. (A) Median-joining network based on the COI sequences: haplotypes marked as H1-H12 (the number of sequences is given in parentheses), the size of the circles is proportional to the number of sequences, the mutational steps values are indicated along the lines; (B) p-distance value based on the COI barcode sequences; (C) Tajima's D neutrality test.

Table 2 .
Primers with their original references used for sequencing of three molecular markers of Paramacrobiotus fairbanksi.

Table 3 .
GenBank accession numbers of sequences obtained in the present study along with the slide numbers of voucher specimens.

Table 4 . Measurements [in µm] and pt values of selected morphological structures of individuals of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 Albanian population mounted in Hoyer's medium (N-number of specimens/structures measured; RANGE refers to the smallest and the largest structure among all measured specimens; SD-standard deviation, pt-ratio of the length of a given structure
to the length of the buccal tube expressed as a percentage).pt values are in italics.

Table 6 .
Measurements [in µm]and pt values of selected morphological structures of individuals of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 Madeira population mounted in Hoyer's medium (N-number of specimens/structures measured; RANGE refers to the smallest and the largest structure among all measured specimens; SD-standard deviation, pt-ratio of the length of a given structure to the length of the buccal tube expressed as a percentage).pt values are in italics.

Table 7 .
Measurements [in µm]and pt values of selected morphological structures of individuals of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 Mongolian population mounted in Hoyer's medium (N-number of specimens/structures measured; RANGE refers to the smallest and the largest structure among all measured specimens; SD-standard deviation, pt-ratio of the length of a given structure to the length of the buccal tube expressed as a percentage).pt values are in italics.

Table 8 .
Measurements [in µm] of selected morphological structures of eggs of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 Albanian population mounted in Hoyer's medium (N-number of specimens/structures measured, RANGE refers to the smallest and the largest structure among all measured eggs; SD-standard deviation).

Table 9 .
Measurements [in µm] of selected morphological structures of eggs of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 Canadian population mounted in Hoyer's medium (N-number of specimens/structures measured, RANGE refers to the smallest and the largest structure among all measured eggs; SD-standard deviation).

Table 10 .
Measurements [in µm] of selected morphological structures of eggs of Paramacrobiotus fairbanksi Schill, Förster, Dandekar & Wolf 2010 Madeira population mounted in Hoyer's medium (N-number of specimens/structures measured, RANGE refers to the smallest and the largest structure among all measured eggs; SD-standard deviation).

Table 12 .
Results of PERMANOVA and post hoc pairwise PERMANOVA comparisons for the first two principal components (PC1 and PC2) of animal measured values; significant post hoc p-values adjusted with the Benjamini-Hochberg correction.

Table 13 .
Results of PERMANOVA and post hoc pairwise PERMANOVA comparisons for the first two principal components (PC1 and PC2) of animal pt values; significant post hoc p-values adjusted with the Benjamini-Hochberg correction.

Table 14 .
Results of PERMANOVA and post hoc pairwise PERMANOVA comparisons for the first two principal components (PC1 and PC2) of animal pt values; significant post hoc p-values adjusted with the Benjamini-Hochberg correction.