The study material consisted of seven myxomycete specimens from various taxonomic groups, collected by R. Hagelstein, T.E. Brooks, I. Ammirati, R. McHugh, S.L. Stephenson, Y.K. Novozhilov, and M. Schnittler (see Table 1) and currently housed in the herbarium of the Institute of Botany and Landscape Ecology at the University of Greifswald, Germany (Suppl. material 1: Sheet A). The specimens were selected to represent an age gradient, with collection dates ranging from 1935 to 1997. Initially, eight specimens were selected; however, one was found to be completely devoid of DNA for unclear reasons and was excluded from the study. Our analysis of this particular specimen is detailed in a separate publication (Leontyev and Schnittler 2025). The selected specimens differed significantly in size and number of fruiting bodies, but all contained friable, uniformly coloured spore masses and lacked visible contamination by fungi or arthropods. Around 1 mm³ of spores were sampled from each specimen; the total spore count was not normalized.

Following the recommendations of Shumskaya et al. (2023), we extracted DNA from the specimens using the filter-based QIAamp DNA Micro Kit (Qiagen, Germany), following the “Isolation of Genomic DNA from Tissues” protocol supplied by the manufacturer. A key feature of this kit is the presence of carrier RNA, i.e. polyA molecules that non-specifically bind short DNA fragments, preventing them from passing through the filter pores. For each sample, 50 µL of DNA solution were obtained. The concentration of the resulting product was measured using a DeNovix DS-11 FX fluorometer (DeNovix Inc. USA) following the TN 143 Broad Range Assay Detailed Protocol (DeNovix Inc. USA). After concentration measurement, the samples were divided in two groups. The first group, where the extracted DNA amount reached 150 ng (sc11426 and sc2833), was directly used for sequencing. For the second group (the rest of the samples, with the lower amounts of DNA), a 2.5 µL aliquot was taken for whole-genome amplification (WGA) using the 4BB TruePrime WGA Kit (4basebio, UK). The amplification process was conducted for 6 hours. Then, 35 µL of genomic DNA solution (first group) or 72 µL of WGA products (second group) with final amount of DNA ranging from 158 to 658 ng (see Table 1), were sent to Eurofins Genomics (Germany) for genome skimming on the Illumina NovaSeq (USA) platform, using the WGS Bacteria/Yeast/Fungi protocol with a sequencing depth of 5 million read pairs (2x 150 bp) per sample.

The demultiplexed, adapter-trimmed paired-end Illumina sequencing data were screened for quality and presence of remaining adapter sequences using FastQC v0.11.9 (Andrews 2010). To obtain a high sensitivity of detection of low-abundance genomic signal expected from old DNA samples, the reads were assembled de novo using a metagenomic assembler metaSPAdes v4.0.0 (Nurk et al. 2017), which is equivalent to running SPAdes with the “--meta” flag. No prior quality-based read trimming or filtering was done before the assembly. The assemblies were run on the High-Performance Computing cluster of the University of Greifswald. For the sample Lindbladia tubulina sc2833, which had the highest input DNA amount and the largest sequencing data volume, this assembly method failed with an out-of-memory error. Therefore, the assembly for this sample was made using the standard assembly mode of SPAdes v4.0.0 (Prjibelski et al. 2020).

For taxonomic profiling of contigs, each resulting assembly was searched against a locally installed NCBI nucleotide database (“core_nt”, updated 31 Jan 2026) using BLASTn v. 2.5.0+ in megablast mode (Camacho et al. 2009). BLAST output was generated in tabular format (“outfmt 6”) including subject TaxIDs, and these TaxIDs were linked to the NCBI taxonomy using the NCBI “taxdump” files.

To estimate per-contig sequencing depth, the original reads of each sample were mapped back to the corresponding assembly using strobealign v. 0.17.0 (Sahlin 2022). Alignment output (SAM, Sequence Alignment Map) was converted to sorted and indexed BAM (Binary Alignment Map) files using samtools v. 1.12 (Danecek et al. 2021). Taxonomic assignments, GC content, contig length, and read depth were integrated and visualized with BlobToolKit v. 4.5.0 (Challis et al. 2020). For each sample, BlobTools2 was used to create a project from the assembly FASTA, BLAST hits, and BAM coverage file, and summary snail and blob plots were produced to assess the distribution of contigs and read support across inferred taxonomic groups at different taxonomic levels.

For every contig where the BLAST search over a local core_nt database found a match to an organism belonging to Eumycetozoa, we additionally conducted a thorough online BLAST search over the NCBI GenBank Nucleotide database (see Suppl. material 1: Sheet B for a list of the best 10 matches with the smallest e-values and Sheet C for a list of genes found).

For all seven samples, we obtained a substantial volume of raw genetic data. The total sequence length was relatively uniform, ranging from 1134 to 2354 Mbp (Table 1). However, the total length of the assembled contigs varied by two orders of magnitude, ranging from 1.7 Mbp in L. epidendrum sc25089 to 162.0 Mbp in T. scabra sc11426. Compared to the genome size of Badhamia polycephala (Miller et al. 2021), this corresponds to an estimated genome coverage between 0.7% and 65%.

Visualization of the assemblies via blob plots (https://doi.org/10.5281/zenodo.20467237) demonstrates significant differences in the data structure between samples. In each specimen, there were a number of contigs (1–19) unambiguously classified as belonging to Eumycetozoa based on BLAST hits (Fig. 1, Tables 1, 2). However, most sequences within the resulting assemblies belonged to the no-hit fraction—sequences with no matches in the NCBI database. We hypothesize (see Discussion) that a substantial proportion of this fraction is composed of myxomycete genetic material, which remains unidentified because genomic data for myxomycetes are still extremely sparse (Takano et al. 2001; Schaap et al. 2016; Kang et al. 2017). In the sample Lindbladia tubulina sc2833 (Fig. 1A, C), the no-hit fraction accounted for 82.3 Mbp and had a unimodal GC content distribution with highest contig density within a 35% to 45% range and a predominant base coverage of 20–40×, matching those of the contigs identified as fragments of myxomycete nuclear genome. In Lycogala epidendrum UARK25089 (https://doi.org/10.5281/zenodo.20467237), the no-hit fraction was smaller (1.4 Mbp) but maintained a characteristic distribution profile, allowing its identification as a residual signal of the nuclear genome. In Trichia scabra sc11426 (Fig. 1B, D), the no-hit contig fraction had a bimodal GC content distribution, with the second peak at around 68%, suggesting heavy contamination by prokaryotic DNA. In samples with low input DNA amount and high fragmentation, such as Cribraria mirabilis sc10428 and Fuligo cinerea sc11445 (https://doi.org/10.5281/zenodo.20467237), the data were represented by a sparse cloud of short contigs with a broad unimodal GC content distribution and the base coverage peaking around one.

In the combined no-hit & Eumycetozoa fraction of all samples except L. tubulina sc2833, the N50 values vary between 207–266 bp (Table 1, Fig. 2) and the statistical mode of read length distribution varies between 56 and 237 bp. This reflects the dominance of short unresolved fragments in the assembly graph that likely resulted from dead ends, short sequence repeats, rare variants, abundant contaminants, low-coverage sequences or sequencing errors. These short and low-coverage contigs are not filtered out from the results by metaSPAdes. The sample L. tubulina sc2833 stands out, with the N50 reaching 4.19 Kbp, which is likely explained by the higher input DNA amount and the different assembly mode applied to it (standard SPAdes mode – see Methods).

Within the examined set of samples, specimen age did not show a statistically significant linear correlation with genomic data quality. The two specimens yielding the highest amounts of DNA (L. tubulina sc2833 and T. scabra sc11426) span widely different ages (33 and 64 years, respectively) and are not the youngest among the studied collections (Table 2). Similarly, no significant correlation was detected between specimen age and DNA yield, read length, or the number of recovered Eumycetozoa contigs across the WGA-treated subset (Fig. 3). While DNA degradation over time is indisputable, these observations suggest that within a timeframe of several decades, initial preservation conditions and sample-specific properties may outweigh chronological age as a limiting factor.

Although the detection of Eumycetozoa within the resulting data remains challenging due to the lack of reference genomes, we succeeded in identifying anchor contigs for all seven specimens that correspond to their preliminary morphological identification (Suppl. material 2). Despite the paucity of these fragments, their presence is critical for specimen barcoding—the primary objective of genome skimming for old herbarium collections. It is important to emphasize that for most specimens the identified Eumycetozoa markers are localized within the main no-hit cloud (Fig. 1; https://doi.org/10.5281/zenodo.20467237). The Sheet B in the Suppl. material 1 lists the ten best BLAST hits for each contig and the detected genes.

Nuclear small subunit ribosomal gene (nucSSU) sequences were identified for all specimens. Their length varied from 157 to 7755 bp, and base coverage ranged from 1 to 127 copies (Table 2). For six out of seven specimens, the detected fragments of nucSSU at least partially covered the ca. 600 bp fragment (Suppl. material 3) conventionally used for DNA barcoding in myxomycetes (Leontyev and Schnittler 2017).

The total number of recognized genes reached 28 (Suppl. material 1: Sheet C). Among them, the eukaryotic elongation factor 1-alpha (EF1a) gene was detected in four specimens, with the contig length exceeding 10 kbp in L. tubulina sc2833. For three specimens, contigs containing the mitochondrial cytochrome c oxidase subunit I (COI) gene were obtained, varying in length (up to 41 kb) and base coverage (14–1506). Sequences including the mitochondrial small subunit ribosomal gene (mtSSU) were recovered for only two specimens, which also differed in length (up to 2310 bp) and base coverage (13–330). Other identified genes like actin, alpha-tubulin, and several tRNAs have not been studied for bright-spored myxomycetes; consequently, the matches found did not allow for species-level identification, as the closest analogues were dark-spored myxomycetes and dictyostelids (Suppl. material 1: Sheet B).

Taxonomic interpretation of non-target contigs revealed a complex metagenomic background, the composition of which varies depending on specimen age, substrate conditions, and archival history (https://doi.org/10.5281/zenodo.20467237). The anthropogenic fraction, primarily represented by human DNA, was recorded in all examined specimens. The proportion of this fraction correlates with the total DNA yield: in libraries with low target material content (C. mirabilis sc10428, L. epidendrum sc25089), anthropogenic noise becomes dominant. In close association with human DNA, a skin microflora complex was identified, including the bacterium Cutibacterium acnes, representatives of the genus Staphylococcus, and lipophilic yeasts of the genus Malassezia. Domestic cattle DNA (Bos taurus) found in several samples likely also belongs to the anthropogenic fraction.

The second metagenomic fraction is conditionally natural, reflecting the substrate inhabitants and ecosystem relationships of the myxomycetes. Arthropod DNA was detected selectively; the most pronounced signal was recorded in specimen L. tubulina sc12075, where oribatid mites of the genus Oppiella, common inhabitants of decaying wood (Maraun and Scheu 2000), were identified. The bacterial component of this fraction is represented by organisms characteristic of forest ecosystems: the cyanobacterium Nostoc and the betaproteobacterium Variovorax (associates of rotting wood) were noted in specimen T. scabra sc11426, while nitrogen-fixing soil bacteria of the genus Bradyrhizobium dominate in L. tubulina sc2833. In the oldest specimen of the series, F. cinerea sc11445, a predominance of putrefactive bacteria of the genus Bacillus was observed. The archaeal component is limited to nitrifying soil forms (Nitrososphaera). The viral fraction is represented by fragments of bacteriophages.