In summer 2019, lichen samples were collected from apple trees growing in six orchards (Locations 1–6) in Hardin and Franklin counties in north-central Iowa, United States (Supplementary Figures S1, S2). One of these (Location 1) had previously been surveyed for tardigrades (Tibbs-Cortes et al., 2020). One sample of lichen was collected from each tree, and three to five trees were sampled at each location. Moss was also present on the sampled apple trees at Location 2, so moss samples were collected from three of these trees. In June 2020, additional lichen samples were taken from the trees at Location 1 to enable comparison across years. All moss and lichen samples were placed in individual brown paper bags, which were stored in a cool, dry room to allow the samples to dehydrate. From each of the 2020 lichen samples, five subsamples of 0.25 g were placed in sterile 1.5 ml tubes and frozen at −20°C for substrate DNA extraction. Table 1 shows a summary of collected samples.

16S rRNA gene amplicon sequencing studies focusing on low-biomass samples are prone to biases from external contamination during sample processing, DNA extraction, library preparation, and sequencing. Therefore, all subsequent tardigrade isolation and DNA extraction steps were carried out using barrier pipette tips (Axygen) and in a sterile work area dedicated to the project. A Bunsen burner was used to create a sterile field for all tardigrade isolations and DNA extractions.

To extract tardigrades, each moss or lichen sample was soaked in glass-distilled water for a minimum of 4 h. Subsamples of this water were then examined under a dissecting microscope, and tardigrades were extracted with Irwin loops (Miller, 1997; Schram and Davison, 2012). The Irwin loop was disinfected by a flame between each collected tardigrade.

Next, isolated tardigrades were washed by immersion in droplets of PCR-grade water treated with diethyl pyrocarbonate (DEPC). Tardigrades were then transferred to a fresh drop of DEPC-treated water. This washing process was repeated for a total of three washes, and Irwin loops were sterilized between each wash. Three to six replicates of 30 tardigrades each were collected from each substrate sample (identified by replicate codes shown in Table 1) and were then stored in DEPC-treated water at −20°C.

The DNeasy PowerLyzer PowerSoil Kit (Qiagen) was utilized for DNA extraction. Substrate samples were first ground with a sterilized pestle before being transferred to the bead tubes, while tardigrades were directly transferred to bead tubes. Bead tubes were then transferred to a Bead Mill 24 homogenizer (Fisher Scientific). Tardigrades were homogenized using a single 30 s cycle at 5.00 speed, and substrate samples were homogenized using three 30 s cycles with 10 s between cycles at 5.50 speed. Homogenized bead tubes were then centrifuged at 10,000 × g (30 s for tardigrades and 3 min for substrate) before proceeding according to the manufacturer’s instructions; the optional 5 min incubation at 2–8°C was performed during steps 7 and 9. Following elution of DNA with 90 μl of elution buffer, DNA quality and concentration of a 1 μl sample was measured using a NanoDrop spectrometer. Extracted DNA was stored at −20°C.

DNA was loaded onto sterile 96 well plates for library preparation and sequencing. In addition to the tardigrade and substrate samples, three types of controls were included to account for contamination. Six tardigrade processing controls (TPC, equivalent to RIDE sampling blank controls) were created by applying the tardigrade extraction and subsequent DNA extraction protocols to blank samples. Ten DNA processing controls (DPC, equivalent to RIDE DNA extraction blank controls) were created by conducting DNA extraction on 100 μl of DEPC-treated water. Finally, ten wells were loaded with DEPC-treated water to form the library processing controls (LPC, equivalent to RIDE no-template amplification controls). Controls and samples were then submitted for library preparation and 16S rRNA gene amplicon sequencing targeting the V4 region at the Iowa State University DNA facility. Library preparation was conducted following the Earth Microbiome Project 16S Illumina amplicon protocol¹ with the following modifications: (1) a single amplification was conducted for each sample rather than in triplicate, (2) PCR purification was conducted using the QIAquick PCR Purification Kit (Qiagen), and (3) all reactions and purification steps were conducted at half volume using a Mantis liquid handler (Formulamatrix) which was cleaned with isopropanol prior to library preparation. Libraries were loaded onto the MiSeq platform at a concentration of ~4pM, and paired-end sequencing was conducted at 500 cycles.

Following sequencing, three paired end samples representing replicates L6_19_Tr2_li2, L6_19_Tr4_li1, and L5_19_Tr3_li3 (Table 1) were removed from the dataset due to poor quality. Raw reads were processed with mothur version 1.43.0. Sequences were screened to remove reads that contained any ambiguities, were shorter than 252 bases, and had homopolymeric sequences greater than 8 bases. In total, 1,157,089 reads were removed from the raw dataset of 7,805,248 reads. Screened reads were then aligned against the SILVA alignment version 138, and reads which aligned outside the region covered by 95% of the alignment were removed. The SILVA database was also used to remove 145,663 chimeric sequences and to classify remaining sequences. De novo OTU clustering was then conducted at a 99% similarity threshold.

R version 4.0.3 running packages decontam (Davis et al., 2018), phyloseq (McMurdie and Holmes, 2013), DivNet (Willis and Martin, 2020), and corncob (Martin et al., 2020), as well as a more efficient implementation of DivNet known as divnet-rs² running in Rust, were used for subsequent analyses. Using decontam, contaminant OTUs were identified and removed based on their relative prevalence in control vs. true samples (prevalence method, threshold 0.25; Davis et al., 2018). Next, OTUs with fewer than 10 reads in experimental samples were removed. From these data, alpha diversity parameters (Shannon and Simpson) were calculated using DivNet and divnet-rs, which use a log-ratio model and share information across all samples to improve diversity estimates (Willis and Martin, 2020). Relative abundance, differential abundance, and differential variability of taxa were calculated in corncob, which uses a beta-binomial model (Martin et al., 2020). Differences were declared significant when False Discovery Rate (FDR)-corrected p values (Benjamini and Hochberg, 1995) were <0.05. Principal Coordinates Analysis (PCoA) was conducted in phyloseq using the default Bray-Curtis distance.

In the cases where OTUs of interest were not classified by mothur to the genus level, BLAST and RDP Classifier (Wang et al., 2007; Camacho et al., 2009) results were used to provide additional information about taxonomic classification. First, the mothur command “get.oturep” was used to generate a FASTA file containing the representative sequence for each OTU; OTUs with fewer than 10 reads in experimental samples were removed. BLAST analysis was performed using BLAST+ v2.11.0. The NCBI 16S RefSeq collection (representing 22,061 taxa; O’Leary et al., 2016) was downloaded and converted into a BLAST database using the “makeblastdb” command. The “blastn” command was then run against this database using the representative sequence FASTA as the query. Results with the 15 lowest E-values were kept for each OTU. The representative sequence FASTA was also entered into the RDP Classifier web tool version 2.11 using 16S rRNA training set 18,³ and the assignment detail for all OTUs was downloaded.

Raw sequencing files are deposited at the Sequence Read Archive under BioProject accession PRJNA801902 at: https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA801902. Mothur output and code used for analysis is available at: https://github.com/LTibbs/tardigrade_microbiome.

The tardigrade community microbiome differed significantly across locations, as shown by the Simpson and Shannon indices, which differed significantly in most pairwise comparisons of locations (Table 2). Across locations, 13 phyla and 44 genera were both significantly differentially abundant and significantly differentially variable. Sixteen OTUs were significantly differentially abundant only, four OTUs were significantly differentially variable only, and three OTUs were both significantly differentially abundant and variable (Supplementary Table S9). These identified differential taxa included the aforementioned top five phyla (Proteobacteria, Bacteroidota, Actinobacteriota, Firmicutes, and Acidobacteriota) and top three genera (Pseudomonas, Bradyrhizobium, and unclassified Enterobacteriaceae) from the experiment as a whole. Despite these differences, the locations clustered together in the PCoA (Figure 1).

The community microbiome of tardigrades extracted from moss did not differ significantly in alpha diversity from that of tardigrades extracted from lichen as measured by the Shannon and Simpson indices (Table 2). PCoA further demonstrates that the overall microbial community did not differ by substrate type (Figure 1). However, between tardigrades collected from moss and those from lichen, five phyla and 11 genera were both significantly differentially abundant and significantly differentially variable, while three OTUs were significantly differentially abundant only (Supplementary Table S9). These included the common phyla Firmicutes and Bacteroidota; of these, Firmicutes were more abundant in moss-associated and Bacteroidota in lichen-associated tardigrades (Supplementary Table S10, Figure 2).

The microbiota of tardigrades was significantly less diverse than that of their lichen substrate, as measured by both Shannon and Simpson indices (Table 2); the tardigrade and substrate samples also formed distinct clusters as shown by PCoA (Figure 1). Between tardigrades and their substrate, 17 phyla, 181 genera, and 101 OTUs were significantly differentially abundant and variable, while 308 OTUs were significantly differentially abundant only and 124 OTUs were significantly differentially variable only (Supplementary Table S9). These differential taxa included four of the top five phyla (all except Proteobacteria) and all three of the three most abundant genera from the experiment as a whole. Remarkably, the relative abundance of Firmicutes was nearly a 1,000 times higher in the tardigrades (20.5%) than in their substrate (0.021%; Supplementary Table S10, Figure 2).

From 2019 to 2020, the tardigrade community microbiome increased in diversity as measured by the Simpson index, though no significant difference was found between the Shannon indices (Table 2). The two years also formed mostly distinct clusters in the PCoA (Figure 1). Between the two years, 44 genera and one OTU were significantly differentially variable and abundant, while two phyla and 26 OTUs were significantly differentially abundant only (Supplementary Table S9). These differential taxa included two of the five most common phyla (Proteobacteria and Actinobacteriota) and two of the three most common genera (Pseudomonas and unclassified Enterobacteriaceae). The unclassified Enterobacteriaceae had a particularly large change in relative abundance, decreasing more than 160-fold from 10.1% in 2019 to 0.062% in 2020 (Supplementary Table S11; Figure 3).

This study examined the microbiota of the full tardigrade community from a particular substrate sample, in contrast to previous surveys that studied isolated species, often from laboratory cultures (Vecchi et al., 2018; Kaczmarek et al., 2020; Mioduchowska et al., 2021; Zawierucha et al., 2022). While the focus of the current study was the bacteria associated with the full tardigrade community in a given sample, it is of course also desirable to identify species-specific tardigrade microbiota, as Vecchi et al. (2018) found that tardigrade-associated bacteria varied among tardigrade species, and we encourage further research on bacteria associated with individual tardigrade species. However, while it is occasionally possible in samples containing only a few tardigrade species to extract and immediately identify members of each for study (Vecchi et al., 2018), many environmental samples contain numerous species, including cryptic species that may be difficult or impossible to distinguish without molecular, life cycle, or other data (Cesari et al., 2013; Guidetti et al., 2016). For example, all three of the tardigrade genera observed in the current survey (Milnesium, Ramazzottius, and Paramacrobiotus) would require additional morphometric or egg observations to identify members to species level (Binda and Pilato, 1987; Kinchin, 1996; Guidetti et al., 2009; Michalczyk et al., 2012). While it would have been possible to culture collected tardigrades in the laboratory to collect eggs and generate enough individuals to both mount for morphometric species identification and to process for bacterial DNA, Vecchi et al. (2018) found that culturing significantly affects the tardigrade microbiome. Therefore, the current study has the unique advantage of better reflecting the tardigrade community microbiome in its natural state compared to studies that focus on cultured tardigrades.

While tardigrade species were not identified in this study, a December 2015 collection effort at Location 1 provides information on tardigrade diversity in the area. From lichen growing on some of the same apple trees used in the current study, the previous study identified Milnesium cf. barbadosense Meyer and Hinton 2012; Mil. burgessi Schlabach, Donaldson, Hobelman, Miller, and Lowman, 2018; Mil. swansoni Young, Chappell, Miller, and Lowman, 2016; and P. tonollii (Ramazzotti, 1956), as well as members of Milnesium Doyère, 1840; Ramazzottius Binda and Pilato, 1986; Paramacrobiotus Guidetti, Schill, Bertolani, Dandekar and Wolf, 2009; and Macrobiotidae Thulin, 1928 not identifiable to species (Tibbs-Cortes et al., 2020). While tardigrade communities are dynamic across both time (Schuster and Greven, 2007, 2013) and space (Meyer, 2006, 2008), the dominant species present in tardigrade communities in a given area can remain remarkably stable across years (Nelson and McGlothlin, 1996; Schuster and Greven, 2007). This suggests that species information from the 2015 survey may be relevant to the current study.

All tardigrade microbiome surveys, including the current study, employed laboratory technique to reduce contamination by washing the tardigrades in sterile water before DNA extraction (Vecchi et al., 2018; Kaczmarek et al., 2020; Mioduchowska et al., 2021; Zawierucha et al., 2022). Working in a sterile environment further decreases contamination; therefore, Zawierucha et al. (2022) extracted tardigrades from substrate in a sterile environment (laminar flow chamber), and we worked in a sterile field created by a Bunsen burner throughout the experiment. Three previous studies included one type each of negative controls recommended by the RIDE standards for low biomass studies (Eisenhofer et al., 2019; DNA extraction blank in (Mioduchowska et al., 2021), sampling blank in (Kaczmarek et al., 2020), and no-template amplification control in (Zawierucha et al., 2022)). Kaczmarek et al. (2020) and Mioduchowska et al. (2021) did not sequence the negative controls; instead, they performed PCR amplification of these controls and determined that no contamination was present because no bands were visible. However, samples without visible bands from PCR can generate sequencing reads (Davis et al., 2018), and this method would not detect contaminants introduced during library preparation or sequencing steps. Zawierucha et al. (2022) removed all OTUs present in the no-template amplification control from analysis. However, low levels of true sequences, especially from high-abundance OTUs, are often present in negative controls due to cross-contamination of samples; these biologically important OTUs would therefore be removed from the analysis (Davis et al., 2018; Karstens et al., 2019). In our study, we included and sequenced all three RIDE-recommended types of negative controls.

No previous survey of the tardigrade microbiota has employed model-based in silico contaminant identification and removal, which we accomplished using the decontam package. Decontam removed 986 OTUs as contaminants, including five that would otherwise have been in the top 10 OTUs in the study by read count. Of course, further improvements are always possible. Contaminants are expected to differ in prevalence among negative control types depending on their point of introduction, but current in silico contamination removal methods treat all negative controls identically (Davis et al., 2018). Future development of a method that leverages the unique information provided by each type of negative control would therefore be desirable. However, by working under a flame, including and sequencing all recommended types of negative controls, and leveraging in silico contaminant identification and removal, we have produced what we expect to be the tardigrade microbiome survey least affected by contamination to date.

We investigated the tardigrade community microbiome across four contrasts. First, we determined that the tardigrade community microbiome varied significantly in structure across locations (Table 2). Vecchi et al. (2018) found that an average of 15.4% of the microbial OTUs in a tardigrade collected from moss or lichen originate from its substrate. Therefore, known impacts of geographical location on microbial communities (Baldrian, 2017; Coller et al., 2019) could have resulted in different microbial communities present in each location to inoculate the tardigrades. Additionally, as the tardigrade microbiota is species-specific (Vecchi et al., 2018), the differences in microbial communities observed across locations may reflect spatial variation in tardigrade communities’ species composition (Meyer, 2008). Of course, these explanations are not mutually exclusive and could both play a role in shaping distinct tardigrade community microbiomes across locations. Vecchi et al. surveyed tardigrades of the same species collected from different locations, but they did not test for differences in the microbiome across locations. However, they did identify an OTU in the genus Luteolibacter that was significantly associated with Ram. oberhaeuseri collected from a location at 34 meters above sea level but not from another location 797 meters above sea level (Vecchi et al., 2018), suggesting that future work may detect differences in the microbiota of the same tardigrade species across locations.

In contrast two, the community microbiome of tardigrades collected from lichen was compared with that of tardigrades collected from moss on the same trees. While a few taxa were significantly differentially abundant and variable across substrates, the two substrates did not differ in alpha diversity and were not separated by PCoA (Figure 1). This similarity in the tardigrade community microbiome was initially surprising, as previous literature has demonstrated significant differences between the microbiota of moss and lichen even on the same tree (Aschenbrenner et al., 2017). However, this similarity across substrates could be due to the presence of similar tardigrade species, as previous studies have failed to demonstrate significant differences between tardigrade communities found in moss and lichen (Young and Clifton, 2015; Nelson et al., 2020). In future surveys, it would be interesting to compare the microbiota of tardigrades from additional substrate types (e.g., soil) and to determine if this similarity in tardigrade microbiome across substrates persists at the species as well as the community level.

Results of contrast three demonstrate that the tardigrade community microbiome is distinct from and significantly less diverse than that of its lichen substrate (Table 2, Figure 1). This result agrees with previous studies that found relatively higher diversity in substrates than in their resident Ecdysozoans, including wild tardigrades collected from moss and lichen (Vecchi et al., 2018), cultured tardigrades (Kaczmarek et al., 2020), and the nematodes Meloidogyne hapla (Adam et al., 2014) and Caenorhabditis elegans (Johnke et al., 2020). This study is the first to demonstrate this trend at the tardigrade community rather than species level. Vecchi et al. (2018) suggested that the lower microbial diversity in tardigrades with respect to their substrates may be due to the small size of tardigrades limiting the biomass and therefore the diversity of their microbiome (small host hypothesis) and/or to selectiveness of tardigrades inhibiting growth of some bacterial species and promoting growth of others (selective host hypothesis). Supporting the selective host hypothesis, earlier work found that tardigrades could be successfully inoculated with some bacteria (Xanthomonas) but not others (Serratia; Krantz et al., 1999). It is possible that some bacteria have co-evolved with tardigrades, becoming permanent residents of the gastrointestinal tract or cuticle, a hypothesis that has been suggested for the Ecdysozoan C. elegans (Zhang et al., 2017). The life cycle of tardigrades poses a unique selective pressure on any permanent residents of the microbiota, as these organisms would also have to survive within the tardigrade during cryptobiosis. This would be especially true for the obligate endosymbiotic taxa Rickettsiales and Polynucleobacter previously observed in tardigrades (Vecchi et al., 2018; Guidetti et al., 2020; Kaczmarek et al., 2020; Mioduchowska et al., 2021), as well as for the Rickettsia identified in the current study (see below).

Contrast four determined that the tardigrade community microbiome is temporally dynamic, changing significantly on the same trees from 2019 to 2020 (Table 2; Figure 1). Again, this may be due to changes in habitat microbiome, as microbiota of other substrates (e.g., soil and litter) are known to vary across years due to changing environmental factors such as nutrient availability (Martinović et al., 2021). This variation may also be due to temporal changes in the tardigrade community composition; although tardigrade species present may remain consistent in a location over years, their relative abundances shift in part due to changes in rainfall, humidity, and temperature (Schuster and Greven, 2007). This temporal variability raises important implications for future studies of the tardigrade community microbiome. For example, the relative abundance of putative phytopathogens differed significantly across years (Supplementary Table S9). Future work could identify temporal variables affecting the ability of tardigrades to act as potential reservoirs of phytopathogens and other bacteria. We also encourage further studies of the tardigrade microbiome to account for temporal changes and to investigate this variation with additional time points to increase resolution.

In this study, the five most abundant phyla were Proteobacteria, Firmicutes, Bacteroidota, Actinobacteria, and Acidobacteria (Supplementary Table S10). All of these except Acidobacteria were previously reported as highly abundant in at least two of the three previous tardigrade microbiome surveys that presented results at a phylum level, with Proteobacteria identified as the most abundant phylum in all cases (Vecchi et al., 2018; Kaczmarek et al., 2020; Mioduchowska et al., 2021). Combined, the tardigrades in these studies represent a diverse set of species, including wild and laboratory-reared specimens isolated from multiple continents, suggesting that the predominance of these phyla is broadly characteristic of the microbiome of Tardigrada, regardless of species or location. These phyla, especially Proteobacteria, are also dominant in the microbiomes of other Ecdysozoans, including soil nematodes (Adam et al., 2014; Dirksen et al., 2016; Elhady et al., 2017), marine nematodes (Arcos et al., 2021), and insects (Colman et al., 2012; Engel and Moran, 2013). The tardigrade microbiota therefore appears similar to that of other Ecdysozoans at the phylum level.

A number of OTUs significantly more abundant in tardigrades than in their substrate in this study belong to taxa previously identified in the tardigrade microbiome. These include members of Enhydrobacter (Vecchi et al., 2018), Enterobacteriaceae (Kaczmarek et al., 2020; Mioduchowska et al., 2021), and Acinetobacter (Vecchi et al., 2018; Mioduchowska et al., 2021; Supplementary Table S9). In this study, OTU 22 was classified as Enhydrobacter, and its abundance in the tardigrade population varied over time, increasing significantly from 2019 (0.21%) to 2020 (2.1%; Supplementary Tables S8, S9). It is possible that Enhydrobacter is common to Ecdysozoan microbiomes, as it is also an abundant taxon in the gut contents of larval wood wasps (Li et al., 2021) and nematodes (Adam et al., 2014). Members of Enterobacteriaceae included OTUs 2 and 20. OTU 20 was further identified as a member of the Escherichia/Shigella complex, but OTU 2 could not be classified to the genus level (Supplementary Tables S2, S3). OTU 2 showed significant temporal variation, decreasing in relative abundance from 10.0 to 0.062% from 2019 to 2020 (Supplementary Tables S8, S9). Enterobacteriaceae is also highly represented in the gut microbiota of insects (Hernández-García et al., 2017; Moro et al., 2021) and nematodes (Zimmermann et al., 2020; Zhou et al., 2022). This suggests that Enterobacteriaceae may be residents of the tardigrade digestive tract. Finally, OTU 16 was a member of Acinetobacter that increased significantly in abundance from 2019 (0.000052%) to 2020 (2.3%) and was one of the three OTUs significantly differentially abundant across substrate type (Supplementary Tables S8, S9). Acinetobacter is associated with the cuticle of the nematodes M. hapla, M. incognita, and Pratylenchus penetrans (Adam et al., 2014; Elhady et al., 2017), suggesting that it may also be associated with the cuticle of tardigrades.

However, many of the tardigrade-associated taxa observed in this study have not been previously reported in the tardigrade microbiome. In fact, the most abundant OTU across all samples in this study (OTU 1) was a member of the genus Bradyrhizobium, which was not previously reported from the tardigrade microbiome. This OTU was also spatially dynamic, differing significantly in abundance across locations (Supplementary Table S9). Bradyrhizobium has been previously observed in the microbiota of plant pathogenic nematodes (Adam et al., 2014; Eberlein et al., 2016) and leaf hoppers (Horgan et al., 2019). This genus has also been found in the lichen microbiome (Bates et al., 2011; Erlacher et al., 2015; Graham et al., 2018), perhaps indicating that tardigrades acquire this bacterium from their habitat. Another tardigrade-associated genus, Micrococcus, was differentially abundant across both locations and years (Supplementary Table S9). This genus has been reported from the cuticles of soil nematodes (Adam et al., 2014) as well as fish parasitic nematodes (Arcos et al., 2021), suggesting that Micrococcus may be associated with the tardigrade cuticle. Another notable tardigrade-associated genus in this study was Nakamurella, represented primarily by OTU 33, which showed differential abundance across locations (Supplementary Table S9). Nakamurella intestinalis has been isolated from the feces of another Ecdysozoan, the katydid Pseudorynchus japonicus (Kim et al., 2017). Nakamurella endophytica and N. flava were identified as endophytes of mangroves and mint, respectively (Tuo et al., 2016; Yan et al., 2020), and N. albus and N. leprariae were originally discovered in lichens (Jiang et al., 2020; An et al., 2021). This suggests that tardigrades could obtain endophytic or lichen-dwelling Nakamurella from their habitat.

Our survey corroborates previous observations of putative endosymbionts of the obligate intracellular order Rickettsiales associated with tardigrades. Three of the four previous surveys of the tardigrade microbiome have detected OTUs of this order (Vecchi et al., 2018; Guidetti et al., 2020; Kaczmarek et al., 2020; Mioduchowska et al., 2021); in addition to unclassified Rickettsiales, these OTUs included members of Wolbachia (Mioduchowska et al., 2021), Anaplasmataceae, and Ca. Tenuibacteraceae (Guidetti et al., 2020). Kaczmarek et al. also detected the obligate intracellular genus Polynucleobacter (2020). In the current survey, we identified two Rickettsiales OTUs. Of these, one was classified by mothur as Wolbachia (OTU 3606), and the other was further classified as Rickettsia (OTU 180) by BLAST and RDP analysis (Supplementary Tables S2, S3). The relative abundance of OTU 180 was significantly higher in tardigrades (0.88%) than in their lichen substrate (0.0012%; Supplementary Table S7) as well as significantly higher in 2020 (1.0%) than in 2019 (0.00026%; Supplementary Tables S8, S9). OTU 3606 was numerically more abundant in tardigrades (0.030%) than substrate (5.1 × 10⁻¹²%), though this difference was not statistically significant (Supplementary Table S7). Taken together, the intracellular nature of Rickettsiales and the higher abundance in tardigrades suggests that OTUs 180 and 3606 are endosymbionts of tardigrades.

The presence of endosymbionts may have implications for tardigrade reproduction and evolution, as members of Rickettsia and Wolbachia are known to manipulate host reproduction in other Ecdysozoans. Wolbachia is well-known for causing parthenogenesis in nematodes and arthropods, as well as feminization of males, cytoplasmic incompatibility, and male-killing (Werren et al., 2008; Kraaijeveld et al., 2011; Correa and Ballard, 2016; Kajtoch and Kotásková, 2018). Similarly, Rickettsia can induce parthenogenesis (Hagimori et al., 2006; Giorgini et al., 2010) and male-killing (Lawson et al., 2001) in arthropods. Parthenogenesis is common in tardigrades (Bertolani, 2001; Guil et al., 2022). Further investigation is necessary to determine if this is due to reproductive manipulators such as Rickettsia and Wolbachia. Future analysis could follow the example of Guidetti et al. (2020) by incorporating FISH to confirm the presence of these and other endosymbionts within tardigrade tissues.

Our analysis also aimed to determine whether wild tardigrades living in apple orchards harbor phytopathogenic bacteria, and in fact, the second most abundant genus overall found in this survey was Pseudomonas, which contains more than 20 known plant pathogens (Höfte and De Vos, 2006). Pseudomonas was significantly associated with tardigrades (relative abundance in tardigrades and substrate of 2.7% and 0.051%, respectively; Supplementary Table S11) and was spatially and temporally dynamic in the tardigrade community microbiome, as relative abundance of Pseudomonas decreased significantly from 2019 (19.6%) to 2020 (3.0%) and differed significantly across locations (Supplementary Tables S9, S11). Pseudomonas was also detected in all four of the previous surveys of the tardigrade microbiome (Vecchi et al., 2018; Kaczmarek et al., 2020; Mioduchowska et al., 2021; Zawierucha et al., 2022), and Vecchi et al. (2018) identified it as part of the core tardigrade microbiome. Pseudomonas is also present in the microbiota of soil nematodes (Adam et al., 2014; Dirksen et al., 2016; Zimmermann et al., 2020) and insects (Hernández-García et al., 2017; Horgan et al., 2019; Xue et al., 2021). Notably, other Ecdysozoans (insects) act as vectors of P. syringae (Stavrinides et al., 2009; Donati et al., 2017), which is one of the most agriculturally damaging Pseudomonas species (Höfte and De Vos, 2006; Xin et al., 2018). However, Pseudomonas is very diverse, containing many non-pathogenic species (Silby et al., 2011; Passera et al., 2019). In fact, some Pseudomonas isolates from wild C. elegans confer resistance to fungal pathogens in their hosts (Dirksen et al., 2016), raising the possibility that Pseudomonas could be similarly beneficial to tardigrades.

Two additional putative phytopathogens were significantly more abundant in tardigrades than their substrate. The first, OTU 261, was identified by mothur as a member of Ralstonia, a genus containing the phytopathogenic R. solanacearum complex. In addition to being found at significantly higher abundance in tardigrades (0.018%) than their substrate (0.00017%; Supplementary Table S7), OTU 261 was temporally dynamic, decreasing significantly from 2019 (0.16%) to 2020 (0.00068%; Supplementary Tables S8, S9). Ralstonia has been previously observed in the tardigrade P. fairbanksi (Mioduchowska et al., 2021) and in nematodes (Eberlein et al., 2016; Elhady et al., 2017). The R. solanacearum complex causes major yield losses in food crops including tomatoes, bananas, and potatoes (Yuliar et al., 2015; Paudel et al., 2020). Two notable members of this complex are spread by insect vectors; the cercopoids Hindola fulva and H. strata act as vectors of R. syzygii, while the Blood Disease Bacterium is spread nonspecifically by pollinators (Eden-Green et al., 1992; Remenant et al., 2011).

The second, OTU 208, was classified by BLAST and RDP analysis to the Erwinia/Pantoea cluster (Supplementary Tables S2, S3), which includes a number of economically important phytopathogens (Kido et al., 2008; Zhang and Qiu, 2015; Dutkiewicz et al., 2016; Shapiro et al., 2016). Erwinia amylovora is of particular note as it causes fire blight in apple trees (Aćimović et al., 2015). This OTU had a significantly higher relative abundance of 0.046% in tardigrades compared to 0.0044% in their substrate (Supplementary Tables S7, S9). While neither Erwinia nor Pantoea have previously been identified in tardigrades, Erwinia has been found in arthropods (Xue et al., 2021) and nematodes (Eberlein et al., 2016). Additionally, multiple phytopathogens in Pantoea and Erwinia are transmitted by insect vectors (Basset et al., 2000; Sasu et al., 2010; Ordax et al., 2015; Walterson and Stavrinides, 2015; Dutkiewicz et al., 2016). However, it is also possible that OTU 208 represents a symbiont in tardigrades, as Erwinia also includes the olive fly obligate gut symbiont Candidatus Erwinia dacicola (Blow et al., 2020).

Additional putative plant pathogens were observed at lower abundances in the tardigrade community microbiome and were not significantly more abundant in tardigrades than their substrate. These include another Ralstonia (OTU 1556) and OTU 1620, which was classified as Pectobacterium by BLAST and RDP (Supplementary Tables S2, S3). Members of Pectobacterium cause soft rot diseases in economically important plants, and some strains are capable of infecting multiple plant species (Ma et al., 2007; Li et al., 2020). Additionally, prompted by previous observation of the tardigrade M. hufelandi acting as a vector of the plant pathogen X. campestris (Benoit et al., 2000), we searched the tardigrade community microbiome for members of Xanthomonas. OTUs 10,409 and 12,281 were classified as Xanthomonas (OTU 10,409 by BLAST and RDP analysis), but were both at extremely low abundance (Supplementary Table S4).

In summary, we observed the presence of multiple putative plant pathogens in the community microbiome of tardigrades isolated from apple orchards. Tardigrades could act as vectors or reservoirs of these putative pathogens, a possibility raised by the previous observation of M. hufelandi as a vector of X. campestris (Benoit et al., 2000). However, a major limitation of this study is the use of only 16S rRNA amplicon sequencing. Because multiple marker genes are required to distinguish among species within Pseudomonas, Ralstonia, Erwinia, and Pantoea (Gomila et al., 2015; Zhang and Qiu, 2015; Palmer et al., 2017; Paudel et al., 2020; Saati-Santamaría et al., 2021), we were unable to identify OTUs in our study to species level. Therefore, we are unable to determine whether the identified OTUs in plant pathogenic genera are themselves phytopathogens. We encourage future analyses of tardigrade-associated bacteria in these groups through techniques such as metagenome sequencing and multilocus sequence typing to clarify this point.