Soils and soil biodiversity are becoming increasingly valued and protected at the policy level1. Soil invertebrates are major components of soil biodiversity, and their activity is important for almost all soil ecosystem services2. For example, soil invertebrates are responsible for up to 50% of the litter decomposition3. They contribute to functional services crucial to humans, such as nutrient cycling, water storage and support above-ground food production through the integration of nutrients in food webs4–6. Furthermore, soil invertebrates play major roles in regulating microbial activity along the plant-soil continuum75. Consistent with their importance in soil ecosystems, they are actively promoted in soil biodiversity conservation frameworks8.
However, soil invertebrates are inherently difficult to study morphologically due to their incredible diversity, huge abundances, and small body size with microscopic morphological details. Though generally tiny, they show a ~ 100-fold variation in body weight, which ranges from nanograms to grams9. There are potentially hundreds of thousands of undescribed species10. Moreover, taxonomic expertise is declining11 and this is particularly problematic for groups where experts have always been rare.
DNA- and RNA-based methods are long promoted to support traditional taxonomy and ecological studies in difficult organism groups. Shotgun metagenomics randomly sequences DNA fragments from a sample, instead of relying on PCR-amplified taxonomic marker genes. Metagenomics is an increasingly feasible approach to record the presence of higher eukaryotes in a diverse range of samples12–14. Since metagenomics can utilise all genomic information for taxonomic identification, it has improved sensitivity and specificity compared to metabarcoding15, and it promises superior quantification of species’ biomass16. Metatranscriptomics in turn records genes which are actively transcribed into RNA and thus drive ongoing biological processes17, informing about the metabolic activity of soil community members18, and functional changes in these communities17.
Comprehensive genome collections are the backbone for metagenomics and metatranscriptomics. If genome databases are available, shotgun metagenomics and metatranscriptomics have shown to provide unprecedented insights17,19, e.g., into vegetation change over glacial cycles15, historic population genomic processes20,21, and kingdom-spanning processes of ecosystem functioning22. Large genome sequencing initiatives like the Earth Biogenome Project23 will provide this data ultimately, but progress so far mainly focused on large, prominent organisms, such as mammals24, birds25, insects26 and plants15. In addition to serving taxonomic identification, broad (many distinct groups) and dense (many species from a group) sequencing of genomes additionally allows identifying common patterns of gene evolution and test the taxonomic generality of hypotheses on genome evolution.
Here, we have generated a large genomic resource to support insights into the structure, activity and functioning of soil invertebrate communities. We sequenced and assembled the genomes of 232 species, representing 14 common soil invertebrate groups (hereafter referred to as “groups”) encompassing 94 families, most of which were lacking whole-genome data thus far (Fig. 1, Table 1), including Collembola (n = 87 species), Oribatida (n = 62), two classes of Myriapoda (n = 23 Diplopoda; n = 19 Chilopoda) and soil Nematoda (n = 18). Genome completeness estimated with benchmarking universal single copy orthologs (BUSCO)27 was 59.78% on average, with an average contig N50 of 28,375 bases (Supp. Table 1).
Table 1
Overview of 232 soil invertebrate genome assemblies. For each taxonomic group we also list the number of species with publicly available genome assemblies (as of June 2022).
Phylum | Taxon group [rank] | Common name | n known species (soil or terrrestrial) | n species (published genomes) | n species (genomes contributed here) |
Annelida | Lumbricidae [family] | Earthworms | 7000 | 2 | 1 |
Annelida | Enchytraeidae [family] | Potworms | 700 | 1 | 7 |
Nematoda | Nematoda [phylum] | Nematodes | 25000 | 73 | 18 |
Tardigrada | Tardigrada [phylum] | Tardigrades | 1150 | 4 | 2 |
Arthropoda | Gamasina [infraorder] | Predatory mites | 40000 | 1 | 4 |
Arthropoda | Astigmata [suborder] | Mites [not soil] | 7 | 1 |
Arthropoda | Oribatida [suborder] | Box mites | 7 | 62 |
Arthropoda | Chilopoda [class] | Centipedes | 3000 | 2 | 19 |
Arthropoda | Diplopoda [class] | Millipedes | 12000 | 3 | 23 |
Arthropoda | Symphyla [class] | Symphylans | 200 | 0 | 2 |
Arthropoda | Pauropoda [class] | Pauropods | 800 | 0 | 1 |
Arthropoda | Isopoda [order] | Pill bugs | 3637 | 5 | 3 |
Arthropoda | Diplura [order] | Diplurans | 1000 | 2 | 2 |
Arthropoda | Collembola [class] | Springtails | 8500 | 35 | 87 |
To demonstrate the relevance of this new genomic resource, we first used the 232 genomes to re-analyse the taxonomic assignment of metatranscriptomic sequences generated from a 2-year sampling of environmental RNA along an elevation gradient28. Such assignments were previously not feasible due to a general lack of soil invertebrate genome data. We found that about 1.21% (854,409 reads) of the classified metatranscriptomic reads (70,728,525) could be assigned to soil invertebrates, in comparison to bacteria (85.5%. 60,474,505), fungi (13.29%, 9,399,611), and plants (0.19%, 134,852). Previous metatranscriptomic studies reported a similar microbial eukaryote to bacteria ratios29,30. We recorded 11 soil invertebrate classes (Fig. 2a), of which the most abundant were nematodes of the class Chromadorea as well as clitellates (comprising both earthworms and enchytraeids). There was a marked dip in soil invertebrate richness along the elevation gradient (Fig. 2b). This is in contrast with findings of hump-shaped elevation - richness relationships in soil invertebrates31. The pattern observed by us might be driven by distinct vegetation covers, although this needs better sampling resolution to confirm. The abundance of soil invertebrates showed statistically significant differences across study years (Fig. 2c, analysis of deviance for multivariate generalised linear model fits, df = 5, dev = 929.87, p = 0.05). This might reflect differences in the amount of precipitation between the study years, as the activity of soil invertebrate communities is known to be sensitive to drought32. However, we did not detect statistically significant changes in community composition along the elevation gradient (df = 8, dev = 359.75, p = 0.15), or among different vegetation types (df = 6, dev = 692.15, p = 0.08). This is consistent with observed high abundances of soil invertebrates at high altitudes, independent of vegetation33. Nonetheless, deeper sequencing efforts might be needed to reveal these effects.
As a second example, we estimated the genome size for 191 species using the assembly-based approach ModEst34. We found a 30-fold range of genome sizes across the groups (Fig. 1), from 79 Mb (the nematode Discolaimus major) to 2.9 Gb (the chilopod Lithobius crassipesoides). Nematoda and Tardigrada had typically small genomes, whereas the genomes of Enchytraeidae were remarkably larger. In addition to between-group variation, some groups also had a wide range of genome sizes among member species. For example, Chilopoda (centipedes) genomes ranged in size from 0.178 to 2.90 Gb, while Oribatida genomes ranged from 0.09 to 1.72 Gb. Repeat content and GC content also varied widely both within and between soil invertebrate groups (Supp. Figure 1).
Next to mutation rate and the nucleotypic theory, numerous further attempts exist to infer causal relationships of genome properties and functional traits (Fig. 3a)35–37. Effective population size should be influenced by habitat generalism, with the expectation that species able to thrive in a wide range of habitats should have larger population sizes and therefore also larger effective population sizes (Ne)38. Ne should be inversely related to body size, as larger populations of small-bodied organisms can be maintained by the same amount of resources in comparison to large-bodied species36. The reproductive mode is known to impact Ne, because the higher the degree of inbreeding, the smaller the expected Ne39. High Ne is frequently hypothesised to contribute to reducing repeats as evolutionary burdens from genomes, as selection is more efficient in larger populations35,37. Repeats are frequently considered to increase genome size40,41. If the repeats themselves are biased in base composition, this should reflect in the overall GC content. Interestingly, GC content is also linked to resource availability42, which may be linked to habitat generalism via higher metabolic flexibility43. Even though most of these observations originate from bacterial studies, ample evidence exists that the environment may influence base composition also in metazoans42, 44–46.
We tested these hypotheses simultaneously with a series of structural equation models (SEMs, Fig. 3a). To parametrize the SEM and connect the 232 new genome assemblies with ecological traits, we first gathered trait data from original literature. Information about habitat preferences was added from the Edaphobase data warehouse for soil biodiversity (https://portal.edaphobase.org/). We focussed on three traits: a) body length as a proxy for body size (minimum female adult body length for nematodes, and mean adult body length for all other taxa), b) reproduction mode, and c) the number of known habitat types where a species occurs, as a proxy of habitat generality (based on CORINE - Coordination of Information on the Environment47). We annotated repetitive elements with species-specific repeat libraries. We estimated effective population size (theta) directly from the genome data by making use of the genome-wide heterozygosity in the reference individual. This proxy measure of effective population size was calculated individually for each genome assembly with at least 4X coverage. Genomic and ecological traits are accessible in Supp. Table 1.
The SEMs supported only two of the hypothesised causal relationships when these were modelled for all taxa together (Fig. 3b): a positive link between repeat content and genome size, and a negative link between repeat content and GC content. Genome size is frequently considered to be driven by repeat content48,49, but with variation in the relationship among higher taxa of vertebrates50. Higher GC content is linked to smaller genome size in many but not all eukaryotic groups48. This link might also originate from the expansion of repeats with low GC content. Our results confirm that the strength and direction of relationships among genome size, repeat content, GC content, and their ecological drivers vary among higher taxa of invertebrates (Supp. Figure 2).
When modelling all soil invertebrate groups together, other hypothesised causal relationships were either statistically insignificant or pointed to the opposite directions than theory predicted (Supp. Figure 2). Most strikingly, high Ne size was linked to higher repeat content which in turn implies larger genome size. This suggests that efficient selection does not universally reduce the evolutionary burden of large genomes and repeat content51. Taxon-specific processes might be more important for shaping repeat landscapes than demography52. Such processes might include epigenetic regulation via repetitive elements and maintenance of chromosome structure53, and modification of gene expression and transcript diversification54. Indeed, SEMs fitted separately to higher taxa (myriapods, oribatid mites, springtails, nematodes) showed marked group-specific differences in the support of causal hypotheses between genomic and ecological traits (Fig. 3c-f, Supp. Figure 2). The link between body size and genome size received statistical support only in myriapods, but with an opposite sign as predicted by the nucleotypic theory36: several smaller-bodied species of the genera Lithobius and Schendyla had considerably larger genomes than larger-bodied species. This indicates lineage-specific expansion or contraction of genomes, reported for diverse eukaryotes50,55. Genome size evolution is complex and the effects of body size are often difficult to disentangle from other co-variables56,57. Habitat generalism was positively linked to effective population size in myriapods and oribatids, but not in springtails and nematodes: generalists may not be as fit as specialists in any particular habitats58,59, and their evolution might depend on differential rates of population evolution compared to rates of environmental change60. Altogether, the SEM results do not support the taxonomic generality of diverse hypotheses on genome evolution, with divergent patterns reflecting the distinct evolutionary history of the different groups.
As a third example, we explored whether shared gene loss might be related to repeated adaptations of phylogenetically distant metazoans to soil conditions. Gene loss is a key process in evolution61,62. Here, the dense taxon sampling for individual groups allowed to differentiate between consistent gene absence across several taxa, which likely indicates gene loss, and the sporadic absence of a gene in individual taxa, which likely represents noise introduced by assembly incompleteness. To further reduce the risk that incomplete gene annotations generate a spurious signal of gene loss, we used a targeted search for orthologs in the un-annotated genome assemblies to determine the presence/absence patterns of genes across taxa. We analysed the presence of 1482 core metazoan gene orthologs. Notably, this revealed that 50 core genes are missing in springtails (n = 78 species), and 97 core genes were not found in the oribatid mites (n = 54 species) (Fig. 4b). Given the large number of investigated taxa in the two groups, it is unlikely that these genes have been accidentally missed. Instead, their absence indicates gene losses early during diversification of the respective groups, which makes it likely that they are connected with the evolutionary emergence of the specific phenotypes of springtails and oribatids, similar to what has been seen for other animals63. Overall, fifteen gene ontology terms were significantly enriched (testFisher < 0.05) among the genes lost involving biological processes such as tubulin metabolism and cellular and subcellular movement (Oribatida). There was a significant loss of genes involved in pyridin-containing compound metabolic processes in springtails (Fig. 4c,d; Supp. Table 2). Pyridine-containing molecules have a considerable spectrum of antimicrobial and antiviral activities64, and associated gene loss might be related to the gain of endogenous antibiotic synthesis ability by many springtail species65. However, it is difficult to connect most functional descriptions of lost genes to adaptations to soil conditions. This highlights the general difficulties with transferring annotations gained from a few model taxa to the breadth of biodiversity, with targeted annotation of specific genes being a solution.
In summary, our large collection of soil invertebrate genomes is a first major step towards a comprehensive DNA- or RNA-based identification of the entire soil biodiversity: they extend the scope of metagenomic or metatranscriptomic studies from microorganisms to metazoans. Further, the taxonomically broad and dense sampling of genomes provides unique insights into genome evolution. Here we could show that no single theory of genome evolution fits all taxa: there are probably no simple overarching explanations for observed variations in genome properties, but interactions of multiple drivers result in divergent genome evolution patterns in different groups, reflecting their unique evolutionary history. Broad genome sampling allows for the identification of group-specific gene loss patterns, highlighting issues and future directions around the functional annotation of genomes from non-model taxa in diverse habitats. Overall, the 232 soil invertebrate genomes demonstrate the importance of genome sequencing efforts for understanding the ecology and evolution of the full scale of eukaryotic biodiversity, and project a future when maximum taxonomic and functional information will be gained from every environmental DNA or RNA fragment.