For bacterial amplicons, the V3–V4 variable regions of the 16S rRNA gene were amplified in a one-step PCR as described by Rausch et al. (2019) using the primer pair 341F-806R (5′-CCTACGGGAGG-CAGCAG-30 and 5′-GGACTACHVGGGTWTCTAAT-30) (Kozich et al., 2013). After verification of the presence of PCR products by gel electrophoresis, normalization (Sequal Prep Normalization Plate Kit; Thermo Fisher Scientific, Waltham, United States) and pooling were performed.

For sequencing of archaeal 16S rRNA genes, a nested-PCR approach was chosen using the PCR reactions and cycling conditions described by Pausan et al. (2019). In brief, a first round of PCR was performed using primer pair 344F-1041R (25 cycles), followed by a second round using the universal primer pair 519F-806R (30 cycles). In-between, PCR products were purified using the MinElute PCR Purification kit (Qiagen; Hilden, Germany). Normalization of final PCR products was done using the Sequal Prep Normalization Plate Kit (Thermo Fisher Scientific, Waltham, United States), and pooling was performed. Archaeal and bacterial amplicon sequencing was conducted on a MiSeq platform (MiSeq; Illumina, San Diego, United States) with v3 chemistry. The settings for demultiplexing were 0 mismatches in the barcode sequences.

Amplicon sequence data were processed using the workflows provided by Rühlemann, M. on GitHub at https://github.com/mruehlemann/ikmb_amplicon_processing (Accessed: 17.01.2024). In short, processing included primer clipping using cutadapt 4.1 (Martin, 2011), followed by quality filtering and trimming with the DADA2 package (version 1.42) (Callahan et al., 2016) and the filterAndTrim()-function with parameters truncLen = c (265, 245), maxN = 0, maxEE = c (2, 2), truncQ = 5, rm.phix = TRUE. All sequences were dereplicated, clustered into bacterial or archaeal amplicon sequence variants (ASVs), and merged using DADA2. Each run was processed individually and then merged. After merging of the sequence tables, chimeras were removed using DADA2, and taxonomy was assigned using the Bayesian classifier and the SILVA database v138.1 NR99 (Quast et al., 2013).

For ease of interpretation, bacterial amplicon sequence variants were renamed to BSVs and will be referred to as such in the following, while archaeal amplicon sequence variants remain ASVs. All data analyses were done using Rstudio V2022.07.01 (RStudio Team, 2023) and R 4.3.2 (R Core Team, 2023). The data were transformed and managed using the phyloseq (McMurdie and Holmes, 2013), microViz (Barnett et al., 2021), and ampvis2 packages (Andersen et al., 2018). ASV and BSV tables were decontaminated by first filtering out all extrinsic and unwanted domains, including chloroplasts, mitochondria, eukaryotes, and bacteria or archaea, respectively (Supplementary Figure S1). They were then statistically decontaminated using decontam (Davis et al., 2018) with a prevalence threshold of 0.5 and all negative controls (prep and sequencing controls) as decontamination reference. Identified contaminants (Supplementary Figures S2, S3) were removed from the dataset, as well as samples with less than 500 reads. In the case of the archaeal dataset, 18 remaining ASVs which occurred in negative controls were removed in addition. Data were transformed into relative abundances for heatmaps of the most abundant genera. Principal coordinates analysis was plotted based on Bray–Curtis dissimilarity matrices of all untransformed ASVs and BSVs with an abundance of 0.01% in at least one sample. Whether samples clustered distinctly by host, sample type or environment (aquatic/terrestrial) was assessed using PERMANOVA and pairwise PERMANOVA with 999 permutations and the packages vegan (Oksanen et al., 2020) and pairwiseAdonis (Martinez Arbizu, 2020). Diversity and richness indices were calculated based on ASV/BSV tables rarefied to 10,000 reads, which covered both domains well (Supplementary Figures S4–S6). Differences in diversity means between terrestrial (C. elegans, D. melanogaster, G. melonella, H. sapiens, M. musculus, and T. aestivum) and aquatic (A. aurita, A. aerophoba, H. panicea, H. vulgaris, and M. leidyi) hosts were calculated using Kruskal–Wallis test with a post-hoc Wilcoxon test. The phylogenetic autocorrelation index Moran’s I (Gittleman and Kot, 1990) was calculated using the R package ape (Paradis and Schliep, 2019). This included the Shannon H′ diversity index and Chao1 richness index of the bacterial and archaeal community, respectively, against species divergence times of a molecular time tree available through TimeTree5 (Kumar et al., 2022). Three hosts were not directly available and had to be substituted with their closest available relatives, namely, H. panicea replaced with Superitida, A. aerophoba replaced with Aplysina fulva, and T. aestivum replaced with Hordeum vulgare. Hosts with no detectable archaea community were dropped from the tree for the analysis of Moran’s I and the archaeal diversity and richness indices. The molecular time tree including all hosts was visualized in iTOL (Letunic and Bork, 2021) for use in the graphical abstract. Bacterial and archaeal indicator analysis for the marine hosts A. aerophoba, H. panicea, A. aurita, and M. leidyi was done using the indicspecies package (De Cáceres and Legendre, 2009). Associations between sequence variants with an abundance >0.1% and host organisms were calculated using a multilevel pattern analysis (multipatt) with the “r.g” function and 999 permutations to determine significance. All indicators from the multipatt analysis with a p < 0.05 were plotted in Cytoscape (3.10.1). Association networks of archaea, bacteria, and archaea with bacteria were calculated using SpiecEasi (Kurtz et al., 2015) with Meinshausen-Buhlmann’s neighborhood selection (mb), a lambda of 100, lambda min ratio of 0.001, and 100 replications of the Bounded StARS method. To reduce the sparsity of the data, only sequence variants that occur with at least 3 reads in 10% of the dataset were considered. To estimate robustness of the analysis, 100 random networks (Erdős-Rényi random graphs) were calculated using the same number of nodes and edges or clustering coefficient of each network (Supplementary Figure S7). The random networks and the community-based networks with shared nodes and edges showed the same degree distributions, which resembled a Poisson distribution across 12 degrees (Supplementary Figure S7). Network characteristics were evaluated using the NetCoMi package (Peschel et al., 2021), as well as basic statistics from igraph (Csárdi and Nepusz, 2006), including degree distribution, diameter and radius, betweenness, connectedness (eigenvectors), and the calculation of hubs.

Mice were maintained and handled according to FELASA guidelines and German animal welfare law (Tierschutzgesetz § 11, permit from Veterinäramt Kreis Plön: 1401–144/PLÖ–004697). The study was conducted in accordance with the Declaration of Helsinki, and ethics approval was granted by the ethics committee at Kiel University for human samples (AZ A103/14).

The archaeal dataset started out with a total of 5,871,779 processed reads and 8,775 ASVs, of which 2,805,728 reads and 2,783 ASVs remained after decontamination and removal of Bacteria, Eukarya, mitochondria, and chloroplasts. The bacterial dataset started out with 6,045,133 reads and 20,113 BSVs. After decontamination and removal of unwanted domains, 5,190,629 reads and 19,386 BSVs remained in the dataset. A detailed overview of overall and per host reads and number of sequence variants can be found in Supplementary Table S1. Identified contaminants include human-associated Methanosphaera and Methanobrevibacter, as well as Nitrososphaeraceae in the marine hosts (Supplementary Figure S3A). As archaeal positive controls are based on human Methanosphaera isolates and cannot be distinguished from the real data, this results in strict removal of Methanosphaera ASVs from human samples. Non-human-associated Methanosphaera remain in the dataset. Some of the identified bacterial contaminants include individual BSVs from host-specific symbionts (e.g., Amylibacter or Synechococcus). Where known these were added back into the dataset (Supplementary Figure S3B).

The archaeome of all host organisms is captured by the most abundant 25 archaea (Figure 1A), which belong to the five archaeal phyla Crenarchaeota, Euryarchaeota, Halobacterota, Nanoarchaeota, and Thermoplasmatota. Marine Porifera, Ctenophora, and Cnidaria harbor mainly Cren- and Nanoarchaeota, including Nitrosopumilus with >60% relative abundance. The low-microbial abundance (LMA) sponge H. panicea harbored a high archaeal diversity and richness in comparison with the other marine hosts and particularly the high microbial abundance (HMA) sponge A. aerophoba (n = 9). Comparison of diversity and richness means using the Kruskal–Wallis test on aquatic and terrestrial hosts showed no significant difference between the groups (Shannon H′: χ² = 0.22, value of p = 0.6389; Chao1: χ² = 1.4843, value of p = 0.2231, Figure 2C). Caenorhabditis elegans (n = 5) host uncultured Nitrososphaeraceae, Cand. Nitrocosmicus, as well as Methanobacterium and Methanosphaera, all at 20–40% relative abundance. Humans harbor predominantly Methanobrevibacter and to a lesser extent Methanomethylophilaceae located in the gut (Figure 1A; Supplementary Figure S5). Skin samples harbor unclassified Crenarchaeota, while oral samples showed the genus Methanosaeta (Figure 1A; Supplementary Figure S5). Triticum aestivum as only plant representative hosts Halobacterota, including members of the Rice Cluster II family and Halococcus, but also unclassified Nitrososphaeraceae (Figure 1A). The diversity and richness among terrestrial hosts are highest in the roots of T. aestivum. No single archaea could be detected in the insect hosts G. mellonella, as well as in the mouse subspecies M. musculus musculus and M. m. domesticus. In D. melanogaster and H. vulgaris, five ASVs remained in one sample each belonging to either Cand. Nitrosotalea or Haloferaceae. These, however, were excluded from the dataset as they could not be replicated. We therefore consider them to lack an archaeal community.

The bacterial host-associated community from the same samples is more diverse than the archaeal community (Figure 2) resulting in only one to five abundant taxa (≥1%) per host among the top 25 bacterial genera and a high proportion of remaining taxa (Figure 1B). The sponge A. aerophoba hosts Chloroflexi of the SAR202 clade, as well as the Cand. Synechococcus spongiarum group, while H. panicea shows a member of the genus Amylibacter with a relative abundance above 60%. The number of taxa is reflected in the bacterial diversity, which is higher in A. aerophoba than in H. panicea. M. leidyi’s main colonizers are Pseudoalteromonas and Flavobacterium (Figure 1B). A. aurita harbors the genus Marivita at 30–40% and an unclassified Gammaproteobacterium at 20%. Differences in A. aurita community composition can be attributed to the polyp and medusa life stages (Figure 3; Supplementary Figure S4). While the unclassified Gammaproteobacterium is most abundant in the polyp stage, together with Vibrio and different Rhodobacteraceae, Marivita is the most abundant genus on the A. aurita medusae, alongside Marinobacter. The most abundant genus shared between the life stages is Formosa. The H. vulgaris bacterial community consists of the genera Rhodoferax and Legionella with less than 20% of other genera. Curvibacter which are often present at relative abundances >70% in H. vulgaris (Franzenburg et al., 2013; Taubenheim et al., 2022) were not observed here. This may potentially be linked to the primer sets used in the studies. While here we use the bacterial V3–V4 region, the previously mentioned studies use the V1–V2 region for their amplicon analysis. Further differences may arise due to a change in the extraction kit. As here the ZymoBIOMICS DNA Microprep Kit was used instead of a chloroform/isoamyl alcohol extraction to optimize for the extraction of archaea. Pseudomonas contribute up to 20% of the C. elegans bacterial community. Other genera reach a maximum of 10%, including, for instance, Acinetobacter and Enterobacteriaceae (Figure 1B). C. elegans has the highest bacterial diversity and richness among the studied hosts (Figure 2B). Insect hosts have the lowest bacterial diversity and richness (Figure 2B). While D. melanogaster hosts mainly Acinetobacter and Lactiplantibacillus, G. mellonella harbor Enterobacter, Enterococcus, and Commensalibacter. The most abundant bacterial taxa in the human dataset each represent a particular subsample (Supplementary Figure S5). The oral dataset harbors the bacteria Streptococcus, Veillonella, and Haemophilus, while Staphylococcus and Cutibacterium are particularly abundant on skin. The stool samples contain the most diverse human-associated community including Bacteroides, Alistipes, and Escherichia-Shigella (Figure 1B; Supplementary Figure S5). While mice also host Bacteroides in their guts, their most abundant genera are made up of different Lachnospiraceae and Muribaculaceae. Wheat, in this case the roots, have a heterogenous and diverse bacterial community where the most abundant taxa are present at comparatively low relative abundances (Figures 1B, 2). The genus Massilia represents the most abundant genus at approximately 20% (Figure 1B). Overall, the bacterial diversity (Shannon H′) between terrestrial and aquatic hosts does not differ significantly (χ² = 0.71, p = 0.40). Bacterial richness (Chao1) is significantly higher among the aquatic hosts (χ² = 4.22, p = 0.04, Figure 2B), matching the pattern observed by Rausch et al. (2019).

Principal coordinate analyses on the archaeal and bacterial community composition were done using Bray–Curtis dissimilarity matrices and all taxa with an abundance higher than 0.01% in at least one sample (Figure 3). All host-specific communities differ significantly from one another except for C. elegans and H. sapiens in the archaeal community (Supplementary Table S2). Significant within-host differences among the marine hosts can be observed between the polyp and medusa life stages of A. aurita, both in the bacterial (p < 0.001) and archaeal community (p < 0.01, Supplementary Table S3). The community of H. panicea shows slight seasonal differences in its bacterial community composition (p = 0.02) but not in the archaeal community (p = 0.14). Comparing terrestrial with aquatic hosts, PERMANOVA indicates that communities differ significantly from one another (Archaea: p = 0.001, Bacteria: p = 0.001).

In addition, an approximate time tree of the closest available hosts was used to test whether phylogenetic autocorrelation (Moran’s I) may influence bacterial and archaeal species diversity and richness. Regarding the bacterial community, Moran’s I calculated with the diversity (Shannon H′) and richness (Chao1) of all hosts determined a positive autocorrelation with diversity (p = 0.02, SD = 0.08) but not richness (p = 0.20, SD = 0.07). A closer look at the phylogenetic autocorrelation within terrestrial and aquatic hosts showed that neither aquatic hosts (Shannon H′: p = 0.59, SD = 0.04; Chao1: p = 0.15, SD = 0.03) nor terrestrial hosts (Shannon H′: p = 0.18, SD = 0.21; Chao1: p = 0.20, SD = 0.17) exhibited a significant Moran’s I. Regarding the archaeal community, there was no phylogenetic autocorrelation across hosts for Shannon H′ (p = 0.37, SD = 0.02) or Chao1 (p = 0.66, SD = 0.03). Notably, the four hosts without detectable archaeal community, namely, H. vulgaris, D. melanogaster, G. mellonella, and M. musculus, were excluded from this analysis. Consequently, Moran’s I could not be calculated for the terrestrial hosts as insufficient data were available. The archaeal community among the aquatic hosts also showed no significant autocorrelation with diversity or richness (Shannon H′: p = 0.14, SD = 0.02; Chao1: p = 0.26, SD = 0.01).

Indicator taxa for the different marine hosts were calculated using the point-biserial indicator value. In the case of the bacterial community, Gammaproteobacteria and Alphaproteobacteria are associated with A. aurita and M. leidyi with similar association strength to each host. Overall, more Alphaproteobacteria occur on A. aurita (n = 36) than on M. leidy (n = 16)(Figure 4A). This includes Marivita, which is closely associated with A. aurita and particularly its medusa stage. M. leidyi is further associated with two Firmicutes belonging to the genera Mycoplasma and Spiroplasma. Only three indicator taxa of the Gammaproteobacteria, including Aliivibrio, are shared between M. leidyi and A. aurita, highlighting that their communities are largely distinct. The same applies to the archaeal community, where no indicators were shared among all the marine hosts. While M. leidyi is associated with a member of the Nanoarchaeota, A. aurita is connected solely to Nitrosopumilus ASVs. The number of associated archaeal phyla is higher in H. panicea and A. aerophoba, where we find four and two Thermoplamatota indicators, respectively. Both are preferentially colonized by abundant Nitrosopumilus, however, separate ASVs, including the highly abundant ASV9 (Nitrosopumilus). Bacterial indicators are highly distinct between the two porifera and mirror the trend of HMA versus LMA we find in the bacterial diversity. Archaeal indicators also reflect the higher diversity we find among H. panicea. In numbers we find that H. panicea is associated with 50 BSVs and 22 ASVs, while in A. aerophoba we find 184 BSVs and 5 ASVs to be indicator taxa (Figure 4). Halichondria panicea resembles A. aurita and M. leidyi more closely regarding phylum affiliation and shares three indicators, two Alphaproteobacteria and one Actinobacteriota BSV, with M. leidyi. The association strength of BSVs to the sponge host is, however, stronger than in the ctenophore and cnidarian. Its archaeal indicators are ASVs, belonging to the Nanoarchaeota and the Nitrosopumilus genus. A. aerophoba does not share any indicators with the other hosts and has the most distinct set of phyla as it is associated with 58 Chloroflexi BSVs and its well-known cyanobacterial symbiont Cand. Synechococcus spongiarum (BSV31). The network of archaeal A. aerophoba ASVs is small but shows similar phyla (Crenarchaeota, Nanoarchaeota, and Thermoplasmatota) to the other sponge.

Archaeal association networks by themselves resulted in a random network with a maximum degree distribution of three and corresponding to the random Erdős-Réiny (ER) network. This is mainly due to the small number of ASVs after the prevalence cutoff (n = 32) (Supplementary Figures S7A,D left). The bacterial association network and the joint network show distinct clusters in comparison with the random ER network (Supplementary Figures S6B,D middle). Bacteria alone result in a network of 199 nodes with 588 edges, while the joint network is based on 231 nodes and 758 edges. Comparing the number of edges resulting from 100 iterations of the ER network, we find that these have a substantially higher number of edges (mean = 4,589, SD = 56) and degrees (mean = 39.83, SD = 5.77) when given the same clustering coefficient as the joint bacterial–archaeal network (Supplementary Figure S7E). This indicates that associations in the data driven network are based on true associations rather than random clustering. The bacterial community provides the main structure of the joint network as it provides the majority of BSVs (n = 199). Archaea (n = 32) are nested in close association with each other within the network. There are a total of six main modules in the network, which, except for one (red) are predominantly composed of bacterial associations (Figure 4C). The modules roughly correspond to the different hosts, namely, M. leidyi (purple), H. panicea (red), A. aurita (orange), and A. aerophoba (blue and yellow). Nodes that were identified as hubs or that have the highest centrality scores are of bacterial origin (Figure 4C; Supplementary Table S4) and belong to the Chloroflexi phylum associated with A. aerophoba (Figure 4C). Three of the five most connected nodes (highest degree) belong to the Archaea, namely, Nitrosopumilus and uncultured SCGC AAA286-E23. They all occur within the module resembling the community of H. panicea (Figure 4C middle; Supplementary Table S4), suggesting that taxa within this module interact closely with one another. Nodes with the highest betweenness represent keystone taxa between hosts. While they are present in both sponge hosts, their edges between the A. aerophoba and the H. panicea modules are predominantly negative correlations. This suggests that they may play different roles between modules. Most other modules are linked by positive correlations with each other, indicating that they are mutually exclusive, similarly to the indicator taxa (Figures 4A,B). The H. panicea symbiont Amylibacter has an average number of interaction partners (degree = 6) including Nitrosopumilus (ASV9), while Cand. Synechococcus spongiarum in A. aerophoba is more closely connected within its cluster (degree = 8) and has no connections to archaea.

Here, we present a multi-host study focusing on both the bacterial and the archaeal microbiome. As typically research focuses on only the bacterial community or individual hosts by themselves, studies covering multiple hosts and target bacteria and archaea are less common. Furthermore, by comparing both domains, we go one step further in identifying potential interaction partners within the microbiome. Challenges in studies, such as the one presented here, may arise due to differences in laboratory methods and the taxonomic databases. As described by Rausch et al. (2019), each host requires specialized protocols for DNA extraction, which have been trialed and tested within each working group to optimize extraction of the microbiota. Extractions are, therefore, not standardized across the different models. These may provide different results when breaking down cell membranes within and across different domains (Borrel et al., 2020). In addition, host evolutionary history may have an effect on the observed microbial community (Youngblut et al., 2019). Using Moran’s I (Gittleman and Kot, 1990), our data showed a weak positive phylogenetic autocorrelation only on the bacterial diversity (Shannon H′). This may be linked to the divergence between our terrestrial and aquatic hosts as within-group diversity was not significantly affected. It should be noted, however, that these data represent an approximation as divergence times were not available for all host species. While species divergence times do not seem to greatly affect our observations of the microbiota, the host environment may have a stronger impact on the community. This may particularly affect hosts, which interact strongly with their environment, including the aquatic hosts which are fully immersed in their environment and experience a strong exchange with the water column, or nematodes such as C. elegans, which recruit microbes from their substrates, such as rotting plant matter, known to harbor highly diverse communities (Johnke et al., 2020). While this was beyond the scope of this study, future in-depth analyses of the hosts should include samples of the immediate surroundings. This will allow for better distinction between persistent host–microbiota and taxa which may be brought in by ingestion or flow-through. Furthermore, the depth of taxonomic assignment may vary substantially between hosts, depending on the volume of research that has been put into isolation and characterization of the microbiota. It should also be noted that different taxonomies exist, which are continuously being updated (Parks et al., 2022). Here, we rely on the SILVA taxonomy, which represents a hybrid taxonomy between the most recent version of the GTDB and the NCBI databases. Discrepancies may therefore arise across taxonomic assignments in comparison with the more frequently updated GTDB. For ease of interpretations, we have therefore added differences in taxonomy in the main phyla to the figures. While host–archaea interactions identified here will need to be verified in the laboratory, they provide a starting point for research into the interactions of this domain with a variety of host models from different branches on the phylogenetic tree of life.

Six of the hosts studied here are terrestrial metaorganisms. The archaeome in terrestrial hosts has been of particular interest in a health-related context as they include humans and models for humans. We find that methanoarchaea such as Methanobrevibacter smithii and, to a lesser extent, Methanomethylophilaceae such as Methanomassilicoccus are most abundant in human stool (Figure 1A; Supplementary Figure S5). These results fall in line with previous descriptions of both the human archaeal and bacterial community (Lloyd-Price et al., 2019; Kurilshikov et al., 2021; Thingholm et al., 2021; Chibani et al., 2022). The archaeome in and on humans is generally well described (Koskinen et al., 2017; Pausan et al., 2019). They are present both in healthy adults, as in this study, but have also been associated with different disease phenotypes. Interestingly it has, however, not been possible to identify a single archaeon which follows a pathogenic lifestyle and has been shown to directly cause disease (Borrel et al., 2020; Runge and Rosshart, 2021). As host–microbiome interactions within the human metaorganism are notoriously difficult to study, microbiome research often uses simpler and/or more controllable host systems. Mice are commonly used as a proxy for human microbiota studies (reviewed in Runge and Rosshart, 2021), and mammals typically harbor methanoarchaea as an essential part of their gut microbial community (Youngblut et al., 2021). This, however, does not hold true for young mice where an archaeome seems to be absent. While humans and mice have a similar diversity and richness in their bacterial communities and share similar bacterial genera in their guts (Figures 1, 2), mice did not have an amplifiable archaeal community. This matches previous data analyses which were not able to detect this domain in mice (Beresford-Jones et al., 2022; Kieser et al., 2022). Studies using multi-host analysis to assess the archaeome also do not include data on mice, suggesting that no sufficient data were found (Youngblut et al., 2021; Thomas et al., 2022). Whether archaea remain below the detection limits of amplicons and metagenomes, are acquired at an older age, or are not present naturally is not known. However, a study comparing metagenomes in young rodents suggests that archaea are neither present in laboratory mice nor in wild populations (Bowerman et al., 2021).

Another, terrestrial host model is C. elegans. The archaeal community, which is less diverse than in humans (Figure 2), has thus far not been studied in this metaorganism. Our data show that C. elegans harbor Methanobacterium and Methanosphaera, similarly to humans (Figure 1A). We find that the archaeal community composition differs noticeably between slug and compost-derived C. elegans. This also applies to the bacterial community (Supplementary Figure S5), which can vary substantially based on their natural habitat in compost or slugs (Dirksen et al., 2016; Johnke et al., 2020; Pees et al., 2021). Easier manipulation of the host and its environment allows for analysis of microbiome-mediated evolution in mesocosms (Petersen et al., 2023), adaptation of microbes to the host (Obeng et al., 2023), or mechanisms of microbiome-mediated immune protection (Kissoyan et al., 2022). Considering that many archaea are anaerobes, this finding suggests that C. elegans provides a host environment which is suitable for anaerobic organisms and their potential function, as well as for aerobic organisms, as previously reported (Dirksen et al., 2016; Johnke et al., 2020). Moreover, our findings also identify C. elegans as future model system to study archaea–host interactions in a controlled setting.

While C. elegans looks like a promising candidate to study archaea, this cannot be said for the two arthropod models D. melanogaster and G. mellonella, which lack a replicable archaeal community. None of our samples passed our detection and filtering thresholds, suggesting that an archaeome is absent or present at only a very low abundance The bacterial diversity and richness are among the lowest found across our host models, corresponding to previous findings (Rausch et al., 2019). In the case of arthropods, it has been found that archaea and specifically methanoarchaea are often limited to detritivores (Tokura et al., 2000; Horváthová et al., 2021), indicating that the archaeome is directly linked to the insect food source. As a consistent core microbiome in D. melanogaster is thought to be absent (Wong et al., 2013), it has been recommended to limit microbiome analyses on flies to observations on microbiome dynamics in response to different interventions, rather than focusing on a description of the microbiota (Fink et al., 2013).

The last terrestrial model in this study is the vascular plant model T. aestivum. We find that its roots harbor the highest archaeal diversity across all tested hosts, an observation which may be attributed to its close association with soil. Plants often rely on highly complex microbial interactions within and around their roots and root nodules for the supply of nutrients from soil (Ofek et al., 2012). While beneficial effects of bacteria on growth and stress resistance are frequently tested, archaea are often overlooked from these analyses (Prudence et al., 2019). It therefore remains to be established whether they are similarly attracted by root exudates as beneficial bacteria are.

The marine hosts covered by our study include one ctenophore, a cnidarian, and two porifera with evolutionarily ancient host–microbiota interactions. They are popular targets for microbe–host interaction studies as their relatively simple body structure facilitates research into these interactions (Weiland-Bräuer et al., 2015). Aquatic hosts often receive less attention than medically relevant model hosts, such as M. musculus, resulting in a bacterial and archaeal community that is less well covered by taxonomic databases. Our data reflect this in so far that many archaeal ASVs remain unclassifiable at the genus level and often above (Figure 1; Supplementary Figure S4). Comparing M. leidyi and A. aurita, we find that both archaeal and bacterial communities are distinct from one another with only three shared indicators, a fact which may be linked to their evolutionary history, as ctenophores are a phylum sister to all other animals, whereas sponges and cnidarians branch off at a later time (Schultz et al., 2023). Our results on the bacterial community composition of M. leidyi and A. aurita resemble previous comparative work on the two phyla (Weiland-Bräuer et al., 2015; Jaspers et al., 2020). Depending on their life stage, they, in this case A. aurita, also show a significant intra-host variability in their bacterial and also archaeal community composition (Figure 3; Supplementary Table S3; Supplementary Figure S4). While polyps are colonized by members of the family SCGC AAA286-E23, medusae harbor mainly Nitrosopumilus, which also feature as archaeal indicators (Figure 4). In the bacterial community, we find a high proportion of the Marivita genus, which is closely associated with the medusae, while polyps harbor an abundant unclassified Gammaproteobacterium. This matches previous observations on A. aurita developmental stages (Weiland-Bräuer et al., 2015), which have shown that development is impaired if the polyps lack their natural microbiota at the right time (Jensen et al., 2023). An interesting next step would be to investigate whether the archaeal component of their microbiota also plays an essential role in polyp development.

Porifera hosts include a high microbial abundance (HMA) sponge from the Mediterranean, A. aerophoba, and a low-microbial abundance (LMA) sponge from the Kiel bight, H. panicea. In agreement with HMA-LMA dichotomy in the bacterial community (Hentschel et al., 2006; Moitinho-Silva et al., 2017), we find bacterial diversity to be lower in the LMA sponge (Figure 2B). It is dominated by Amylibacter (Figure 4A), whereas the HMA sponge shows a more complex community composition (Figures 1B, 4A). Amylibacters’ closest relative in the NCBI database is Cand. Halichondribacter symbioticus which has been shown to populate the mesohyl of H. panicea (Knobloch et al., 2020). We find that this symbiont has a medium number of interactions with other taxa (Figure 4), concurring with results by Schmittmann et al. (2022). Among its six interaction partners, we find one persistent colonizer of the HOC36 order (BSV70) as well as ASV9 assigned to the ammonia oxidizing Nitrosopumilus genus (Figures 4A,B). Together with Cand. Halichondribacter, they are involved in the removal and conversion of sponge waste products, such as ammonia (Knobloch et al., 2020). The low number of interactions with bacteria suggests that they are more reliant on exchange with the host than on the remaining microbiota. In opposition to LMA sponges, symbionts in HMA sponges are more diverse and their abundances are more evenly distributed. Commonly, they include the phylum Chloroflexi (Bayer et al., 2018), Cand. Synechococcus spongiarum, or Cand. Poribacteria (Slaby et al., 2017). The association networks show dense clustering around Chloroflexi hub nodes, as well as the Cand. S. spongiarum. This suggests that they interact more strongly with each other, for instance during the cycling of dissolved organic matter from Cand. S. spongiarum to Chloroflexi (Usher, 2008; Bayer et al., 2018).

Remarkably, we find that archaeal diversity is lower in the HMA sponge than in the LMA sponge (Figure 2), which is also reflected in the number of indicator taxa found per sponge species (Figure 4B). Archaeal diversity and richness, thus, show exactly the inverse trend to the common HMA-LMA dichotomy known for associated bacteria (Gloeckner et al., 2014; Moitinho-Silva et al., 2017). This could stem from different control mechanisms for bacteria and archaea in HMA and LMA sponges. If A. aerophoba has more control over its archaeal community composition, it may harbor archaeal members more specific for this host, as is suggested by the negative associations between Aplysina-associated ASVs and the larger clusters linked to the other three marine hosts (Figure 4). However, similarly to bacterial communities, both HMA and LMA sponges are known to harbor sponge-specific archaea that differ from the surrounding seawater (Chaib De Mares et al., 2017; Steinert et al., 2020; Polónia et al., 2021) and that are adapted to the sponge-host lifestyle (Haber et al., 2020; Wang et al., 2022). Therefore, both A. aerophoba and H. panicea are expected to actively control archaea and to be dependent on these symbioses. With a less diverse bacterial microbiome, H. panicea might be more dependent on these functions and metabolites provided by the associated archaeal community and thus select for a more diverse archaeal community. It should be noted, however, that a study comparing A. aerophoba to other LMA sponges found similar archaeal diversity between species (Chaib De Mares et al., 2017). Furthermore, the relative abundances and diversity we discussed here may differ from absolute abundances. Archaeal abundances quantified by qPCR for several sponge species including A. aerophoba correspond to the HMA-LMA dichotomy (Bayer et al., 2014). At this point, we cannot distinguish whether the patterns we observe here are specific to the sponge species or related to their bacterial HMA-LMA status. Future studies could therefore benefit from inclusion of the surrounding seawater, as well as the assessment of microbial ratios through qPCR or metagenome sequencing.

The marine hosts analyzed in this study were found to harbor predominantly archaea from the Nitrosopumilus genus, which are ubiquitous in the marine environment and have previously been observed in association with sponges (Bayer et al., 2008; Chaib De Mares et al., 2017; Steinert et al., 2020). Nitrosopumilus have been shown to be essential players as ammonia oxidizers in the marine nitrogen cycle (Francis et al., 2007; Wright and Lehtovirta-Morley, 2023). We find that they are associated with other nitrogen cyclers, such as Marivita which can reduce oxidized nitrogen compounds during denitrification (Zheng et al., 2019) and are affiliated with A. aurita. A further interaction partner is Nitrosomonas (BSV1013) which are thought to contribute to nitrification in the sponge mesohyl (Hentschel et al., 2006; Bayer et al., 2008). Potential ammonia sources, such as cyanobacteria, including the A. aerophoba symbiont Cand. Synechococcus spongiarum (BSV31), do not associate with archaea in our analysis (Figure 4C). This role may fall to the sponge as it can provide ammonia by itself (Bayer et al., 2008). Nitrosopumilus have also been found to be able to produce oxygen from nitrate under oxygen limited conditions (Kraft et al., 2022). This may provide an additional advantage for survival in the sponge mesohyl, which undergoes diurnal fluctuations in oxygen as reviewed by Zhang et al. (2019).

Other abundant archaea associated with ctenophores, cnidaria, and sponges belong to the Woesearchaeales and individual Thermoplasmata affiliated with the prevalent Marine Group II (Figure 4). Little is currently known about the uncultured order of the Woesearchaeales. Recent research on metagenome assembled genomes shows that they predominantly occur in anoxic environments where they may provide hydrogen to methanoarchaea in a syntrophic relationship (Liu et al., 2018). Woesearchaeota and the SCGC AAA286-E23 have also been detected in the seagrass phyllosphere (Vogel et al., 2021) but are only rarely detected in the marine water column. While they may be adapted to anaerobic niches within the marine hosts, they associate mostly with Nitrosopumilus and aerobic Gammaproteobacteria in our dataset. Marine Group II archaea may supply nutrients to our marine hosts as they are thought to be involved in the initial breakdown of organic matter, peaking after algal blooms (Zhang et al., 2015; Jain and Krishnan, 2021). However, no member of this group has been isolated to date. As taxonomic resolution of Woesearchaeales and Marine Group II remains inexact and undergoes regular reconstruction (Parks et al., 2022), and isolates are still amiss; further research will be needed to discover their taxonomy and role in the metaorganism.

Most hosts in this study show a substantially higher bacterial diversity and richness than in the archaeal community. Grouping the hosts by either terrestrial or aquatic environment, we find that only bacterial richness is slightly significantly higher in the aquatic hosts. This trend in the bacterial community is far less pronounced here in comparison with the cross-model study by Rausch et al. (2019), where terrestrialization was hypothesized to represent a key evolutionary event for the diversification of host-associated microbiota. Differences in our study may arise particularly due to H. panicea, whose high archaeal diversity is responsible for a change in overall group means and differs from our initial expectations. In cases such as this, the low bacterial diversity is linked to the high abundance of a single taxon, which persists in the host and is often passed on directly to its offspring (Schmittmann et al., 2022; Carrier et al., 2023). Variability in symbiont presence/absence and abundance can occur among different host populations but also different sampled body parts. For instance, D. melanogaster communities often show Wolbachia, which are transmitted through their germline (Simhadri et al., 2017). In this case, we used a Wolbachia-free strain explaining the absence of the symbiont; however, the use of extraction protocols in which only the intestine or even only the feces are analyzed can also substantially reduce Wolbachia signals as the majority of Wolbachia endosymbionts are found in the organs of the germline (Fink et al., 2013). Invertebrates, including ctenophores and cnidarians, show species and even population specific differences regarding the presence or absence of highly abundant bacterial taxa (Fraune and Bosch, 2007; Weiland-Bräuer et al., 2015; Jaspers et al., 2020). While we find recurring taxonomic groups in the marine hosts, particularly among the archaea, their ASVs tend to cluster separately from each other. In the case of Nitrosopumilus associated with A. aerophoba or the other marine hosts, association is even negative, suggesting that these may be sponge-specific variants which have co-evolved with their hosts (Theis et al., 2016; Chaib De Mares et al., 2017).

An aspect that may significantly influence host diversity, regardless of its aquatic or terrestrial origin, may be the level of domestication. Hosts addressed in this study generally show very different levels of domestication. They range from freshly caught M. leidyi and sponges to wild-derived but lab-maintained mouse subspecies, and long-standing laboratory populations of H. vulgaris. Hosts, which have a long history of domestication or rearing in the laboratory, such as lab-reared mice, often show a decrease in their bacterial diversity (Rosshart et al., 2019; Bowerman et al., 2021). The same has been observed in plants which often experience a loss of bacterial diversity in association with domestication (Hassani et al., 2020; Özkurt et al., 2020). Here, this trend may also apply to the microbial community of H. vulgaris as it represents a long-standing laboratory model and completely lacks a detectable archaeal community. Further research is, however, needed to understand whether archaea are similarly affected by domestication than host-associated bacteria.