Mobile genetic elements shape microbial diversity and functions in thawing permafrost soils

Microbiomes are increasingly recognized for driving biogeochemical and nutrient cycling across diverse systems^1,2,3. Metagenomic advances have sharpened our understanding of microbial diversity in these systems, from the low resolution of methods such as denaturing gradient gel electrophoresis and 16S amplicon sequencing (reviewed in ref. ⁴) to metagenome-assembled genomes (MAGs) recovered from shotgun sequencing⁵. Recovering MAGs typically uses metagenome binning approaches that collapse diversity within a population, omitting variants from the consensus⁶. New tools are bringing metagenomic strain-level diversity into focus in intensively studied systems such as the human microbiome, but still struggle in complex or understudied systems such as soils^7,8,9. Yet, within-species variation can have profound effects on ecosystem outcomes (for example, ref. ¹⁰). Such microdiversity can arise through both mutation and gene flow¹¹ (reviewed in ref. ¹²). Understanding the constraints on these processes is a critical step towards understanding the pace and scope of adaptive change accessible to microbiomes in changing environments.

Mobile genetic elements (MGEs) are a major driver of gene flow and genomic rearrangement^13,14. Recombinase genes enable MGEs to move between genomic loci and can be used to distinguish between several MGE types¹⁵: insertion sequences and transposons, phage and phage-like elements, conjugative elements and mobility islands, and integrons. These types differ in several respects, including their ability to move between cells unaided and the likelihood that, when they move, they carry genes that are not required for MGE mobility but whose expression may alter host metabolism or function (‘cargo genes’). MGEs can influence host phenotypes and niches through their cargo and/or the genomic context of MGE insertion (Fig. 1a), causing gains or losses of function and/or regulatory shifts^16,17 that should scale and diversify with the number of such movements (Fig. 1b). The wider an MGE’s host range is, the greater is its range of partners—both hosts and other MGEs—for exchanging genetic modules¹⁸.

Yet, so far, MGEs have not been accessible to comprehensive characterization in a natural, non-host-associated community. With comparatively few reference genomes available in such communities, analysis depends on de novo assembly, but the mobility of MGEs tends to break assembly graphs^19,20; where binning succeeds, the loss of strain-level variation in consensus sequences hides differences in MGE presence, absence or genomic context. Moreover, few long-term projects sequence deeply enough or often enough to track MGE activity. These structural challenges have made MGEs a persistent unknown in environmental microbiology. Targeted studies have grappled with some aspects of MGE biology, for example, transfers of a given MGE-associated trait (such as antimicrobial resistance^15,21) or MGE type^22,23,24, or MGE identification in the genomes of isolates¹⁵. But so far, there is no comprehensive survey of MGEs in a complex natural system, and thus no quantitative basis for bounding their potential impacts. Addressing this gap requires deeply sequencing the same microbiome repeatedly over short and long timescales.

Long-term ecological study sites provide an opportunity to capture such datasets. The model permafrost peatland Stordalen Mire, Sweden, lies in the rapidly changing discontinuous permafrost margin, where rising temperatures convert the intact palsa habitat to partially thawed bog and fully thawed fen, with concomitant increases in methane efflux from the soil^25,26,27 (Fig. 2a). Sampling from active-layer soil in these three habitats over almost a decade of peak growing seasons (2010–2017) has yielded 615 short-read metagenomes²⁸, of which 50 have paired short-read metatranscriptomes. Previous analyses of a subset of this data have identified key taxa and processes²⁶ and have pointed to the role of functional redundancy and dispersal in allowing each habitat’s microbiome to resist environmental change^26,29. MGE-mediated gene flow could be synergistic with both, but so far, even bounding the potential impacts of MGEs on a natural system has been intractable.

Here, we use the Stordalen Mire dataset to develop and apply the needed analytic approach. We use paired short- and long-read sequencing to determine the expected yields in MGE recovery from short-read data, and use these rates to obtain conservative bounds on the per-genome abundance of several MGE types across different phyla. We assess the functions of genes affected by MGEs, revealing that MGE-mediated gene gain and loss reaches many functions in core nutrient metabolism, with implications for community metabolism and ecosystem outputs. Finally, we present three complementary metatranscriptomic and metagenomic metrics of MGE mobility at different timescales and across MGE types. Collectively, these analyses shed light on the scale and potential impacts of MGE activity in a soil microbiome and provide a set of approaches that will enable quantitative investigation of MGEs in short-read metagenomic data from complex natural systems.

MGE recombinases from permafrost thaw gradient microbiomes

Stordalen Mire is a climate-critical permafrost peatland, extensively studied for decades to assess the interplay of vegetation, microbiota and climate along a three-stage thaw gradient of palsa, bog and fen habitats (least to most thawed, respectively). Before our work, ~5.8 terabases of metagenomic and metatranscriptomic sequencing data had been generated from 8 years of field sampling focused on the active layer across habitats (Fig. 2a and Supplementary Table 1). We first identified and quantified all major MGEs in these data via recently developed hidden Markov models (HMMs) and identification protocols that leverage recombinases as MGE hallmark genes¹⁵ (Supplementary Table 2). Importantly, because this approach relies on recombinase detection, it does not fully capture plasmid diversity, but focuses instead on integrative MGEs. We screened 131.5M genes predicted from 30.7M short-read metagenomic contigs, identifying 2,107,204 MGE recombinases (1,130,776 unique protein sequences) on 1,790,672 contigs (Fig. 2b–d, Extended Data Figs. 1a,b and 2a,b, Supplementary Fig. 1 and Supplementary Table 3, item A). On average, this represented 17.3 MGE recombinases per million base pairs (Mbp) of microbiome sequence data, dwarfing the representation found for even the most abundant Kyoto Encyclopedia of Genes and Genomes (KEGG) gene family, the stress-response sigma factor rpoE (2.4 copies per Mbp) (Fig. 2c). This is consistent with previous analyses of cultured microorganisms¹⁵ and reports of recombinase prevalence and abundance in microbiomes³⁰.

We next classified these abundant recombinases into four previously-established types¹⁵: (1) integrons, (2) insertion sequences and transposons (‘IS_Tn’), (3) phage/phage-like elements (‘phage’) and (4) conjugative elements and mobility islands (‘CE’). This classification facilitates the evaluation of ecological impacts, as these types are known to differ in length, mode of propagation and gene content. Prior work evaluating 187 metagenomes (available in 2010) and >80K bacterial and archaeal isolate genomes^15,30 found IS_Tn was dominant in highly abundant MGE pools. In Stordalen Mire, we also found that IS_Tn was the most commonly recovered type (58.6%), distantly followed by phage (16.1%), CE (2.8%) and integron (1.2%) (the remaining 21.3% were ambiguous) (Methods and Extended Data Fig. 1a). Ecologically, this rank-order pattern of recombinase types remained remarkably constant across time (8 years of peak growing season sampling) and space (within and between the study site’s three habitat types) (Fig. 2d and Extended Data Fig. 1b,c). The same pattern largely persisted across phyla for both MAG- and contig-level taxonomic assignments (available for 216,720 and 804,456 MGE recombinases, respectively, assigned to a total of 89 phyla with 93.7% concordance for the subset shared) (Fig. 2a,b, Extended Data Fig. 3a,b and Supplementary Table 3, item A). MGE distribution broadly followed host lineage relative abundances, with most MGEs in ecosystem-dominant bacterial and archaeal phyla such as Acidobacteriota, Actinomycetota, Pseudomonadota, Verrucomicrobiota and Halobacteriota²⁶ (Fig. 3 and Extended Data Figs. 3 and 4). Echoing isolate-based observations¹⁵, the virtual absence of integrons in the abundant Actinomycetota (Extended Data Fig. 3, ~0.2 per genome in this lineage compared with ~1.0 on average for all lineages) stands out against a background of low cross-phylum variation in the frequency of each MGE type in our Stordalen microbiomes.

Next, we wondered how per-genome MGE frequency, previously measured in isolates, would compare with estimates from natural populations, as MGEs may be lost during laboratory cultivation owing to reduced selective pressure in clonal isolates. We quantified MGEs per genome in Stordalen microbiomes via two approaches. First, we used a contig-based estimate, comparing identified MGE recombinase counts with total genome counts estimated using a standard single-copy marker gene approach for each sample^31,32. This revealed a median of 36 MGE recombinases per genome (Fig. 2e). Separately, we considered only the subset of MGEs that were binned into medium/high-quality MAGs (≥70% completeness, ≤10% contamination) (Supplementary Table 3, item B). Scaling per-MAG MGE counts by MAG completeness revealed a median of 23 MGE recombinases per genome (Fig. 2e). Comparing the two approaches, we note that only 15.9% of MGE-bearing contigs were binned, even though Stordalen MAGs capture a remarkable 73% of sequence reads at the genus level^26,29. Thus, we interpret the MAG-based estimate as a lower bound, as it empirically captures only a subset of the contig-observed MGEs, and dynamic rearrangements in MGE-bearing genomic regions are a known challenge for short-read metagenome assembly and binning^20,33.

To empirically quantify the impact of assembly and binning challenges on short-read counts, we performed long-read sequencing on six samples and evaluated MGE recombinase recovery throughout assembly and binning. In total, these data suggested that short-read assembly and binning recovered 35–60% of MGE recombinases, with IS_Tn recombinases disproportionately under-recovered (Fig. 2f, Extended Data Figs. 5 and 6, and Supplementary Text). Using this benchmarking to scale our MAG-based counts, we obtain a revised estimate of 38–66 MGE recombinases per genome in situ. We interpret this number as comparable with the contig-based estimate and consistent with the hypothesis that isolates might lose MGEs during laboratory cultivation.

MGE impacts range far beyond biotic social interactions

We next wondered whether MGEs affected genes of ecological importance, and how these genes fit into microbiome activity in Stordalen Mire. Prior genome-resolved metagenomic analyses^26,29 show that intercellular interactions and public goods matter in Stordalen Mire: cross-feeding relationships between carbon generalists and methanogens are persistent and strongly associated with methane flux, and extracellular degradation of polysaccharides is a bottleneck in carbon flow. Beyond microorganisms, abundant viral community members harbour auxiliary metabolic genes that probably mediate exopolysaccharide breakdown²⁷. Such outward-directed ‘social’ functions are expected to be affected by MGEs^34,35,36, as are some stress responses³⁷, but these are only a subset of microbial activities. Each cell must also manage its own routine ‘domestic’ processes (for example, genome maintenance, gene expression and nutrient homeostasis), functions that shape physiology and ecosystem outputs but that have been little studied in the mobilome so far. Recently observed plasmid-borne genes affecting domestic activities such as nitrogen fixation³⁷ highlight the potential for MGEs to influence nutrient cycling directly, but the prevalence of such functions among MGE-affected genes remains unknown.

Thus, we asked which Stordalen Mire functions could be confidently detected as mobilome-affected (that is, encoded in MGEs or interrupted by them) and how far beyond biotic interactions the mobilome’s direct influence extends. To ensure that host gene activities were not misattributed to MGEs, we adopted a highly conservative approach for curating MGE boundaries (see the Supplementary Text section ‘Conservative MGE detection and impacted gene curation’). Notably, this stringent approach excluded all but ~0.8% of our 2.1M MGE recombinases from functional gene analysis, precluding the investigation of MGE-mediated changes in metabolic capacity over time. Annotation identified putative functions for cargo genes within conservative boundaries of 5,288 MGEs (900 CEs, 551 integrons and 3,837 phages) and for host genes interrupted by 11,295 IS_Tns (Methods and Supplementary Fig. 1). Among these, manual curation identified 3,631 MGE-affected genes with high-confidence functional predictions across 1,199 distinct KEGG annotations in diverse categories including carbohydrate metabolism, transport, antiviral defense, signalling, transcription, post-translation modification, lipid metabolism, nucleotide metabolism, DNA repair and antibiotic metabolism (Fig. 4, Supplementary Tables 3 (items C and F) and 4, and Supplementary Text). Here, we highlight two ecosystem-critical processes: carbon cycling and salt stress response.

We uncovered MGEs affecting central carbon metabolism genes (6 carried and 45 interrupted) (Supplementary Table 4), including those encoding enzymes of the tricarboxylic acid (TCA) and reductive TCA cycles and the pentose phosphate shunt. Alcohol dehydrogenases (1 carried and 27 interrupted) and pyruvate metabolism enzymes (5 interrupted) were well represented (Supplementary Table 4); because these enzymes’ substrates can be metabolic branchpoints, MGE activity may alter the partitioning of fluxes between pathways with different energy yields. Intriguingly, four of the interrupted alcohol dehydrogenases were in Actinomycetota, where methanotrophy sometimes uses members of this enzyme family^38,39. We also identified two direct connections to methane metabolism: the key hydrogenotrophic metabolism gene formate dehydrogenase⁴⁰ was carried by a phage-like MGE from Halobacteriota (Methanoregula sp003136075), a group of methanogens active at Stordalen Mire^25,26,29, and an IS_Tn-interrupted methane monooxygenase (pmoC) subunit in a Pseudomonadota (unbinned Alphaproteobacterial contig). Thus, even within our highly conservative subset of MGE-affected genes, we find evidence that MGEs may alter Stordalen methane fluxes both directly (by changing genomic potential for methanogenesis and methanotrophy) and indirectly (by changing the distribution of carbon-processing pathways that control the flow of substrates to methanogens).

Salt stress is common in permafrost-affected soils owing to brine film formation around soil when water freezes⁴¹. Cells can respond by modulating membrane stability and permeability, or by synthesizing and/or transporting osmolytes⁴². We found MGEs impacting genes that may contribute to each of these responses. Membrane-tuning cyclic fatty acids vary in abundance between the microbiota of nearby Sphagnum (bog) and Carex (fen) peats⁴³, suggesting that the adaptive value of cyclopropane-fatty-acyl-phospholipid synthase (CFAS; two carried and three interrupted; one in bog and four in palsa) may vary with thaw stage. MGEs also affected genes involved in the synthesis or transport of numerous osmolytes, including carnitine (gamma-butyrobetaine hydroxylase, one interrupted), trehalose (malto-oligosyl trehalose synthase, five interrupted), sucrose (sucrose-6F-phosphate phosphohydrolase, one carried), proline (iminopeptidase, one interrupted; proline-specific peptidase, three carried) and glutamate (formiminoglutamate hydrolase, one interrupted; glutamate ABC transporter ATP-binding protein, two carried). However, the involvement of several of these compounds in central metabolism complicates interpretation.

Critically, even though our stringent analysis criteria and the limitations of functional gene annotation in understudied environments together reduced the set of MGE-affected functions analysed to genes affected by 0.02% of the total MGEs identified, we still find abundant examples of mobilome-affected social and domestic functions beyond previous observations⁴⁴. Our analysis adds examples of MGE-impacted functions that influence carbon flow through this climate-critical ecosystem, both directly (exopolysaccharide degradation and central carbon metabolism) and indirectly (via nitrogen, sulfur and/or iron availability to support carbon flux). Moreover, even the small subset of MGE-impacted genes identified here may have large regulatory impacts⁴⁵. For example, we identified many MGEs affecting sigma factors (1 carried and 45 interrupted) of the stress-response sigma-24 (rpoE) (Fig. 4) and housekeeping sigma-70 (rpoD) families (Supplementary Table 4). Such master regulators’ large regulons could indirectly amplify MGE impacts. Despite prior work’s focus on MGE impacts on biotic interactions, in our analysis, ~55% of MGEs affect at least one domestic function, and ~51% of MGE-affected functions are domestic. If this conservative subset is representative of the Stordalen mobilome as a whole, then as sequencing methods, functional annotation and databases improve, we would expect ~1.1M of the full ~2.1M MGE recombinases to influence domestic functions, with important implications for nutrient cycling, as well as microbiome activity, ecology and evolution.

MGE recombinase mobility and activity

To assess the ecological and evolutionary change in Stordalen Mire MGEs, we evaluated their mobility and activity from three different viewpoints: (1) snapshots of metatranscriptomic MGE recombinase expression as a proxy for mobilization underway, (2) genomic neighbourhood change as a marker of completed MGE movement between insertion sites and (3) partial carriage of an MGE by the host population as the result of evolution and selection over longer timescales. Integrating results across these analyses can help highlight individual MGE dynamics. For example, high recombinase expression together with low rates of genomic neighbourhood change for a given MGE could suggest a strong preference among insertion sites. Our overall findings are as follows.

First, we assessed MGE recombinase transcriptional activity by read mapping from 50 available metatranscriptomes against contigs from paired metagenomes (Supplementary Table 5), taking ribosomal protein-encoding gene coverage from read recruitment as a baseline for constitutive transcriptional activity of actively growing cells⁴⁶ (Supplementary Table 6). This revealed detectable transcriptional activity in ~3.4–56.4% (median 27.8%) of ribosomal protein-encoding genes and ~0.04–9.1% (median 1.2%) of MGE recombinase genes (Fig. 5a and Supplementary Tables 3 (item G) and 7), making MGE recombinases the second-most-transcribed non-ribosomal genes in Stordalen Mire (Fig. 5a). We found roughly equal levels of transcriptional activity for IS_Tn and phage recombinases, approximately five- to tenfold higher than for integron and CE recombinases (Fig. 5b). Ecologically, IS_Tn recombinase transcription increased significantly with thaw from palsa to fen (P ≤ 0.01, Wilcoxon rank-sum test) (Fig. 5c and Extended Data Fig. 7b), suggesting an increased MGE activity in the more-disturbed fen. Taxonomically, transcriptionally active MGE recombinases were detected across hosts spanning 30 phyla, with activity largely proportional to MGE and phylum abundance (Fig. 5d and Extended Data Fig. 7c). The archaeal phylum Halobacteriota was an exception, although neither abundant nor MGE recombinase-rich it was the most active phylum for MGE recombinase transcription. Given that Stordalen Mire Halobacteriota are known to be methanogens that are active and selected for in thawed fen and bog habitats^25,47, we speculate that their disproportionate MGE recombinase activity might be adaptive.

Second, for MGE recombinases observed multiple times, we assessed the genomic context (300-bp upstream and downstream) surrounding each instance. We interpreted changing context as a proxy for intra-/inter-genomic mobility, either of the recombinase itself or of other recombinase-proximal genes (Methods; see examples of changing neighbourhoods in Extended Data Fig. 8). From the original ~2.1M MGE recombinases, we identified 74,257 MGE recombinase clusters (100% average amino acid identity (AAI)) observed more than once in our dataset. Among these, 0.1–10% showed neighbourhood variation, with rates depending on MGE type (Fig. 5b and Supplementary Table 3, item H). Neighbourhood variation was much more common among IS_Tn recombinases (~10%) than among phage (~1%), integron (0.2%) and CE (0.1%) recombinases. This was consistent among years and habitats, with a larger variation detected across diverse host phyla (up to 37% in Desulfobacterota_G) (Extended Data Fig. 9). Mindful that hypervariable islands in cellular genomes can often be impacted by multiple MGE types, we quantified co-occurring IS_Tn insertions, which could cause the neighbourhood to change while the focal MGE remains stable. Fewer than one in five of our changed-genomic-neighbourhood examples was due to co-occurring IS_Tn insertions (Extended Data Fig. 8). Because our approach only detects neighbourhood variation immediately flanking an MGE recombinase, this is a lower-bound estimate of mobility at hypervariable multiple-MGE genomic islands among MGE context changes at Stordalen.

Third, we sought to assess what fraction of a population contains any given MGE recombinase by evaluating read-mapping-based coverage of the recombinase against that of its surrounding genomic context (Methods). Where coverage drops significantly within an MGE recombinase, we infer that only a subset of the host population encodes the MGE recombinase. Overall, across MGE types, we found 8,033 recombinase clusters (100% AAI) that had ≥20 kbp of genomic context to support coverage analysis. Of these, ~5–23% showed coverage drops consistent with only a subset of the cellular population having the MGE (Fig. 5b and Supplementary Table 3, item I). These findings were again broadly consistent across years, habitats and well-sampled lineages (Extended Data Fig. 10b–d). In spite of spatio-temporal consistencies, there were notable variations in intrapopulation penetrance between MGE types, with IS_Tn and phage types particularly prone to varying across individuals within a host population (Extended Data Fig. 10d). Because MGE sizes vary, and read mapping for lesser-sampled taxa is noisy, these numbers represent conservative lower bounds on active MGE penetrance into cellular population genomes. Even so, the extent of partial carriage observed highlights MGEs’ role in what amounts to large ongoing in situ knockout experiments: among 372 IS_Tn elements that both interrupt a protein-coding gene and were amenable to partial carriage analysis, 81 (~22%) variably interrupt or leave intact host protein-coding genes. Together, these results suggest that MGE penetrance within member genomes of any given population differentially provide individuals access to genetic diversity, while shielding the population from loss of function due to MGE insertion^48,49.

MGEs have long been known from cultivated isolate genomes to be important drivers of strain-level variation and key functional adaptations, particularly in biomedical contexts related to antimicrobial resistance^15,21. However, virtually no data exist for diverse MGE types in natural settings. Our analyses presented here for Stordalen Mire microbiomes illuminate the mechanisms and scope of MGE impact on a climate-sensitive ecosystem. We find that, well beyond biotic interactions, the mobilome affects functions intrinsic to a wide range of physiologically and ecologically important processes such as carbon cycling and salt stress response. Further, active-layer soil in all three habitats shows persistent MGE activity/mobility, which could be an important force driving host strain variation and adaptation in this thawing permafrost soil system. Strikingly, the IS_Tn elements that heavily dominate recovered MGE recombinases in our dataset are also the most active by two of our three measures (that is, transcription and genomic neighbourhood change) (Fig. 5b), highlighting the need for future studies to examine context-dependent MGE effects such as modulation of host gene expression downstream of MGE insertion sites. This persistent activity occurs against a backdrop of overall MGE stability, agreeing with a companion study showing the same microbiomes are stable across the years studied here²⁹. Importantly, the limitations of short-read sequencing and assembly could drive systematic underestimates of MGE-mediated variation, through low recovery either of variable MGEs or of rare or fast-disappearing subpopulations/lineages experiencing non-selective excision/recombination events. To assess this, we confirmed and quantitatively evaluated the bias of short-read versus long-read metagenomic MGE recovery, which helped establish our identified mobilome identifications and ecosystem impacts as a lower-bound estimate. However, the framework and analytics established here are readily transferable to long-read sequencing datasets for future MGE studies.

As improved methods for studying the genomic potential of microorganisms in nature come on-line, the field has progressed from community-level taxonomic and functional descriptions towards species-resolved links between taxa and their functions^26,31. Continued advances in sequencing and analytical approaches are increasing this resolution to reveal intrapopulation variation, early speciation and selection⁷, with recent MGE typing efforts¹⁵ and further analytical advances made here finally bringing a quantitative understanding of the mobilome within reach, even in complex natural systems. This work benchmarks sequencing artefacts and provides empirical baselines and inferred impacts for all major integrative MGEs in a natural soil system. These findings suggest that ~1% of MGE recombinases are active, that genomes contain 36–66 MGE recombinases on average and that more than half of MGEs impact genes affecting basic cell physiology and ecosystem outputs. As advancing technologies bring the mobilome more fully into focus, given the toolkit and baseline findings presented here, we predict that mapping the mobilome will help reveal subtle alterations that could serve as the signal for early detection of systemic change.

Data collection, metagenome assembly and binning, gene annotation and taxonomy assignment

Our study site, Stordalen Mire, is a long-term ecological research site in northern Sweden. Large sampling and sequencing efforts since 2010 have yielded ~4.9 Tb of metagenomic sequence data from 615 readsets from 408 bulk soil samples that span space (across the site’s three habitat types: n = 188 palsa, n = 200 bog and n = 227 fen; average of ~8 Gb per sample) and time (summer peak growing season, 2010–2017) (Supplementary Table 1) with details described in another parallel study²⁹, in previous studies^26,27,50 and summarized below. The IsoGenie Consortium and EMERGE Biology Integration Institute have developed extensive resources from Stordalen Mire metagenome data, including assembly across all samples to produce 30,710,180 contigs (≥1,500 bp) and 36,419 medium- to high-quality MAGs (completeness ≥70% and contamination ≤10%; CheckM2 v1.0.2)⁵¹ that clustered into 2,048 populations (95% average nucleotide identity (ANI)) per community standards⁵². Select samples have also been sequenced using metatranscriptomics to assess gene expression (50 spatio-temporally diverse samples; average of ~18 Gb per sample). For this work, we additionally performed paired long- and short-read DNA sequencing on six samples, chosen to be representative of the larger dataset, as a control to assess short-read MGE recovery.

The protocol used for sample preparation and sequencing for short-read metagenomes has been described previously^26,27,50. In brief, nucleic acids were extracted according to a published protocol⁵³, except for the MAG source samples from substrate addition incubations (‘201607_SubstrateInc_9to19’), which were extracted with the DNeasy PowerSoil Kit (cat. no. 12988-100, Qiagen). Libraries were prepared with either KAPA Hyperprep, TruSeq or Nextera XT kits for short-read metagenomes; SMRTbell Express Template Prep Kit 2.0 for long-read metagenomes; or TruSeq Stranded mRNA kit for metatranscriptomes. Reads were sequenced with either Illumina HiSeq (2000 or 2500), NovaSeq 6000 or NextSeq 500 for short-read metagenomes; PacBio Sequel IIe for long-read metagenomes; or NovaSeq 6000 for metatranscriptomes.

Short-read sequencing data from individual samples were trimmed by Trimmomatic (v.0.36)⁵⁴, assembled by SPAdes with the ‘–meta’ option (v.3.12)⁵⁵ and then binned by MetaBAT2⁵⁶, GroopM2⁵⁷ and MaxBin2⁵⁸ (Cronin.v2), or an ensemble approach²⁹. Those sequenced at Joint Genome Institute (JGI) were processed by JGI’s standard assembly and binning pipeline. In brief, reads were quality controlled, filtered and trimmed using the BBTools package⁵⁹, assembled with SPAdes with the ‘–meta’ option⁵⁵ and binned with MetaBAT2⁵⁶. The bins produced from this effort were combined with 9,406 bins from ref. ²⁶ and 11,120 bins from stable isotope probing experiments with field peat soil and labelled litter (sequenced with NovaSeq at JGI and assembled with MEGAHIT v1.1.3)⁶⁰ and then filtered using 70% completeness and 10% contamination CheckM2 (v1.0.2)⁵¹ cut-offs to form the 36,419 MAGs (consisting of 4,193 from Cronin.v1, 27,423 from Cronin.v2, 847 from JGI, 1,543 from ref. ²⁶ and 2,413 from stable isotope probing). Next, bins obtained from the same metagenome were de-duplicated using Galah (95% ANI, commit e9e59a8)⁵², and all bins across the entire dataset were further dereplicated at 95% ANI using Galah (commit e9e59a8), producing 2,048 population level clusters. Further, the contigs encoding a recombinase were annotated with DRAM (–min_contig_size 1500)⁶¹ and taxonomically assigned at phylum level with CAT⁶², requiring a minimal contig length of 3 kb. The contigs in MAGs were assigned to the MAG taxon (obtained with GTDB-tk)⁶³, as taxonomic assignments based on whole MAGs are typically more reliable than single-contig assignments.

MGE recombinase identification

Recombinase identification was on the basis of a framework recently developed to robustly identify and classify MGE recombinases from microbial isolate genomes¹⁵. Here, we applied the same recombinase HMM profiles and suggested cut-offs to identify MGEs in metagenomic assemblies. In brief, open reading frames (ORFs) were predicted from assembled contigs by prodigal (v2.6.3) using the ‘-p meta’ option⁶⁴; then, recombinases were identified among the predicted protein sequences by hmmsearch (v3.3.2)⁶⁵ using the 68 recombinase subfamily HMMs (Supplementary Table 2) and the ‘–cut-ga’ option. Identified recombinases were assigned to the respective MGE recombinase type corresponding to the HMM with the best score, following previous recommendations¹⁵.

Annotation-based inclusion/exclusion criteria

MGE recombinase candidates were further quality controlled, as recommended previously¹⁵, by using eggnog-mapper (v2.1.9, default settings)⁶⁶. First, genes annotated as invertase, replication protein, DNA repair, DNA binding, homologous recombination, viral exonucleases or reverse transcriptase (interpreted to be spurious hits to other host genes related to DNA manipulation) were removed, except for the ones also annotated with the Pfam catalytic domains of ‘transposase’ or ‘helicase’ that were retained after manual curation. This removed 0.68% of the identified possible recombinases (Supplementary Fig. 1). Second, sequences annotated by eggnog-mapper as CRISPR, defense response to virus or maintenance of DNA repeat elements were removed to filter out spurious hits thought to be similar to Cas family recombinases, but not considered as an MGE recombinase. This removed 0.46% of the identified possible recombinases. Third, we discarded predicted recombinases that contained only a helix–turn–helix (DNA binding) domain, as well as relaxase hits for which Pfam annotates a helix domain, but without the catalytic domain required for recombinase activity. This removed 0.07% of the identified possible recombinases. Fourth, recombinases matching a ‘cell’ type recombinase HMM (and thus unlikely to be encoded by an MGE) were removed. This removed 2.52% of the identified candidate recombinases, and yielded a final total of 2,107,204 MGE recombinases (Supplementary Table 3, item A). Last, some recombinase subfamilies/HMMs include recombinases from different MGE types, which prevents assignment for that HMM to a specific MGE type. These recombinases with ambiguous origins accounted for 20.51% of the identified recombinases. They were included in the analyses of total MGE recombinases, but excluded in analyses of individual MGE types, which focused on four recombinase-based MGE types: ‘IS_Tn’ (insertion sequence and transposon), ‘CE’ (conjugative elements including conjugative element and mobility islands), ‘phage’ (phage and phage-like elements) and ‘integron’ (integron) (Supplementary Fig. 1). Note that IS_Tn recombinases integrated into the host chromosome, as opposed to into larger MGEs, cannot be distinguished using our approach and data, as defining the boundaries of the MGE is challenging in short-read assemblies, except for a small fraction of MGEs with curated boundaries (see below). To assess the diversity of these MGE recombinases, we considered recombinase sequences as ‘redundant’ when they clustered with other sequences at 100% AAI using CD-HIT (v4.8.1)⁶⁷ with ‘-c 1 -g 0 -G 1 -l 30 -d 0 -M 0’. This yielded 1,130,776 clusters that we considered as unique recombinases.

Identification of cargo genes/interrupted genes associated with MGEs

Once MGE recombinases were identified, we surveyed their genomic context to assess potential cargo and/or interrupted genes. This was done using MGE-specific approaches as follows.

Integron

All integron-containing contigs were screened for cargo genes (also known as ‘gene cassettes’) using integron-finder (v2.0.2)⁶⁸ with default parameters. In brief, integron-finder identifies integrons by using integron-specific integrase HMMs to identify integron integrases (intI), and a covariance model to identify the integron attachment site (attC), as the two bounding regions of an integron sequence. Only integrons encoding intI and at least one attC were treated as ‘complete’ by integron-finder and were retained. There are two types of integron structures that are known in terms of the orientation of intI and attC: (1) a regular type with intI and attC on opposite strands and facing outwards and (2) another type with intI and attC on the same strand, called an inverted integrase integron⁶⁸. Of the 686 ‘complete’ integrons we identified, 681 integrons also met these orientation expectations and the 5 that did not were not considered further. The gene cassettes delimited by intI and attCs in these integrons were treated as cargo genes associated with integrons⁶⁹, resulting in integron regions with various lengths (minimum 983 bp, median 2,310 bp, mean 3,092 bp and maximum 14,195 bp). Problematically, most contigs (97.4% or 25,953 out of 26,634) encoding Khedkar’s ‘integron’ type recombinases were not considered complete by integron-finder. To test this, we also replaced the integron integrase HMM in the integron-finder tool with that from the study by Khedkar et al.¹⁵ and found it produced nearly identical results and did not yield additional instances of predicted complete integrons. One explanation could be that attC, not integron integrase, is the limitation, as the flanking attC sequences are identified using a database of known attC sequences and these could miss novel attC sequences used by integrons in natural microbiomes. Given this, we interpret the large discrepancy between Khedkar’s and integron-finder’s integrons to represent one of the following scenarios: (1) degraded or domesticated integrons that lack attC sites^70,71, (2) integrons that lack attC sites owing to fragmented assembly or (3) intact and probably complete integrons that have novel attC sites. Thus our analyses will underestimate the extent of cargo genes carried by integrons in Stordalen Mire.

CE

CEs are made of three key components—a relaxosome (includes the CE-type recombinase), a conjugation couple protein and a type IV secretion system (T4SS)—that can be used as hallmark genes that are more conserved than other CE genes⁷². To identify cargo genes associated with CE, the CE hallmark genes were identified in CE recombinase-encoding contigs by conjscan implemented in MacSyFinder (v2.0rc7) with ‘–models CONJScan/Plasmids all–db-type gembase’ and ‘–models CONJScan/Chromosome all–db-type gembase’ with gene protein sequences as input^72,73. Conjscan specifically establishes eight types of CEs using hallmark genes, which are core across all eight CE types, and accessory genes that are diagnostic for specific CE types. CEs identified by conjscan from Stordalen Mire sequence data were those that matched a relaxase, a couple protein (an ATPase that couples MGE with the T4SS), a T4SS ATPase and an additional structural gene from T4SS. For the purposes of identifying cargo genes, matches of one model, MOB, were removed, as it only requires the relaxase and thus does not identify the boundaries to establish a complete CE. As a result, we identified 1,017 near-complete CEs. As recommended, we required the distance between the hallmark and accessory genes to be ≤30 genes, except for the relaxase, which could be as much as 60 genes away from any other hallmark and accessory genes⁷². In a small number (13 of 1,017) of instances conjscan identified >1 CE in a region. For these, we manually separated the identified hallmark and accessory genes into subgroups according to gene proximity and required each subgroup to include all hallmark genes (relaxase, conjugation couple protein and T4SS ATPase) and at least one accessory T4SS gene. Once such ‘complete’ CEs were identified, the genes between the leftmost and rightmost hallmark and accessory genes were assumed to be cargo genes carried by the CE (except for hallmark and accessory genes themselves). This resulted in CE regions with various lengths (minimum 7,907 bp, median 16,290 bp, mean 23,383 bp and maximum 81,148 bp). Overall, 1.6% of the 62,660 CE-type recombinase genes identified by our analyses were considered near-complete CE by conjscan. Again, this was not a function of improved HMMs as we tested replacing the relaxase HMMs in conjscan by those from Khedkar et al.¹⁵ and found it produced similar results. Instead we interpret this large discrepancy to suggest that most CEs in these data are degraded, fragmented from poor assembly or cannot be formally detected as complete because of limitation and bias in hallmark gene databases.

Insertion sequences/transposons

To identify IS_Tn recombinase-interrupted genes, up to 3,000 bp upstream and downstream of IS_Tn recombinases were collected from each contig, requiring a minimal length of at least 50 bp on each end and at least 500 bp across the two ends combined (with 1,063,800 IS_Tn recombinases in total meeting the requirements), and then aligned by minimap2 (v.2.22, default settings)⁷⁴ to the coding sequence of representatives of protein clusters (PCs) obtained by MMseqs2 at 95% identity (v4.8.1)⁷⁵ from all contigs (≥1,500 bp) used in this study. The corresponding coding sequence was considered interrupted if the following four conditions were met: (1) if upstream and downstream alignments combined covered ≥50% of the coding sequence of the reference PC; (2) if there is a gap between upstream and downstream alignments of the coding sequence, the gap is ≤10% of the coding sequence length; (3) if there is an overlap between upstream and downstream alignments on the coding sequence, the overlap is ≤10% of the shorter alignment; and (4) if there are multiple alignments either upstream or downstream of the coding sequence, the pair with largest coverage of coding sequence is chosen. This identified 5,342 PCs putatively interrupted by a total of 12,332 IS_Tn-type recombinases.

Phage/phage-like elements

All contigs encoding ‘phage’ recombinases were screened for temperate phage by VirSorter2 (v2.2.3) following VS2 SOP step 1 (v3.0)^76,77, and the predicted phage sequences trimmed by checkV (default setting, version 0.7.0)⁷⁸. Next, the checkV trimmed sequences encoding a ‘phage’ recombinase were considered as phage/phage-like regions. To identify cargo genes, that is, host(-like) genes encoded on phage/phage-like elements, we opted for a conservative approach in which we reduced the dataset for ‘cargo gene inspection’ to include only phage/phage-like regions with one ‘phage’ recombinase, and use the leftmost and rightmost checkV, called ‘viral gene’ or ‘phage’ recombinases, as the boundaries. This identified 4,277 temperate phages with ‘phage’ recombinases and at least one extra gene other than the two at the ends. Meanwhile, 347 phage/phage-like elements were removed for having more than one ‘phage’ recombinase, in order to avoid cases where multiple contiguous temperate phages or degraded ones yielded a single prophage prediction, which could lead to host genes being mistaken as cargo genes. This resulted in 3,930 conserved phage/phage-like regions (minimum 1,386 bp, median 7,293 bp, mean 13,779 bp and maximum 352,302 bp).

Annotations of cargo/interrupted genes

Insertion sequence or transposon interrupted genes and genes within other MGEs were functionally annotated using KEGG and Pfam databases in DRAMv for phages/phage-like elements or DRAM for other MGEs (v1.2.3,–min_contig_size 1500)⁶¹, and further evaluated by comparison with best hits in NCBI NR (downloaded on 4 October 2020) using DIAMOND (v2.0.8.146, blastp bitscore ≥50)⁷⁹. A subset of the genes of specific interest to the Stordalen Mire thawing permafrost system were selected and manually curated for display (Supplementary Table 4) including transport, carbon/nitrogen/sulfur cycling, signalling, antiviral defense, salt stress response and genome maintenance and expression.

Evaluation of short-read assembly and binning artefacts for recovering MGE recombinases using long reads

Assembly and binning of paired long reads and short reads

Soil samples (six total) were collected at the Stordalen Mire site from 2019 from two depths (centimetres below ground) across three habitats (palsa, bog and fen). DNA extraction and Illumina NovaSeq sequencing was performed as described above. In addition, DNA sequencing was performed from the same DNA extractions using the PacBio Sequel IIe HiFi sequencing platform at JGI. PacBio data were processed at JGI to form filtered circular consensus sequencing reads. Assemblies were generated through two methods: short-only and hybrid short-read plus long-read (herein, hybrid). Short-only was assembled with metaSPAdes v3.15.4⁵⁵ using Aviary v0.5.3⁸⁰ with default parameters. Hybrid assembly was performed using Aviary v0.5.3 with default parameters. This involved a step-down procedure with long-read assembly through metaFlye v2.9-b1768⁸¹, followed by short-read polishing by Racon v1.4.3⁸², Pilon v1.24⁸³ and then Racon again. Next, reads that did not map to high-quality metaFlye contigs were hybrid assembled with metaSPAdes v3.15.4 and binned out with MetaBAT2 v2.1.5⁵⁶. For each bin, the reads within the bin were hybrid assembled using Unicycler v0.4.8⁸⁴. The high-coverage metaFlye contigs and Unicycler contigs were then combined to form the assembly fasta file. Genome recovery was performed using Aviary v0.5.3 with samples chosen for differential abundance binning by Bin Chicken v0.4.2⁸⁵ using SingleM metapackage S3.0.5⁸⁶. This involved initial read mapping through CoverM v0.6.1⁸⁷ using minimap2 v2.18⁷⁴ and binning by MetaBAT v2.1.5⁶, MetaBAT2 v2.1.5⁵⁶, VAMB v3.0.2⁸⁸, SemiBin v1.3.1⁸⁹, Rosella v0.4.2⁹⁰, CONCOCT v1.1.0⁹¹ and MaxBin2 v2.2.7⁵⁸. Genomes were analysed using CheckM2 v1.0.2⁵¹ and clustered at 95% ANI using Galah v0.4.0⁵².

Assessment of MGE recombinase recovery by short-read assembly and binning

We considered the high-quality hybrid MAGs as the best available references and used these to assess the recovery of MGE recombinases in short-read contigs and MAGs. A total of 17 high-quality MAGs (completeness ≥95%, contaminant ≤5% and at most 100 contigs) from hybrid assembly were selected and paired with 17 MAGs recovered from short-read assembly through the 95% ANI clusters (the highest ANI were chosen if there were many) for comparison. The short-read MAGs included 5 medium quality (≥70% completeness, ≤10% contamination) and 12 high quality (≥95% completeness, ≤5% contamination).

Evaluation of these 17 MAGs were carried out as follows. First, each pair of hybrid MAG and short-read MAG were compared in terms of total number of MGE recombinases, as well as unique and shared MGE recombinases at 100% AAI by CD-HIT (v4.8.1) with ‘-c 1 -g 0 -G 1 -l 30 -d 0 -M 0’ (ref. ⁶⁷). Second, we evaluated the impact of sequencing depth, assembly and binning to assess which step was most problematically leading to under-detection of MGE recombinase identification in short reads. Initially, we mapped short-read assembled contigs (sheared into 500-bp fragments) to hybrid MAGs from the same sample using CoverM v0.6.1, requiring each read to be 75% aligned with 95% identity. Third, these data were used to assign MGE recombinases in the 17 hybrid MAGs into one of three categories: (1) those recovered by short-read MAGs (termed ‘covered’); (2) those not fully covered by short-read contigs, that is, those missed by short-read assembly (termed ‘missed assembly’); and (3) those fully covered by short-read contigs, that is, those recovered by short-read assembly but missed by short-read binning (termed ‘missed binning’). Fourth, the coverage of each gene was generated by dirseq v0.4.3²⁶ using the BAM file from above and used to assess the relative abundances of each category as a metric for the impact on MGE-representation in the short-read assemblies.

With a better understanding for which step was failing most often for MGE recombinase recovery in short-read assemblies, we sought next to evaluate what factor(s) might be causing the failed capture events. First, we evaluated sequencing depth, as low coverage might cause breakages in short-read assembly graphs. This was done by mapping trimmed short reads to hybrid MAGs using CoverM v0.6.1, requiring each read to be 75% aligned with 95% identity, calculating the read coverage of each gene by dirseq v0.4.3 and comparing coverages among the three categories of MGE recombinases. Second, we compared the genomic neighbourhood number among the three categories, reasoning that perhaps indels were an issue. Genomic neighbourhood number was calculated by clustering genomic neighbourhoods of each unique MGE recombinase (see below) at 90% ANI. To test this, nucleotide diversity for each base position was calculated as 1 − [(number of ‘A’ bases/total bases)² + (number of ‘C’ bases/total bases)² + (number of ‘T’ bases/total bases)² + (number of ‘G’ bases/total bases)² + (number of ‘deletion’ bases/total bases)² + (number of ‘insertion’ bases/total bases)²], adopting the nucleotide diversity definition from inStrain⁷ with the addition of deletion and insertion variants. The variant count at the base position was generated by pysamstats (v1.1.2; https://github.com/alimanfoo/pysamstats) with the ‘–type=variation’ option. Third, we assessed the nucleotide diversity of each MGE recombinase, reasoning that perhaps elevated nucleotide diversity could lead to short-read assembly graph breakages. To test this, for each of the 17 MAG pairs, MGE recombinases only identified in hybrid MAG were compared with short-read assembled contigs to track whether they were assembled by short reads, and then whether they were binned in MAGs from other MAG clusters. The tracking of MGE recombinases in contigs was accomplished with ‘MMseqs map’ (v13.45111), with default settings, and then filtered for hits with 100% identity in the protein sequence⁷⁵.

MGE recombinase per genome estimation from contigs and MAGs

To estimate the MGE recombinase number per genome from contigs, the median count of ribosomal genes was first used to estimate the number of genomes in each assembly. To do this, ribosomal genes were extracted from contig annotations by DRAM⁶¹ using KEGG IDs (Supplementary Table 6), and the average number of recombinase genes per genome was calculated by dividing the total number of recombinases in the sample by the estimated total number of genomes. We next linked individual recombinase-encoding contigs to MAGs built from the same sample. To that end, we used exact string searching to identify MGE-recombinase-encoding contigs in the dereplicated MAGs from the same sample. Through this process, 106,017 out of a total 665,700 (15.9%) MGE-recombinases-encoding contigs ≥3 kbp could be linked to a MAG. To account for partial MAG recovery, the number of MGE recombinases per genome estimated from MAGs was calculated as the number of MGE recombinases in a given MAG divided by MAG completeness.

Genomic neighbourhoods of MGE recombinases and mobility evaluation

MGE mobility was assessed by comparing the genomic neighbourhoods of individual recombinases. For this analysis, we limited the input data to recombinases that occurred ≥2 times in the total dataset (assembled contigs per sample). This was achieved by clustering the recombinases at 100% AAI using CD-HIT v4.8.1 with ‘-c 1 -g 0 -G 1 -l 30 -d 0 -M 0’ (ref. ⁶⁷), only keeping the clusters with at least two recombinases. For each MGE recombinase, 300 bp upstream and downstream of the recombinase were considered as its immediate genomic neighbourhoods, excluding those near the contig ends (<300 bp). For each recombinase cluster, a pairwise comparison of neighbourhoods was performed among all member recombinases using both neighbourhoods (that is, upstream versus upstream and downstream versus downstream) using BLAST+ (v2.8.1+) with ‘-task blastn’⁹² and a bitscore ≥50. The ANI for each neighbourhood comparison was calculated as the ‘sequence identity from blastn output’ multiplied by the ‘alignment coverages’. For each comparison, we labelled the neighbourhood with the higher ANI value as ‘high_ani’ and the neighbourhood with lower ANI as ‘low_ani’. For some MGEs, we expect the neighbourhood on the MGE side of the recombinase to be part of the MGE, and thus remain identical after integration into a new locus, so we used ‘low_ani’ values (that is, the lower ANI between the upstream and downstream region of the pair) to evaluate whether the neighbourhood of an MGE changed or not. An MGE recombinase pair was considered to be in different genomic locations if ‘low_ani’ was <1%, and a recombinase cluster was treated as active if one pairwise comparison showed different neighbourhoods. This approach only used immediate neighbourhoods (within 300 bp) as a requirement, as neighbourhoods with longer ranges would remove more recombinases from this analysis owing to the short average contig length. However, as a sensitivity analysis, we also tested the above with 600-bp and 1,000-bp neighbourhoods. These analyses showed that the size of the genomic neighbourhood (300–1,000 bp) had very little impact on the fraction of MGE recombinases identified as mobile via this changing genomic neighbourhood metric (Extended Data Fig. 9).

MGE recombinase activity by metatranscriptome

MGE activity was assessed via documenting which MGE recombinases were detectably transcribed in metatranscriptomes. Leveraging 81 paired metatranscriptomes and metagenomes, we evaluated MGE recombinase transcriptional activity by mapping metatranscriptome reads to contigs assembled from metagenomes using transcriptM v0.3.1⁹³ with ORF predictions from DRAM (see above) as the ‘–gff’ input, which also removed tRNA, rRNA, tmRNA and PhiX from the reads. To ensure the quality of RNA sequencing data (reverse stranded), we only analysed 58 pairs with ≥80% of reads mapped to the coding reverse strand of ORFs in contigs, based on dirseq v0.4.3²⁶ with BAM and GFF files from above as inputs (Supplementary Table 5). We only considered contigs ≥10 kbp to avoid spurious ORF predictions from short contigs. Next, we defined genes with ≥90% of gene length covered as expressed using pydirseq v0.1.0-alpha with BAM and GFF files from above as inputs⁹⁴. To normalize for the differences in sequencing depth and overall activity of the microbial community across samples, we used the average percentage of ribosomal genes (Supplementary Table 6) detected as active as the baseline activity to normalize other genes (termed ‘normalized active ratio’), as ribosomal genes are expected to be constitutively expressed in active cells. This ‘normalized active ratio’ thus provides an estimation of the number of genes detected as expressed within a functional category (for example, recombinases) relative to an estimated total number of genomes detected as transcriptionally active in the same sample.

MGE recombinase partial carriage by host population

As a result of MGE movement within the host population and an equilibrium of the co-evolution of MGEs and their hosts⁹⁵, we noticed that some MGE recombinases were present in only a fraction of the contigs assigned to a given host population (termed ‘partial carriage’). To quantitatively evaluate this, trimmed metagenomic reads were mapped to contigs from the same sample using CoverM v0.6.1 with default settings⁸⁷, the resulting BAM file was then filtered (requiring ≥90% of each read aligned with ≥95% identity) and used as the input for ‘samtools depth’ (v1.17) with the ‘-a’ option to get the read depth coverage at each position⁹⁶. The first and last 150 bp of each contig were not considered, as read depth coverage near the ends are expected to be artificially lower. We defined partial carriage as cases where the coverage of the MGE recombinase was at least one standard deviation lower than the larger one (requiring a minimal coverage of 1×) of the average coverages of the upstream and downstream contig regions (Cohen’s d ≥ 1). To ensure enough host regions were considered for reliable mapping-based estimation of genome coverage, the MGE recombinase total set was filtered to 23,647 MGE recombinases with ≥20 kb upstream and downstream. Note that sensitivity analyses revealed similar patterns when analysing 2.5-kbp instead of 20-kbp regions (Extended Data Fig. 10), which suggests that this conservative cut-off did not substantially influence the detection of partial carriage in these data.

We expect that this method should give us a lower-bound estimate of partial carriage, as it is not able to detect cases for which (1) the hosts have low read coverage (rare), (2) the carriage ratio is lower than but close to 1, (3) host genome coverage has a wide range of variation in different genomic locations due to strain variation within the population or (4) the MGE recombinase is on the edge of a contig.

Statistics

All statistics and visualization were performed in R v4.2.1⁹⁷. The R packages and versions used are available via GitHub⁹⁸. For statistical significance between groups, the Wilcoxon rank-sum test was used for two groups and Kruskal Wallis test was used for >2 groups as implemented in ‘compare_means’ in ggpubr package v0.6.0, with P values adjusted by Benjamini–Hochberg adjustment for multiple comparisons.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.