Limitations of Common Molecular Markers in Fungal Biodiversity Analysis and the Benefits of Their Synergistic Use

TL;DR

Using only ITS markers for studying fungal diversity can be inaccurate, but combining multiple markers improves accuracy and provides a better picture of fungal communities.

Contribution

The study demonstrates that combining ITS1, ITS2, and protein-coding gene markers improves fungal diversity and abundance estimates compared to single-marker approaches.

Findings

ITS2 detected 75.3% of species, but combining ITS1, ITS2, and elongation factor 1-alpha gene detected 92.2%.

Multi-marker approaches correct biases in relative sequence abundance estimates.

Two markers failed to produce sufficient data for reliable diversity estimates.

Abstract

High‐throughput sequencing of the Internal Transcribed Spacer (ITS) regions is the primary method for estimating fungal diversity from environmental DNA. However, reliance solely on ITS markers is complicated by its high variability in sequence length and the presence of multiple variants within a single genome, which can bias diversity estimates. This study compares the performance of the ITS1 and ITS2 molecular markers with five alternative markers targeting protein‐coding genes using a defined mock community of 413 fungal species from 32 orders. We show that none of the markers provided a reliable estimate of fungal diversity and two failed to produce sufficient data. While the ITS2 marker showed the highest proportion of individual species detections at 75.3%, a multi‐marker approach significantly improved this result. Combining the top three markers (ITS1, ITS2 and elongation factor 1‐alpha gene) increased this number to 92.2%. This approach also provided a crucial corrective for relative sequence abundance (RSA) estimates, as combining markers helps to offset the individual biases of each of them resulting in a more balanced picture of the community composition. Although the ITS region remains a cornerstone in fungal metabarcoding studies, we recommend combining both ITS1 and ITS2 markers with at least one alternative DNA marker to achieve more accurate species richness and RSA estimates. The broader application of these alternative markers will depend on the continued development of comprehensive, curated reference databases.

Figures

(A) Bar plot showing the proportion of degenerate bases in primer sequences used for PCR amplification of the seven tested fungal barcode markers. (B) Bar plot summarising the bioinformatic processing outcomes of Illumina MiSeq sequence data from the fungal mock community for the seven tested markers. Categories of filtered and retained reads include: Low quality—reads with a mean PHRED score < 30; Putative chimeras—reads identified as chimeric by the UCHIME algorithm; contamination—reads assigned to species not included in the mock community; non‐target—reads not matching the expected marker region or assigned to non‐fungal taxa; low similarity—reads that matched mock community species but with < 98.5% sequence similarity; retained—reads that passed all filters and were retained in the final dataset.

(A) Bar plot showing the number of species recovered by each individual marker. (B) Bar plot showing cumulative species recovery for the best‐performing combinations of one to five markers. The combinations shown are: one (ITS2), two (ITS2 + ef1‐α), three (ITS2 + ef1‐α + ITS1), four (ITS2 + ef1‐α + ITS1 + rpb2b), and five (all tested markers). (C) Matrix layout plot visualising the overlap in species detections. The vertical bars show the number of species detected by specific marker combinations (intersections). Filled circles below the vertical bars specify which markers are part of each intersection.

(A) Bar plot showing the increase in species recovery achieved by combining ITS2 with each of the other four tested DNA markers: ITS1, mtLSU, rpb2b and ef1‐α. (B) Venn diagram showing the contribution of four individual markers, ITS1, mtLSU, rpb2b and ef1‐α, and their combinations to additional species detections when used in combination with ITS2.

(A) Bar plot showing the number and relative proportion of the input specimens by fungal order and division in the mock community. (B) Bar plots showing the number of amplicon sequences grouped by fungal order and division in the mock community and the difference in their relative proportions compared to the original input specimens for all tested markers. Order abbreviations: Ag—Agaricales, Bo—Boletales, Ca—Capnodiales, Ch—Chaetothyriales, Di—Diaporthales, Eu—Eurotiales, Gl—Glomerellales, He—Helotiales, Hy—Hypocreales, Le—Leotiomycetes incertae sedis, Mi—Microascales, Mu—Mucorales, On—Onygenales, Op—Ophiostomatales, Pe—Pezizales, Pl—Pleosporales, Po—Polyporales, Ru—Russulales, So—Sordariales, Xy—Xylariales.

Tables

TABLE 1. Percentage of species recovered from the fungal mock community by each barcode marker. Colour coding: Green—complete species recovery; yellow—partial recovery; red—no recovery. ‘n’ indicates the number of species evaluated within each fungal order. Only orders with at least three species (n ≥ 3) are shown.

Order	ITS1	ITS2	mtLSU	ef1‐α	rpb2b	Combined
Ascomycota
Eurotiales (n = 4)	75%	75%	75%	0%	75%	100%
Helotiales (n = 4)	75%	75%	100%	75%	0%	100%
Hypocreales (n = 20)	70%	50%	65%	85%	80%	100%
Microascales (n = 7)	57%	71%	43%	43%	29%	86%
Onygenales (n = 5)	60%	100%	80%	0%	40%	100%
Pezizales (n = 3)	100%	100%	0%	67%	67%	100%
Pleosporales (n = 13)	54%	77%	62%	77%	23%	92%
Xylariales (n = 9)	56%	89%	33%	78%	89%	100%
Basidiomycota
Agaricales (n = 179)	74%	79%	21%	44%	30%	93%
Boletales (n = 17)	41%	53%	59%	71%	35%	94%
Polyporales (n = 33)	76%	73%	24%	82%	58%	100%
Russulales (n = 53)	72%	72%	13%	25%	72%	94%
Mucoromycota
Mucorales (n = 13)	62%	69%	38%	77%	8%	85%

TABLE 2. The 20 most abundant species detected individually by each of the five tested barcode markers, along with their relative sequence abundance (RSA) and the corresponding number of fungal specimens (SP) in the mock community.

No.	ITS1	RSA	SP	ITS2	RSA	SP	mtLSU	RSA	SP	ef1‐α	RSA	SP	rpb2b	RSA	SP
1	Entoloma lividum	74.86%	1	Heterobasidion annosum	16.07%	3	Clonostachys rosea	19.12%	1	Akanthomyces muscarius	11.74%	5	Heterobasidion annosum	10.57%	3
2	Entoloma sericatum	18.84%	1	Epichloe typhina	6.56%	3	Melanoleuca humilis	8.96%	4	Russula helodes	8.26%	2	Melanoleuca monticola	7.19%	3
3	Heterobasidion annosum	1.48%	3	Heterobasidion parviporum	5.74%	3	Flavocillium bifurcatum	7.00%	1	Amanita citrina	6.09%	1	Heterobasidion parviporum	5.43%	3
4	Schizophyllum sp.	1.15%	2	Flavocillium bifurcatum	4.36%	1	Schizophyllum sp.	6.00%	2	Phoma sp.	5.02%	2	Russula nitida	4.44%	1
5	Heterobasidion parviporum	0.30%	3	Amanita submembranacea	3.61%	7	Melanoleuca grammopodia	5.50%	3	Flavocillium bifurcatum	4.97%	1	Russula amethystina	4.13%	2
6	Laetiporus sulphureus	0.27%	5	Laetiporus sulphureus	3.59%	5	Hymenopellis radicata	4.82%	5	Russula subfoetens	4.48%	17	Schizophyllum commune	3.70%	6
7	Cucurbitariaceae sp.	0.17%	1	Melanoleuca luteolosperma	2.90%	1	Melanoleuca monticola	3.00%	3	Didymella macrostoma	3.65%	1	Phoma sp.	2.68%	2
8	Amanita submembranacea	0.17%	7	Melanoleuca grammopodia	2.90%	3	Verticillium sp.	2.88%	1	Laetiporus sulphureus	3.37%	5	Russula subfoetens	2.53%	17
9	Schizophyllum commune	0.17%	6	Schizophyllum commune	2.87%	6	Mucidula mucida	2.65%	14	Simplicillium lamellicola	2.75%	1	Tricholomopsis flammula	2.33%	3
10	Flavocillium bifurcatum	0.16%	1	Kretzschmaria deusta	1.92%	1	Kuehneromyces mutabilis	2.61%	1	Russula sphagnicola	2.68%	3	Gamszarea sp.	2.26%	2
11	Ovicillium variecolor	0.11%	2	Melanoleuca malenconii	1.80%	3	Serpula himantioides	2.37%	1	Gamszarea sp.	2.60%	2	Hebeloma sordidum	2.19%	6
12	Russula subfoetens	0.10%	17	Melanoleuca monticola	1.75%	3	Isaria farinosa	1.92%	1	Russula maculata	2.48%	1	Akanthomyces muscarius	2.10%	5
13	Armillaria ostoyae	0.10%	4	Serpula himantioides	1.43%	1	Melanoleuca stridula	1.89%	2	Melanoleuca grammopodia	2.01%	3	Melanoleuca humilis	2.00%	4
14	Melanoleuca grammopodia	0.10%	3	Lycoperdon echinatum	1.36%	1	Gamszarea sp.	1.86%	2	Akanthomyces attenuatus	1.81%	1	Russula mutabilis	1.94%	1
15	Melanoleuca humilis	0.10%	4	Entoloma lividum	1.34%	1	Graphium fimbriisporum	1.75%	1	Russula foetens	1.80%	6	Ovicillium variecolor	1.89%	2
16	Melanoleuca monticola	0.07%	3	Russula subfoetens	1.26%	17	Tricholomopsis flammula	1.65%	3	Isaria farinosa	1.58%	1	Russula illota	1.81%	4
17	Psilocybe serbica	0.06%	6	Trametes sanguinea	1.22%	2	Phoma sp.	1.63%	2	Russula integra	1.28%	3	Schizophyllum radiatum	1.75%	2
18	Isaria farinosa	0.05%	1	Pyrenochaeta sp.	1.20%	1	Cordyceps confragosa	1.63%	1	Mycetinis alliaceus	1.25%	3	Heterobasidion abietinum	1.65%	3
19	Mucidula mucida	0.05%	14	Melanoleuca humilis	1.19%	4	Schizophyllum commune	1.60%	6	Schizophyllum sp.	1.11%	2	Agaricus altipes	1.57%	1
20	Phoma herbarum	0.04%	1	Kalmusia longispora	1.11%	1	Melanoleuca pseudopaedida	1.53%	1	Psilocybe cyanescens	1.02%	3	Russula sororia	1.55%	2

Full text

Introduction

1

The advent of high‐throughput sequencing (HTS) technologies has greatly enhanced our ability to study fungal communities in various environments (Baldrian et al. 2022, 2023). In recent decades, metabarcoding using HTS of PCR‐amplified short DNA regions from mixed DNA samples has become the most widely used technique in fungal ecology (Nilsson et al. 2019; Oliveira and Azevedo 2022; Niskanen et al. 2023). This technique allows detection of multiple fungal taxa in a single environmental DNA sample, yet results remain highly sensitive to methodological biases including barcode region selection, PCR primer design, PCR amplification and downstream bioinformatics (Tedersoo et al. 2022; Hakimzadeh et al. 2024).

The most universal and widely used barcode region for fungi is the multi‐copy, non‐coding Internal Transcribed Spacer (ITS) region of the nuclear ribosomal DNA (Větrovský et al. 2020; Kauserud 2023). This region is located between the small (18S) and large (28S) ribosomal subunit genes, and it includes the ITS1 and ITS2 subregions separated by the 5.8S gene. The high copy number of ITS often facilitates amplification from degraded or low‐biomass samples (Bellemain et al. 2010; Toju et al. 2012; Nguyen et al. 2015), while its rapid evolution supports species‐level resolution in many fungal groups (Schoch et al. 2012; Kõljalg et al. 2013).

Despite these advantages, the ITS region also has several limitations. ITS copy number and length vary across fungal taxa (Lofgren et al. 2019; Wilson et al. 2023), and intra‐genomic divergence can reach up to 11.8% within a single species (Cedeño‐Sanchez et al. 2024). In metabarcoding, stochastic PCR amplification of divergent sequence variants may lead to ambiguous species identification or even omission of species with unusually long or short ITS sequences (Bellemain et al. 2010; Tedersoo et al. 2015; Badotti et al. 2017; Castaño et al. 2020). These limitations indicate the necessity for alternative fungal barcodes that would offer more consistent amplification and reliable species‐level resolution.

There is a long tradition of use of protein coding DNA regions for fungal phylogeny (Matheny et al. 2007; Raja et al. 2011; Lücking et al. 2020) which may serve as alternative fungal barcodes. These markers are typically present in only one or a few copies per genome, lack the length variation characteristic of the ITS region and possess higher discriminatory capacity, particularly for resolving closely related species within complex genera (Stielow et al. 2015; Xu 2016; Swenie et al. 2023). Their use is significantly limited by the absence of universal primers, as testing on single‐template DNA consistently shows low taxonomic coverage and amplification success (Schoch et al. 2012; Stielow et al. 2015). In addition, all protein coding primers exhibit a high degree of degeneracy, which increases the risk of amplification of non‐target DNA templates, which can bias community composition and reduce taxonomic resolution (Tedersoo et al. 2015, 2018; Tedersoo and Lindahl 2016). Only a few studies tested primers for alternative barcode markers in metabarcoding. The most promising application is for taxa where ITS has low discriminatory power. This is the case for Fusarium, where elongation factor 1‐alpha gene (ef1‐α), used with genus‐specific primers, is a viable marker for metabarcoding (Boutigny et al. 2019). The same marker outperformed ITS1 in both taxonomic coverage and relative quantification when tested on a mock community comprising several clinical fungi from three fungal phyla (Weaver et al. 2024). ITS and the beta tubulin gene (TUB2) proved to be complementary in the analysis of endophytic environmental DNA (Theologidis et al. 2023). Finally, gene encoding the second largest subunit of RNA polymerase II (rpb2), studied using a single primer pair, yielded low taxonomic coverage in the study of (Rué et al. 2023), while it showed broad taxonomic coverage, comparable to outperforming ITS marker in others (Větrovský et al. 2016; Shapkin et al. 2025).

Nevertheless, these studies were either based on environmental DNA samples or mock community covering phylogenetically narrow groups of fungi, limiting conclusions about their general utility for fungal metabarcoding. Extensive study that compares the performance of multiple universal primers for alternative fungal markers in metabarcoding using a large and taxonomically diverse fungal mock community is still lacking.

The presence of individual taxa detected in fungal metabarcoding is often quantified using their relative sequence abundance (RSA), defined as the proportion of total amplicon sequences attributed to each taxon. However, the accuracy of taxon presence quantification based solely on amplicon sequencing results remains a subject of ongoing debate (Piñol et al. 2015; Deagle et al. 2019; van der Loos and Nijland 2021; Mächler et al. 2021). Biases in taxon abundance estimates derived from the ITS region have been demonstrated in several studies (Bakker 2018; Lamb et al. 2019; Castaño et al. 2020), whereas information on such biases for alternative fungal barcode markers remains very limited (Větrovský et al. 2016).

For the comparison of well performing molecular markers, we selected both ITS1 and ITS2 and five protein‐coding regions. The first marker is located within the mitochondrial large subunit ribosomal DNA locus (mtLSU). Similarly to nuclear rDNA genes, this is a multi‐copy gene but with lower mutation rates (Zeng et al. 2004; Lang 2014). Previous studies showed that the application of partial mtLSU as a marker for environmental metabarcoding has decent performance in the identification of arbuscular mycorrhizal fungi (Huang et al. 2014). The other two markers are within the rpb2 gene. The first rpb2 barcode (rpb2a) is located between the 6th and 7th conservative regions of the rpb2 gene (Matheny 2005). This marker has already proved its applicability for fungal metabarcoding (Větrovský et al. 2016). The second rpb2 marker (rpb2b) is in a shorter region between the 5th and 6th conserved partitions of rpb2 and was more suitable for the Illumina sequencing option (Reeb et al. 2004). Ef1‐α was included because it was introduced as a secondary fungal barcode region known for its ability to discriminate between species and for its high amplification success rates across a broad range of fungi (Rehner and Buckley 2005; Stielow et al. 2015). Lastly, we selected one marker from the minichromosome maintenance complex component 7 gene (mcm7). This marker is also known for its good amplification success rates and species discrimination across major groups of fungi (Schmitt et al. 2009; Raja et al. 2011; Schoch et al. 2012). However, this gene has not yet been tested in fungal metabarcoding prior to this study.

We aimed to evaluate and compare the performance of seven fungal DNA barcode markers, ITS1, ITS2, mtLSU, rpb2a, rpb2b, ef1‐α and mcm7, in metabarcoding. Using a large, taxonomically diverse mock community covering fungal phyla Ascomycota, Basidiomycota, Mucoromycota and Zoopagomycota, we assessed each marker's amplification success, species recovery, and accuracy in estimating relative taxon abundances. The overall goal was to identify optimal marker combinations for improved fungal diversity assessment. A central objective was to test whether highly degenerate primers, which were not described as universal, have potential for future use in metabarcoding studies. Based on the current knowledge, we expect that, despite previous experience with single‐template DNA, the primers targeting protein‐coding regions may perform more universally in a multi‐template setting. Furthermore, we hypothesize that these low‐copy markers may provide more accurate estimates of the proportional representation of taxa within the community compared to the multicopy ITS markers.

Materials and Methods

2

Mock Community Composition and Processing

2.1

The mock community consisted of an equimolar mixture of genomic DNA extracted from 676 fungal specimens obtained either from pure cultures (613) or field collections (63). All these fungi belonged to the group Eumycota: 103 Ascomycota strains, 554 Basidiomycota strains, 18 Mucoromycota strains and 1 Zoopagomycota strain, corresponding to 413 different species, 202 genera, 96 families, 32 orders and 11 classes. Altogether, 71 genera were represented by several species, and 125 species were represented by more than one specimen.

All fungal specimens were identified based on morphology and full‐length ITS region sequences. We received biological materials from the Culture Collection of Fungi, Department of Botany, Faculty of Science of Charles University in Prague; the Culture Collection of Basidiomycetes, Institute of Microbiology CAS; collections of the Laboratory of Fungal Genetics and Metabolism, Institute of Microbiology CAS; collections from the Department of Forest Protection and Wildlife Management at the Faculty of Forestry and Wood Technology, Mendel University in Brno; collections from the Department of Environmental Geology and Geochemistry, Institute of Geology CAS; and from the herbarium of the National Museum in Prague. The origin of each fungal strain is documented in Table S1.

Fungal DNA was obtained either from 7 to 14 days old cultures or from fresh fruiting bodies using UltraClean Microbial DNA isolation kit (MoBio Laboratories Inc., USA). DNA from each fungal specimen was quality‐checked by PCR to ensure that only amplifiable, high‐quality DNA was included in the mock community. DNA concentrations for each specimen were then quantified with a NanoDrop 1000 Spectrophotometer (NanoDrop, USA) and pooled in equimolar proportions. Primer pairs used for amplification of each DNA marker are listed in Table S2. PCRs were performed in triplicates independently for each sample and marker combination. Each 25‐μL reaction mixture contained 12.7 μL H_2_O, 2.5 μL 10× Buffer for PerfectTaq Plus DNA Polymerase, 3 μL 10 mg/mL BSA, 0.5 μL PCR Nucleotide Mix (10 mM), 2 μL forward primer (final concentration, 10 pmol/μL), 2 μL reverse primer (final concentration, 10 pmol/μL), 2 μL template DNA and 0.3 μL PerfectTaq Plus DNA Polymerase (Final Concentration, 1.5 U/μL) (5 PRIME, Germany). A similar PCR program was used for all DNA markers: initial denaturation (94°C, 5 min); denaturation (94°C, 45 s), annealing (see Table S2, 30 s), elongation (72°C, 30 s)—repeated 37 times; final elongation (72°C, 10 min). Markers rpb2a and rpb2b required additional Mg^2+^ to maximise the yield. The quality and yield of PCRs were electrophoretically verified on agarose gel. PCR triplicates were pooled together and purified using ZR‐96 DNA Clean & Concentrator kit (Zymo Research Corporation, USA). Post PCR samples were quantified by Quant–iT PicoGreen dsDNA kit (Invitrogen, USA) using Bio‐Rad CFX96 Touch Real‐Time PCR Detection System (Bio‐Rad, USA) then equally pooled based on the estimated number of amplicons. Ligation steps followed the Low Sample Protocol of Illumina TruSeq DNA PCR‐Free Sample Preparation Guide (Illumina Inc., USA). Final libraries were cleaned‐up with AMPure XP Beads (Beckman Coulter Inc., USA) and quantified using the KAPA Library Quantification Kit Illumina platforms (Kapa Biosystems Inc., USA) following the manufacturer's instructions. Sequencing was performed in triplicates on Illumina MiSeq platform using the MiSeq Reagent Kit v3 (600 cycle) for 2 × 300‐bp reads (GeneTiCA s.r.o., CZE) in two independent sequencing runs.

Bioinformatics Analysis

2.2

We analysed Illumina MiSeq amplicon sequence data using the SEED 2 pipeline (Větrovský et al. 2018) which integrates internal functions and external software. Forward (R1) and reverse (R2) reads were joined using the utility Fastq‐join (Aronesty 2013) with a minimum overlapping length of 40 bp and less than 15% of mismatches in the overlap region. There was no overlap between forward and reverse reads obtained from the longer regions of rpb2a and mcm7 markers. For these markers we analysed forward reads only. Reads with a PHRED quality score (Q‐score) mean below 30, ambiguous bases or a mismatch in the tag were excluded. Adaptors and tags were trimmed. ITS1 and ITS2 reads were truncated to contain only the highly variable region using the utility ITSx v. 1.0.11 (Bengtsson‐Palme et al. 2013). ITS1 and ITS2 reads shorter than 30 bp, and mtLSU, ef1‐α, rpb2a, rpb2b, mcm7 reads shorter than 200 bp were excluded. Putative chimeric reads were identified and excluded using de novo strategy in the UCHIME (Edgar et al. 2011) algorithm implemented in the USEARCH (Edgar 2010).

Reference Sequence Dataset Generation

2.3

For taxonomic assignment of amplicon sequence data resulted from sequencing of the mock community we used reference sequences of individual fungal specimens. The reference sequences were obtained by Illumina MiSeq sequencing of 64 small subsets of fungal specimens included in the mock community. To avoid difficulties with recognising sequences of closely related species, we prepared small subsets of 10–11 phylogenetically unrelated fungi preferably from different families. Sample preparation, sequencing and bioinformatics analysis were done the same way as described above. Reads were grouped into operational taxonomic units (OTUs) at 97% similarity threshold using the UPARSE clustering algorithm implemented in the USEARCH (Edgar 2010). We selected the most abundant sequence variant from each OTU and assigned it to the family as the reference sequence based on best BLASTn match in the NCBI database (for mtLSU, ef1‐α, rpb2a, rpb2b and mcm7 markers) and the UNITE database (for ITS1 and ITS2 markers). Fungal sequences from the NCBI were downloaded on 14.04.2024 using queries listed in Table S3. The UNITE database was downloaded as the general FASTA release version 10.0 (as of 04.04.2024). In subsets containing more than one species from a single fungal family (Russulaceae, Strophariaceae, Tricholomataceae), species identifications were also supported by sequence comparison in alignments and construction of phylogenetic trees using PhyML 3.0 (Guindon et al. 2010).

The resulting reference sequence collections included ITS1 sequences from 650 specimens, ITS2 from 646, ef1‐α from 409, rpb2b from 332, mtLSU from 212, rpb2a from 195 and mcm7 sequences from 106 specimens (Supplementary Sequence Data). To enable a fair comparison of markers, we ensured that all datasets contained an identical set of orders. Any order that was missing from even one marker's reference collection was removed from all of them. This standardisation resulted in the exclusion of 52 specimens across 12 orders (Table S4). Consequently, our final analyses were performed using a consistent set of 20 orders across all marker datasets, with our reference collections for the ITS1 and ITS2 regions together providing the most complete coverage for species identification.

Taxonomic Assignment of Reads From Mock Community

2.4

OTU clustering based on fixed similarity thresholds for species richness analysis is a widely used practice in environmental DNA studies. However, this method often leads to over‐splitting or over‐merging of sequences from individual taxa (Ryberg 2015; Tedersoo et al. 2022). This issue is particularly pronounced when comparing different groups of fungi (Bonin et al. 2023; Wilson et al. 2023), and different genomic regions (Shapkin et al. 2025).

To avoid clustering‐based biases, every unique sequence variant from the mock community was identified to species by comparing with the reference sequence collections using the BLASTn algorithm with thresholds set at ≥ 98.5% sequence similarity and ≥ 95% sequence coverage. Based on the sequence similarity within the reference dataset, this threshold is sufficient to distinguish all taxa present in the mock community. Higher‐rank classification was assigned according to UNITE (Abarenkov et al. 2024). For species included in the mock community, we accepted detections even if supported by a single annotated read.

Reads that did not meet the initial high‐confidence assignment were additionally compared against fungal sequences deposited in public repositories (NCBI or UNITE) then categorised using a broad, low‐confidence cutoff of 70% sequence similarity, which served to separate potentially relevant signals from unrelated noise. Reads that best matched mock community species were labelled as ‘low similarity’. Those that best matched species not included in the mock community were classified as ‘contamination’. Finally, any read failing to meet this cutoff was labelled as ‘non‐target’.

Final Statistical Analyses and Data Visualisation

2.5

All data were imported into R software environment version 4.2.0 (R Core Team 2020) for statistical computing and graphics. Data manipulation and preparation for analysis was done using the package dplyr version 1.1.4 (Wickham et al. 2023). To test the correlation between the numbers of fungal specimens and amplicon sequences we performed Spearman's rank correlation using the cor_mat() function from the package rstatix version 0.7.2 (Kassambara 2023). All barplots were prepared using the package ggplot2 version 3.4.4 (Wickham 2016). For matrix layout visualisation of detected species, we used the package UpSetR version 1.4.0 (Conway et al. 2017).

Results

3

Quality of Sequencing Data and Amplification Success

3.1

Illumina MiSeq sequencing of the studied fungal mock community resulted in the following numbers of assembled paired‐end reads: 458,269 for ITS1, 153,396 for ITS2, 356,491 for mtLSU, 580,047 for ef1‐α, and 542,610 for rpb2b. Reverse and forward reads of rpb2a and mcm7 did not overlap; therefore, we used only forward reads for these two markers. In total, 128,415 for the first one and 859,445 for the second were obtained.

Read filtering based on the average Q‐score of all bases removed 18.8% of rpb2a and 12.6% of mcm7 forward reads but only 0.05%–0.3% of assembled reads for ITS1, ITS2, mtLSU, ef1‐α and rpb2b markers. Truncation of ITS1 and ITS2 regions using ITSx (Bengtsson‐Palme et al. 2013) removed less than 0.1% of reads in both datasets. A 200 bp minimum length filter removed 42.2%–67.1% of reads from the low‐copy markers (rpb2a, rpb2b, ef1‐α and mcm7), while 30 and 200‐bp filters removed only 0%–2.0% of reads from the multi‐copy markers (ITS1, ITS2 and mtLSU). The BLASTn search on the discarded short reads against all available sequences in the NCBI database revealed that the majority (e.g., ~97% for rpb2b) were not specific to any organism. The UCHIME's de novo algorithm identified as chimeric, and removed 1.95% of ITS1 reads, 12.3% of ITS2, 16.3% of mtLSU and 6.0% of ef1‐α but only 0.34%–0.59% of rpb2a, rpb2b and mcm7 reads.

The BLASTn search of filtered sequence data for each marker against collections of reference sequences and fungal sequences available in public databases revealed that 0.84% of ITS reads, 7.6% of ITS2, 0.2% of mtLSU, 5.5% of ef1‐α, 12.4% of rpb2a, 21.1% of rpb2b, and 15.8% of mcm7 reads received no hit to any sequence. These reads potentially belonged to non‐target sequences or non‐fungal organisms. The average proportion of reads that best matched sequences for species not included in the mock community was below 0.15% across all markers.

Despite available reference sequence data for fungi included in the mock community, a substantial fraction of reads remained taxonomically unassigned because they did not meet the assignment criteria (minimum similarity threshold of 98.5%). This reads with low‐similarity to reference sequences represented 19.7% of ITS1, 29.5% of ITS2, 20.4% of mtLSU, 12.5% of ef1‐α, 11.4% of rpb2a, 15.2% of rpb2b, and 4.3% of mcm7 reads. Their similarity was too low for the taxonomic assignment, but they were similar to reference data more than 70%, which is the threshold for filtering out non‐target reads. Most of these reads were within 97% similarity to the reference sequences. This suggests that they likely resulted from technical errors (e.g., sequencing, PCR or chimeric sequences missed by UCHIME) or biological variation (e.g., pseudogenes or paralogs) which would be masked by centroid sequences in traditional OTU clusters.

Numbers of reads sorted by each filtering step are listed in Table S5 and proportionally visualised in Figure 1B. The rpb2a and mcm7 datasets retained very low numbers of target sequences, representing less than 10% of the total species diversity from the mock community. Therefore, we excluded these two markers from further evaluation.

(A) Bar plot showing the proportion of degenerate bases in primer sequences used for PCR amplification of the seven tested fungal barcode markers. (B) Bar plot summarising the bioinformatic processing outcomes of Illumina MiSeq sequence data from the fungal mock community for the seven tested markers. Categories of filtered and retained reads include: Low quality—reads with a mean PHRED score < 30; Putative chimeras—reads identified as chimeric by the UCHIME algorithm; contamination—reads assigned to species not included in the mock community; non‐target—reads not matching the expected marker region or assigned to non‐fungal taxa; low similarity—reads that matched mock community species but with < 98.5% sequence similarity; retained—reads that passed all filters and were retained in the final dataset.

Taxonomic Coverage

3.2

After filtering out specimens which lacked references, the mock community comprised 372 species from 20 orders across three fungal divisions (Table S6). Among the tested markers, ITS2 achieved the highest species recovery with 75.3% of species detected (280 species), followed by ITS1 with 69.6% (259 species), ef1‐α with 51.1% (190 species), rpb2b with 43.0% (160 species), and mtLSU with 30.4% (113 species) (Figure 2A). Only 6.7% of species (25 species) were detected by all five markers, while 5.6% (21 species) were not detected by any marker. Several species were detected by only a single marker: 4.0% (15 species) were exclusive to ITS2, 3.2% (12 species) to ef1‐α, 2.1% (8 species) to ITS1, 1.1% (4 species) to rpb2b, and 0.8% (3 species) to mtLSU (Figure 2C).

(A) Bar plot showing the number of species recovered by each individual marker. (B) Bar plot showing cumulative species recovery for the best‐performing combinations of one to five markers. The combinations shown are: one (ITS2), two (ITS2 + ef1‐α), three (ITS2 + ef1‐α + ITS1), four (ITS2 + ef1‐α + ITS1 + rpb2b), and five (all tested markers). (C) Matrix layout plot visualising the overlap in species detections. The vertical bars show the number of species detected by specific marker combinations (intersections). Filled circles below the vertical bars specify which markers are part of each intersection.

Species recovery varied across individual markers and fungal orders (Table 1). The ITS1 and ITS2 markers detected the highest numbers of species within the orders Microascales, Pezizales and Agaricales. The ef1‐α marker achieved the highest species recovery in Boletales, Polyporales and Mucorales, but failed to detect any species from Eurotiales and Onygenales. The mtLSU marker was the only one to recover all species from Helotiales, but did not detect any species from Pezizales. The rpb2b marker recovered the highest numbers of species from Xylariales and Russulales, the same as ITS2, but the fewest number of species from Mucorales and failed to detect any species from Helotiales.

Among the five markers, ITS2 recovered the highest proportion of species and was the only marker that successfully amplified at least some species from every fungal order represented in the mock community. Given its broad taxonomic coverage and strong performance, we selected the 75.3% species recovery achieved by ITS2 as the baseline for evaluating the additive contribution of the other markers to total species recovery in combined marker approaches.

Combining ITS2 with ef1‐α increased species recovery by 11.8% to 87.1% (324 species), with ITS1 by 10.5% to 85.8% (319), with rpb2b by 9.4% to 84.7% (315), and with mtLSU by 6.7% to 82.0% (305; Figure 3A). Although many species detections overlapped among these markers, combining more markers together further improved species recovery (Figure 3B). The most effective three‐marker combination was ITS2 + ITS1 + ef1‐α, which yielded 92.2% species recovery (343 species). Adding rpb2b as a fourth marker further increased recovery to 93.8% (349 species). The highest recovery of 94.6% (352 species) was achieved by combining all five markers (Figure 2B).

(A) Bar plot showing the increase in species recovery achieved by combining ITS2 with each of the other four tested DNA markers: ITS1, mtLSU, rpb2b and ef1‐α. (B) Venn diagram showing the contribution of four individual markers, ITS1, mtLSU, rpb2b and ef1‐α, and their combinations to additional species detections when used in combination with ITS2.

Notably, 5.6% of the species in the mock community were not detected by any of the five markers during metabarcoding, although they were successfully amplified by at least some of them during the preparation of reference amplicon subsets by Illumina sequencing.

Relative Sequence Abundance

3.3

After the filtering and taxonomy assignment steps described above, the following numbers of amplicon sequences, belonging to species included in the mock community, remained in datasets for each marker: 77.3% for ITS1 (354,396 sequences), 48.3% for ITS2 (75,286), 61.0% for mtLSU (217,514), 20.9% for rpb2b (113,118), and 11.9% for ef1‐α (68,719) (Table S7).

To evaluate how well each marker recovered the structure of the community, we compared ranks of RSA per fungal order to the ranks of the specimen counts per order. The specimen counts served as a proxy for the true relative DNA composition of the community due to equal input of DNA per specimen (Figure 4A).

(A) Bar plot showing the number and relative proportion of the input specimens by fungal order and division in the mock community. (B) Bar plots showing the number of amplicon sequences grouped by fungal order and division in the mock community and the difference in their relative proportions compared to the original input specimens for all tested markers. Order abbreviations: Ag—Agaricales, Bo—Boletales, Ca—Capnodiales, Ch—Chaetothyriales, Di—Diaporthales, Eu—Eurotiales, Gl—Glomerellales, He—Helotiales, Hy—Hypocreales, Le—Leotiomycetes incertae sedis, Mi—Microascales, Mu—Mucorales, On—Onygenales, Op—Ophiostomatales, Pe—Pezizales, Pl—Pleosporales, Po—Polyporales, Ru—Russulales, So—Sordariales, Xy—Xylariales.

Among the tested markers, ITS1 showed the strongest correlation with expected structure of the mock community (Spearman's r = 0.94, p < 0.01). However, this high correlation coincided with an overrepresentation of Agaricales, which represented 96% of ITS1 sequences, and all other orders were underrepresented. The ITS2 and ef1‐α markers both produced RSA profiles closely aligned with the mock community structure (r = 0.91 and r = 0.81, respectively; p < 0.01). In contrast, rpb2b and mtLSU showed weaker correlations (r = 0.74 and r = 0.71, respectively; p < 0.01). The lower correlation observed for mtLSU was primarily due to poor representation of species from Polyporales and Russulales, while rpb2b showed strong overrepresentation of Russulales. Despite these marker‐specific differences, some RSA deviations were consistent across datasets. Hypocreales and Pleosporales were overrepresented and Agaricales strongly underrepresented by all markers except ITS1 (Figure 4B).

At the species level, RSA was highly uneven and showed only a weak correlation with the number of fungal specimens per species (Table S8). To explore this further, we examined the 20 most abundant species recovered by each marker (Table 2). In the ITS1 dataset, most sequences were concentrated in just two species: Entoloma lividum (RSA 74.9%) and Entoloma sericatum (RSA 18.8%), each of which was represented by a single fungal specimen. These two species were represented by only 1.7% of sequences in the ITS2 dataset and were absent from all other marker datasets. The most abundant species in the ef1‐α dataset was Akanthomyces muscarius (RSA 11.7%), while mtLSU was dominated by Clonostachys rosea (RSA 19.1%). In the ITS2 and rpb2b datasets, the most abundant species was Heterobasidion annosum, with RSA values of 16.1% and 10.6%, respectively.

The overlap in the 20 most abundant species detected by each marker varied. There were 11 shared species between the ITS1 and ITS2 datasets, 5 species shared among the multi‐copy markers (mtLSU, ITS1 and ITS2), and 4 species shared between the low‐copy markers ef1‐α and rpb2b. No single species was among the top 20 most abundant across all five markers.

Less than 2.0% of species detected by the mtLSU (one species) and rpb2b (three species) markers were singletons (species represented by a single sequence). In contrast, the ITS1 dataset contained 13.1% (34 species), ef1‐α—5.2% (10 species), and ITS2—4.6% (13 species) singletons. This indicates that sequencing depth for mtLSU and rpb2b was likely sufficient, while deeper sequencing may improve species recovery for ITS1, ef1‐α and ITS2.

Discussion

4

Taxonomic Coverage Across Markers

4.1

Successful amplification of a molecular marker is a fundamental requirement for any set of PCR primers in fungal metabarcoding. Primers must not only amplify across a broad taxonomic range but also do so consistently and proportionally, without favouring certain lineages or overlooking others (Ihrmark et al. 2012). Low‐copy gene markers with greater resolution power compared to ITS markers could provide an additional perspective on sequenced community composition. However, their limitation lies in the trade‐off between taxonomic universality and amplification specificity, dictated by the number of degenerate positions in primer design.

Among the tested markers, mtLSU was amplified using the most specific primer pair (Zeng et al. 2004). Originally designed for detection of common airborne fungi, these non‐degenerate primers showed decent taxonomic coverage across most orders included in the mock community. However, the efficiency of their amplification was considerably biased towards the Ascomycota group, which is consistent with their original design because they were developed using an alignment of almost exclusively sequences from that group. This taxonomic bias limits their utility in broader community profiling. To enable wider application of the mtLSU marker in fungal metabarcoding, future efforts should focus on designing less specific primers that achieve more balanced amplification across major fungal lineages.

Primer pairs selected for amplification of ITS1 and ITS2 markers represent the established standards in fungal barcoding, with exception of ITS1FI2 (Bellemain et al. 2010; Větrovský et al. 2020). The ITS1FI2 primer was designed specifically for deep sequencing on the Illumina platform (Schmidt et al. 2013) but has not been widely adopted in fungal metabarcoding. These markers showed the highest taxonomic coverage, and ITS2 was the only one to amplify at least one species from each order. This superior performance of the ITS2 subregion aligns with a recent study by Winand et al. (2025), who found that ITS2 generally provided higher precision with comparable species recall than ITS1 in an extensive comparison using 37 mock communities. However, the debate is not fully settled, as other studies indicate that ITS1 may be more suitable for studying fungal communities in particular environments (Mbareche et al. 2020; Rué et al. 2023).

For amplification of ef1‐α and rpb2b markers we used primer pairs which were not originally designed as matching forward‐reverse combinations but were paired in our study to generate amplicons of suitable length (~400 bp) for Illumina sequencing. Both primer sets contain a substantial number of degenerate positions, intended to broaden their taxonomic range across Ascomycota and Basidiomycota (Liu et al. 1999; Rehner and Buckley 2005). Both markers lacked taxonomic coverage in Ascomycota and the rpb2b also in Mucoromycota (Table 1), but overall amplification efficiency across orders included in the mock community was closer to those of ITS markers compared to the other alternative markers tested in this study.

Primer pairs used for amplification of rpb2a and mcm7 markers target longer regions and they contain most degenerate positions among the tested primers (Liu et al. 1999; Schmitt et al. 2009; Figure 1A). In our study, we failed to obtain a sufficient amount of sequence data for evaluation of these markers using these primers. The same was reported in another study conducted on the Illumina MiSeq platform using the rpb2a marker (Rué et al. 2023). In contrast, in the metabarcoding study conducted on a sequencing machine which can produce longer reads, this marker amplified using the same primer pair showed more accurate and less biased amplification compared to the full ITS region (Větrovský et al. 2016). This suggests that the performance of these markers was rather limited by the length of the amplicons they generate (typically exceeding 600 bp), which poses a challenge for Illumina short‐read platforms. Therefore, in accordance with Větrovský et al. (2016), the rpb2a remains a promising marker to use in metabarcoding studies where longer read sequencing is possible. For the mcm7, this study represents the first documented attempt to use this marker in fungal metabarcoding and requires more research to confirm the utility.

The markers ITS1, ITS2 and mtLSU yielded a significantly higher proportion of quality‐filtered reads and target sequences, compared to markers ef1‐α, rpb2a, rpb2b and mcm7. The likely explanation is the difference in primer degeneracy, which was lowest for ITS1 and mtLSU and highest for the protein‐coding genes, especially mcm7 (Figure 1). PCR with degenerated primers with low amplification success under single‐template DNA conditions is capable of much broader taxonomic coverage, or even significant amplification of non‐target genes, in environments with multi‐template environmental DNA. This behaviour of degenerate primers was reported in other studies using ef1‐α and rpb2 genes in fungal metabarcoding (Větrovský et al. 2016; Boutigny et al. 2019; Shapkin et al. 2025). The main drawback of degenerate primers for sequencing of protein‐coding genes in metabarcoding on the mock community was the significant proportion of amplified non‐target sequences. In our study, we were able to separate non‐target sequences from the sequences of species included in the mock community and still recover fungal community structure at sufficient sequencing depth, because of reference sequence data used. However, in environmental DNA metabarcoding, this issue can substantially complicate bioinformatics analysis and potentially distort diversity estimates (Piñol et al. 2015; Tedersoo and Lindahl 2016).

An important challenge remains the selection of an appropriate OTU clustering threshold for species diversity assessment in fungal metabarcoding. In this study, a 97% clustering threshold was utilised during the generation of reference sequence datasets from small subsets of phylogenetically unrelated fungi. This approach worked well because the strains in each subset were initially selected to cover broad taxonomic ranks with few closely related taxa. However, the impact of OTU clustering threshold selection must be considered carefully in environmental studies, where the occurrence of closely related species is significantly more likely (Shapkin et al. 2025; Noffsinger et al. 2025).

Furthermore, it is important to note that our conclusions regarding marker universality are primarily driven by the Basidiomycota‐rich composition of our mock community (554 of 676 strains). Consequently, conclusions regarding underrepresented groups, such as Mucoromycota and Zoopagomycota, should be considered tentative given the limited number of strains included in the analysis.

Multi‐Marker Approach for Better Taxonomic Coverage

4.2

Fungal metabarcoding offers standardised protocols for PCR and HTS of the ITS region and pre‐set data processing pipelines for sequence data analysis (Tedersoo et al. 2022; Hakimzadeh et al. 2024). The main limitation is currently considered to be the availability and quality of reference sequences essential for the accurate annotation of amplicon sequences to known species (Hleap et al. 2021; Bik 2021). Our results showed that even the most complete collections of reference sequences for the ITS region do not guarantee recovery of the full species diversity present in a fungal community using a single primer pair/marker approach. ITS2 recovered the highest number of species, outperforming ITS1 and all alternative markers but reached only 75.3% recovery of species included in the mock community (Figure 2). The loss of species recovery when using only one subregion of the full ITS region was described in multiple previous studies (Bazzicalupo et al. 2013; Wang et al. 2015; Yang et al. 2018; Shapkin et al. 2025). Combining both ITS1 and ITS2 markers in our study increased species recovery by a substantial 10.5%, highlighting the advantage of using multiple primer pairs targeting the ITS region in fungal metabarcoding, but still did not reach complete species recovery.

Reference coverage for alternative markers was uneven and significantly lower than for ITS region. Moreover, a key limitation emerged with the rpb2b marker: while it is more suitable for Illumina sequencing due to its shorter length, it is located outside the portion of rpb2 gene commonly used in taxonomic studies. As a result, reference sequence availability for rpb2b is substantially lower compared to other markers. Despite these limitations, our results showed that combining ITS2 with one additional alternative marker can lead to a notable increase in species recovery, ranging from 6.7% to 11.8%, depending on the marker (Figure 3). To fully leverage the utility of alternative markers in fungal metabarcoding, development of more comprehensive and curated reference databases for targeted DNA regions are essential.

Inclusion of a third marker provided a further improvement of 0.8% to 5.3%, while adding a fourth or fifth marker resulted in only marginal gains of 0.1% to 1%. The highest additional value was achieved using the ITS2 marker with ef1‐α marker (87.1%) or with ef1‐α and ITS1 markers (92.2%). We suggest that combining two or three markers is an optimal number to maximise species detection in fungal metabarcoding.

This finding aligns with a recent study by Commichaux et al. (2024), which reached a similar conclusion for bacterial markers, pointing to the broader limitations of single‐marker approaches in metabarcoding across microbial domains. The benefits of synergistic multi‐marker approach were also shown in metabarcoding of other organisms (Arulandhu et al. 2017; Lekberg et al. 2018; Ogier et al. 2019; van der Loos and Nijland 2021; Espinosa Prieto et al. 2024). For instance, species recovery in metabarcoding of marine communities increased from 77% with a single marker to 89%–93% when combining markers (van der Loos and Nijland 2021). The incorporation of an additional marker in metabarcoding of plants led to a 7%–10% increase in species recovery, and using a combination of four markers and five primer pairs showed 97% recovery of species in a mock community (Espinosa Prieto et al. 2024).

Despite the benefits of a multi‐marker approach, significant technical limitations remain regarding its application to environmental DNA. Currently, different markers located at distinct genomic positions must be sequenced in independent runs from the same DNA template. This prevents the direct concatenation of sequences from a single individual or organism into a multi‐locus genotype. Identification across multiple ranks thus depends heavily on the existence of comprehensive reference sequences. However, even without single‐molecule linkage, parallel multi‐marker sequencing provides a critical corrective for community‐level biases, allowing researchers to trace and monitor target species with greater accuracy than a single‐marker approach would allow.

Reliability of RSA‐Based Estimation of Fungal Community Structure

4.3

A realistic representation of community structure in fungal metabarcoding is essential for the accurate interpretation of ecological patterns (Op De Beeck et al. 2014; Nguyen et al. 2015; Camacho‐Sanchez 2024). In our study, all five tested markers produced a broadly consistent and reliable quantitative profile of the mock community at the order level, with only a few major discrepancies observed among them. Primer pairs targeting ITS showed the most balanced amplification across all orders included in the mock community (Figure 4). This consistency suggests that higher taxonomic ranks are relatively robust to methodological biases. However, at the species level, read distributions were highly deviating across all markers, with remarkable marker‐specific variation in amplification efficiency and taxa representation. This is well documented by the list of the most dominant taxa, which showed little overlap between markers (Table 2).

This discrepancy between expected taxa abundances and their RSA estimates is a common theme in fungal metabarcoding. Bakker (2018) and McTaggart et al. (2019) similarly observed significant over‐ and under‐representation of taxa in their fungal mock communities mixed at equimolar DNA proportions, which they attributed to a combination of variable ITS region copy numbers and taxon‐specific PCR amplification biases. The challenge of quantification was further highlighted by Winand et al. (2025) and Shapkin et al. (2025), who also concluded that using read counts to estimate community structure performed poorly, with experimental RSA often deviating substantially from the expected. Although some studies have attempted to mitigate this problem by applying normalisation factors that account for species biomass, the correction of RSA estimates remains a difficult task with many hidden sources of bias (Luo et al. 2023).

In our study, the most striking was the high amplification of Entoloma species (> 90% of total sequences) by the ITS1 primers. This pattern was repeated in all three PCR runs, which strongly suggests the presence of preferential amplification bias. This is likely due to the design of the ITS1FI2 forward primer, as no similar issues were documented for the ITS2 reverse primer prior to this study (Bellemain et al. 2010; Toju et al. 2012; Monard et al. 2013; Op De Beeck et al. 2014). However, our in silico analysis of ITS1FI2 matching against sequenced fungal genomes did not confirm a preference for Entoloma over other Agaricomycetes fungi. The explanation for this phenomenon is therefore unclear. These results do not support the suggestion that alternative markers may provide a more accurate picture of fungal community structure (Větrovský et al. 2016). Instead, such heterogeneity of recovered species‐level profiles from the same mock community using different markers points to the dominant impact of accidental events, such as primer‐template compatibility, differences in amplification efficiency, and the stochastic nature of PCR amplification (Deagle et al. 2019; van der Loos and Nijland 2021; Mächler et al. 2021; Shapkin et al. 2025).

Considering these limitations, we tend to agree with the conclusion of Castaño et al. (2020), that metabarcoding should be regarded as a semi‐quantitative approach, particularly at the species level. Nevertheless, the application of multiple barcode markers can help to offset the individual biases of each marker and yield a more accurate and balanced picture of community structure.

Conclusions and Recommendations

4.4

While our findings confirm the dominance of the ITS region among fungal molecular markers, they also support the growing consensus that no single marker or primer pair can capture the full breadth of fungal diversity in metabarcoding. When using short‐read sequencing platforms where the forward and reverse reads do not overlap to form a single contig of the full ITS region, we suggest analysing the data from parallel sequencing of both the ITS1 and ITS2 subregions. Furthermore, we recommend incorporating at least one alternative fungal barcode marker to achieve higher species recovery and wider taxonomic resolution.

We demonstrated that the presumed lower specificity and degeneracy of the primers does not preclude their use in PCR on a multi‐template matrix of environmental DNA. Among the tested alternative markers, we recommend ef1‐α with the EF1‐1577F/EFgr primer pair for short‐read sequencing, and rpb2 with the bRPB2‐6F/bRPB2‐7R primer pair when longer‐read sequencing is available. Despite both primer sets demonstrating near‐universal PCR amplification and strong species recovery, their use in metabarcoding requires further methodological optimization to reduce sequence data losses due to non‐target amplification.

The broad application of ef1‐α and rpb2 genes in fungal metabarcoding remains limited, primarily due to their lower representation in public reference databases compared to the ITS region. Nevertheless, our results demonstrate that incorporating these genes alongside ITS can significantly improve species recovery, even in the context of uneven reference coverage. We recommend adopting a synergistic multi‐marker approach as a worthwhile and good practice in fungal metabarcoding, both to achieve more comprehensive community characterisation and to contribute to the long‐term enrichment of global reference datasets by sequences of less used markers.

The conclusions of this study are primarily applicable within a short‐read Illumina metabarcoding framework. It remains unclear how directly the reported marker performance can be extrapolated to metabarcoding studies that employ long‐read sequencing platforms, such as PacBio or Oxford Nanopore, relying exclusively on conventional Sanger‐derived reference databases. Future research should explore platform‐specific effects on reference accuracy and sequence resolution to further validate these multi‐marker synergies.

Author Contributions

Conceptualisation: Tomáš Zelenka, Petr Baldrian, Miroslav Kolařík; methodology: Vasilii Shapkin, Tomáš Zelenka, Tomáš Větrovský, Martin Kostovčík; formal analysis: Vasilii Shapkin, Tomáš Zelenka, Tomáš Větrovský, Martin Kostovčík; resources: Tomáš Zelenka, Ivana Eichlerová, Lucia Žifčáková, Jan Borovička, Michal Tomšovský, Miroslav Kolařík; writing – original draft: Vasilii Shapkin, Tomáš Zelenka, Miroslav Kolařík; writing – review and editing: Vasilii Shapkin, Tomáš Zelenka, Petr Kohout, Slavomír Adamčík, Petr Baldrian, Miroslav Kolařík; supervision: Miroslav Kolařík; funding acquisition: Slavomír Adamčík, Petr Baldrian, Miroslav Kolařík.

Funding

This work was supported by the Czech Science Foundation (grant no. 21‐17749S, 25‐17553S), Slovak Research and Development Agency (grant no. APVV‐24‐0473), Ministry of Education, Youth and Sports of the Czech Republic (project AdAgriF, no. CZ.02.01.01/00/22_08/0004635, and project Talking microbes, no. CZ.02.01.01/00/22_008/0004597) and Czech Academy of Sciences (Strategie AV21 project VP33 MycoLife—the world of fungi).

Disclosure

Benefit Sharing Statement: Fungal specimens used in this study were obtained from institutional collections within the Czech Republic (Table S1). Specimens originating from the Czech Republic were collected in accordance with national legislation. For specimens originating from other countries, we have exercised due diligence to confirm that these genetic resources were acquired in full compliance with the Nagoya Protocol and the national legislation of the provider countries. Benefits from this research are shared on a non‐monetary basis through the publication of these results, and the public deposition of all sequence data.

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Data S1: Supporting Information.