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 locations making it a true cosmopolitan species. The ‘Everything is Everywhere’ (EiE) hypothesis assumes that microscopic organisms have unique features that help them to inhabit many different environments, meaning they can be considered cosmopolitan. In the present work we report four new populations of Pam. 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 Pam. fairbanksi at the genetic and morphological levels. The p-distances between COI haplotypes of all sequenced Pam. 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, and propose epigenetic factor as a main cause rather than temperature, predation risk and food availability. Furthermore, we also discuss differences in the potential distribution of two Paramacrobiotus species.


Introduction
The Phylum Tardigrada currently consists of ca. 1,500 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 microplacoid in pharynx (richtersi group), smooth lunules under all claws, granulation on all legs, and eggs with reticulated conical processes without caps or spines. Pam. 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 20 th century 11,12 suggests that microorganisms and small invertebrates should have a cosmopolitan distribution. Microscopic organisms are often considered cosmopolitan species, as, the presence of specific adaptations allows them to inhabit most environments. These adaptations include a) the possibility of easy passive dispersion (by wind, rivers, sea currents, other animals, etc.), b) the presence of very resistant spore stages (which include cysts, eggs or cryptobiotic individuals) that help to survive extreme conditions, and c) the presence of asexual or parthenogenetic reproduction, allowing for rapid increase in the number of individuals 12,13,14,15,16 . Cosmopolitism was strongly suggested for many tardigrade species in the past, however, the suggestion was later undermined (e.g. ref 17,18,19 ). At present, we have strong and compelling evidence of a wide distribution of some tardigrade taxa, which means that we return to the concept of cosmopolitism of at least some species of tardigrades (e.g. ref 8,9,20,21,22 ) which can support the hypothesis 'Everything is Everywhere' (EiE) for tardigrades. According to Gąsiorek et al 21 , "a species may be termed as cosmopolitan if it was recorded in more than one zoogeographic realm". There are four tardigrade species known from more than one zoogeographic realm, i.e., Echiniscus testudo 23 (Doyère, 1840), Milnesium inceptum 24 Morek, Suzuki, Schill, Georgiev, Yankova, Marley & Michalczyk, 2019, Pam.
Furthermore, two parthenogenetic species from the genus Paramacrobiotus, i.e., Pam. fairbanksi and Pam. gadabouti, are contenders as they show a wide distribution that supports the hypothesis EiE.
In the present paper we compare different populations of Pam. fairbanksi from all known localities of this species in Albania, the Antarctic, Canada, Italy, Mongolia, Poland, Portugal (Madeira) and the USA. We also discuss genetic and morphological differences between them and consider the general distribution of Pam. fairbanksi.

Sample processing
Four moss samples from trees and rocks were collected in 2018 (Mongolia) and 2019 (Albania, Canada and Madeira) (for details, see Table 1, Figure 1). The samples were packed in paper envelopes, dried at room temperature and delivered to the laboratory at the Faculty of Biology, Adam Mickiewicz University in Poznań, Poland. Tardigrades were extracted from the samples and studied following the protocol of Stec et al. 25  Additionally, we used morphometric and genetic data of Pam. fairbanksi populations from the Antarctic, Italy, Spain, the USA and Poland 9 .   Figure 8). The world map is from https://www.wpmap.org/blank-world-map-with-antarctica/blank-world-map-jpg/ and the figure was prepared in Corel Photo-Paint 2021.

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. 26 .
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 27,28 . 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.
All measurements are given in micrometers [μm]. Structures were measured only if their orientation was suitable. Body length was measured from the anterior extremity to the end of the body, excluding the hind legs. Buccal tubes, claws and eggs were measured according to Kaczmarek & Michalczyk 29 . Macroplacoid length sequence is given according to Kaczmarek et al. 30 . The pt ratio is the ratio of the length of a given structure to the length of the buccal tube, expressed as a percentage 31 . The pt values always provided in italics. Morphometric data were handled using the "Parachela" ver. 1.8 template available from the Tardigrada Register 32 .

Genotyping
Before genomic DNA extraction, each specimen of Pam. fairbanksi was identified in vivo using light microscopy (LM). To obtain voucher specimens, DNA extractions were made from individuals using a Chelex® 100 resin (Bio-Rad) extraction method provided by Casquet et al. 35 with modifications described in Stec et al. 25 . We sequenced three molecular markers, which differ in effective mutation rates: two nuclear fragments (18S rRNA and 28S rRNA) and one mitochondrial fragment (COI). All DNA fragments were amplified according to the protocols described in Kaczmarek et al. 9 , with primers listed in Table 2. Alkaline phosphatase FastAP (1 U/μl, Thermo Scientific) and exonuclease I (20 U/μl, Thermo Scientific) were used to clean the PCR products. Sequencing in both directions was carried out using the BigDyeTM terminator cycle sequencing method and ABI Prism 3130xl Genetic Analyzer (Life Technologies). 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  Analyses were performed with 1 000 replicates.

Statistical analysis
We All analyses were carried out using the R software program 49 . 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 50 . 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" 51 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 52,53 . The presence of a structure in the PCA data was tested using a randomization approach on the eigenvalues and statistics according to Björklund 54 Table 1 and Pam. gadabouti from Kayastha et al. 22 . The coordinate list is provided in SM.01 and the R script for ENM in SM.02.

Morphometric comparison of different Pam. fairbanksi populations
No significant differences were shown by the ANOVA test performed on BL between the studied populations (df = 7; F = 7.832; p = 0.902; N = 106; Table 7, Table 9, Table 11, Table 13; Figure 2A). However, significant differences were found on BTL between different populations (df = 7; F = 5.633; p = 0.010; N = 106; Table 7, Table 9, Table 11, Table 13), where the buccal tube of the specimens from Mongolia was significantly longer than in specimens from the Albanian and Canadian populations (p = 0.002 and p = 0.005 respectively; Figure 2B).
The randomization test in PCA demonstrated that only the first two PCs explained greater variation than was anticipated by 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 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 ( Figure 5) explained 90% of the overall variation (83.7% for PC1 and 6.7% for PC2) and for animals' pt indices ( Figure 6) explained 65% of the overall variation (46.3% for PC1 and 18.7% for PC2). PCA of egg measurements ( Figure 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, Table 4, Table 5, Table 6).

Genetic comparisons and phylogeographical analyses of different populations of the Pam. fairbanksi
The COI sequences of Pam. 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 Pam. fairbanksi from Mongolia, and haplotype H4 was found in populations from Canada and Madeira (for details see Table 3 and Figure 8A, B). No stop codons, insertions or deletions were observed. The translation was successfully carried out with the -2 nd reading frame and the invertebrate mitochondrial codon table. The p-distances between COI haplotypes of all sequenced Pam. 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%) ( Figure 8B). In total, twelve haplotypes (H1-H12) of COI gene fragments were identified after comparing all available COI sequences of Pam. 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 Pam. fairbanksi from Albania, Canada, Madeira and Mongolia were 917-1547 bp-long (Table 3) and no nucleotide substitution was found (although a single "N" was identified, i.e. software was unable to identify this base). Compared with the data available in GenBank sequences of Pam. fairbanksi (sequences were alignment 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:

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.,  Figure 9A) to good (yellow areas on the map in Figure 9A) with the most suitable habitats in the northern hemisphere. Pam. gadabouti shows maximal suitability around areas with a Mediterranean climate, although it also has wide distribution ( Figure 9B).  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.

Discussion
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 favor 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 Pam. 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 66 . Various studies combining COI variation and phenotypic plasticity were conducted throughout different invertebrates' phyla, including tardigrades 9,21,66,67,68,69 , 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 67 .
Furthermore, several studies uncovered data incongruence between mitochondrial and nuclear markers, e.g., for earthworms 70 or corals 71 , suggesting that occasionally COI may fail as a barcode marker due to hybridization events. Many studies have already shown that Phenotypic plasticity, in morphology and other aspects of phenotype, such as life history traits, is seen as an advantage for thriving in heterogeneous environments (e.g. ref 74 ), which tardigrades' habitats clearly are. Furthermore, phenotypic plasticity has been widely observed in other invertebrates like corals (e.g. Pseudopterogorgia bipinnata Verrill 1864 75 ), scallops (e.g. Pecten maximus Linnaeus, 1758 76 ), marine invertebrates, gastropods (e.g. Littorina littorea Linnaeus, 1758 62 ), rotifers (Keratella tropica (Apstein, 1907) 77 ) and many more. No concordant genetic variation was observed, but a large and discrete differentiation of morphotypes was present and was always associated with external environmental factors such as temperature, predation risk and food availability 78.79,80,82,82 .

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' 83 . Guidetti et al. 8 concluded that the difference in the dispersal potential of tardigrades is associated with the two types of reproduction, 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 84 (for example, the freshwater ostracod Eucypris virens (Jurine, 1820) 85  Similarly, Kihm et al. 88 proposed epigenetic factors as a main cause for variability in tardigrade Dactylobiotus ovimutans 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 96 . Therefore, it is also likely that the morphological differences that we observed in the present study might be due to epigenetic factors, as the studied populations were cultured under controlled laboratory conditions.

"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 97 . It effectively predicts habitat suitability for rare and poorly studied taxa 98,99 . Pam. 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 Pam. fairbanksi and compared it with other widely distributed species of the genus Paramacrobiotus, i.e., Pam. gadabouti. Paramacrobiotus fairbanksi already reported from various continents exhibit a cosmopolitan distribution covering different types of environments, whereas Pam. gadabouti, although also potentially cosmopolitan, has a clear affinity to areas with a Mediterranean climate. Its distribution is poorly known due 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 is relatively low and, in the future, when the number of records of these species will be higher, a distribution pattern may look different.

Data Availability
The datasets generated and/or analysed during the current study are available in the GenBank repository (all accession numbers listed in Table 2: ON911917-18, ON872386, ON872380-81, ON911919, ON872387, ON872382, ON911920-21, ON872388, ON872383, ON911922-23, ON872389 and ON872384-85). The data of all sequences will be available for public access within a few days.