The field experiments were conducted in the Aoba-no-Mori Park and were carried out with the permission of the park administrator.

Stemflow was collected from a Japanese apricot (Prunus mume) in Aoba-no-Mori Park, Chiba City, Japan (35.5984N, 140.1395E), on a rainy day on July 1, 2023, from 08:30 to 17:00. Stemflow sampling was conducted in accordance with Sakata et al. (2023), with a modification to the gauze collection process. A round rubber rope (IT-9149, Itodai Seiko) and gauze (Insert Care Gauze No. 1, 63-1452-99, Hakujuji Co., Ltd.) were wrapped around the tree trunk. A silicone funnel (Silicone Funnel Icho, 07,438, Yamazaki) was attached to the rubber rope and gauze, with a hose connected to a backflow prevention unit. A 1-L plastic bag (DP16-TN1000, Cowpack LTD) was connected to this unit to collect stemflow from the tree trunk until the rain stopped.

In this modified method, the gauze was retrieved from the tree trunk and placed directly into the collected stemflow within a single-use plastic bag for filtration, rather than being soaked in Milli-Q water in a separate pre-sterilized 50 mL syringe as described in Sakata et al. (2023). The collected stemflow (approximately 1 liter in a plastic bag) underwent gravity filtration following Oka et al. (2022). The plastic bag’s cap was removed and replaced with a custom-made connecting attachment (Oka et al. 2022), which was linked to the inlet port of a Sterivex filter cartridge. A plastic tube was connected to the outlet port, and the plastic bag was suspended at approximately 2 meters to facilitate gravity filtration. Before filtration, the plastic bag was shaken approximately 10 times to ensure homogenization.

The stemflow was successively filtered into 15 Sterivex filter cartridges (pore size 0.45 μm; Merck Millipore, MA, USA), with each cartridge receiving approximately 50 mL of filtrate, as measured using a plastic measuring cup. For negative controls, purified water was filtered through three additional Sterivex filter cartridges, each receiving 50 mL. After filtration, 1–2 mL of RNAlater (Thermo Fisher Scientific, DE, USA) was added to each cartridge through the inlet using a disposable pipette (Nihon Medical Science, Osaka, Japan). The Sterivex filter cartridges were stored in a freezer at –20 °C until DNA extraction.

The 15 frozen filter cartridges along with 3 frozen cartridges designated as negative controls were left at room temperature for approximately 10 minutes. The remaining RNAlater solution in the cartridges was then removed following the method described by Minamoto et al. (2021), with slight modifications. To remove RNAlater, a 2 mL tube was placed inside a 50 mL conical tube, and the outlet port of each cartridge was connected to the 2 mL tube. The capped conical tube was centrifuged at 6,000× g for 1 minute. After centrifugation, an aspirator (QIAvac 24 Plus, Qiagen, Hilden, Germany) was used to ensure complete removal of any remaining liquid from the cartridge. The lichen fragments trapped on the filters were either processed without homogenization or homogenized before eDNA extraction using the DNeasy Plant Mini Kit (Qiagen, Hilden, Germany).

Without bead-beating homogenization (NB): This method was applied to five cartridges, which proceeded directly to DNA extraction without homogenization.

Bead-beating homogenization (BB): One gram of zirconia beads (φ 0.5 mm; AsOne, Tokyo, Japan) was added to 10 filter cartridges from the inlet using folded weighing paper. Then, 100 μL of AP1 buffer and 1 μL of RNase A stock solution were added to each filter cartridge containing zirconia beads. Five of these cartridges were attached to the Vortex Adapter 24 (Qiagen, Hilden, Germany) and vortexed for 3 minutes.

Frozen bead-beating homogenization (FB): The remaining five cartridges with zirconia beads were inserted into a customized tube block for an Automill (Tokken, Chiba, Japan) and soaked in liquid nitrogen for approximately 5 minutes. The frozen filter cartridges in the tube block were then homogenized at 2,000 rpm for 3 minutes using the Automill.

DNA extraction was performed from all 18 cartridges, including three designated as negative controls, using the DNeasy Plant Mini Kit (Qiagen, Hilden, Germany). After processing, either with or without homogenization, each cartridge was filled with 600 μL of AP1 buffer and 6 μL of RNase A stock solution. The cartridges were vortexed for 1 minute and incubated at 65 °C for 5 minutes while rotating on a rotator (Mini Rotator ACR-100, AS ONE, Tokyo, Japan); this process was repeated once. The cartridges were then brought to room temperature for approximately 10 minutes, after which 260 μL of P3 buffer was added, and the samples were placed on ice for 5 minutes. The lysate was extracted following a modified version of Miya et al. (2016). Specifically, each cartridge’s liquid content was transferred to a 2.0 mL tube using a centrifugation method. A 2 mL tube was placed inside a 50 mL conical tube. The inlet port of the cartridge was connected to the 2 mL tube, and the capped conical tube was centrifuged at 6,000× g for 1 minute to collect the lysate. From this point onward, DNA extraction followed the standard instructions provided in the DNeasy Plant Mini Kit.

The targeted lichen eDNA fragments were amplified using a two-step PCR process to prepare paired-end libraries for metabarcoding analysis on the MiSeq platform (Illumina, CA, USA) (Miya et al. 2015).

The first round of PCR (1^st PCR) was conducted with 35 cycles in a reaction volume of 12 μL, containing 6 μL of 2× PCR buffer for KOD FX Neo, 2.4 μL of 2 mM dNTPs, 1.4 μL of multiplexed primers (5 μM; ITS-PHLichenF, ITS-PHLichenR, ITS-PALichenF, and ITS-PALichenR; Sakata et al. 2023), 0.24 μL of KOD FX NEO DNA polymerase (Toyobo, Osaka, Japan), 1 μL of Milli-Q water, and 1 μL of template, with two technical replicates. The amplified DNA region targets a partial sequence of the nuclear ribosomal internal transcribed spacer 1 (ITS1) and a partial sequence of the 5.8S ribosomal RNA gene, with an expected fragment length of 145–174 bp. The thermal cycling profile, after an initial 2-minute denaturation at 94 °C, was as follows: denaturation at 98 °C for 10 seconds, annealing at 55 °C for 30 seconds, extension at 68 °C for 1 minute, and a final elongation at 72 °C for 7 minutes. Duplicate 1^st PCR products were pooled in a 1.5 mL tube and purified using SPRIselect (Beckman Coulter, Brea, CA). The concentration and product size distribution of purified PCR products were measured using a TapeStation 4200 (Agilent, Tokyo, Japan). The concentration of the pooled 1^st PCR products was adjusted to 0.1 ng/μL using Milli-Q water. Samples with concentrations less than 0.1 ng/μL were used as templates without dilution.

The second round of PCR (2^nd PCR) was conducted with 10 cycles in a 15 μL reaction volume containing 7.5 μL of 2× KAPA HiFi HotStart ReadyMix, 0.88 μL each of 2^nd-PCR forward and reverse primers (Miya and Sado 2019) with each different combination of indices (5 μM), 3.88 μL of Milli-Q water, and 1.88 μL of template. The thermal cycling profile, after an initial 1 minute denaturation at 95 °C, was as follows: denaturation at 98 °C for 20 seconds, and annealing and extension combined at 72 °C for 30 seconds, and a final elongation at 72 °C for 5 minutes. Indexed 2^nd PCR products were pooled, and the target-sized DNA was excised using E-Gel SizeSelect 2% (Thermo Fisher Scientific, MA, USA) with the E-Gel Pre-cast Agarose Electrophoresis System (Thermo Fisher Scientific). The concentration of the library was quantified using a Qubit dsDNA HS Assay Kit and a Qubit fluorometer (Thermo Fisher Scientific). The concentration of the library was adjusted to 4 nM using Milli-Q water. The library was sequenced using a MiSeq Reagent Kit v2 with 150 bp × 2 paired-end sequencing.

All raw DNA sequence data and associated information were deposited in the DDBJ/EMBL/GenBank database and are available under accession number DRR624633–624650.

Raw MiSeq reads were pre-processed and analyzed using PMiFish ver. 2.4 (https://github.com/rogotoh/PMiFish; Miya et al. 2020). Data pre-processing followed several steps implemented in USEARCH version 10.0.240 (Edgar 2010). Forward and reverse paired-end reads were merged using the “fastq_mergepairs” command. During this process, low-quality tail reads were discarded based on a cut-off threshold set at a quality (Phred) score of 2. Additionally, reads that were too short (<64 bp) after tail trimming and paired reads with too many differences (>5 positions) in the aligned region (approximately 70 bp) were also removed. Primer sequences were removed from the merged reads using the “fastx_truncate” command. Primer-trimmed sequences underwent further quality filtering using the “fastq_filter”. Low-quality reads with an expected error rate > 1% and shorter than 50 bp were removed using the “fastq_filter”. The filtered reads were dereplicated using the “fastx_uniques” command. During denoising, the “minsize_parameter” was set to 4, and singletons, doubletons, and tripletons were excluded from further analysis to avoid false positives, as suggested by Edgar (2010). The dereplicated reads were denoised using the “unoise3” command with a “minsize_parameter” of 4 and “unoise_alpha” parameter of 2.0, which removed putative chimeric and erroneous sequences prior to taxonomic assignment.

Reference sequences for taxonomic assignment (ITS) were collected to comprehensively represent the lichen species visually identified within the investigated plum grove. They were obtained primarily from lichen specimens that we ourselves identified and whose target nucleotide sequences we directly determined. These ITS sequences were obtained from materials collected in Chiba Prefecture, the same location as the study site, except for two Japanese local sites (Nagasaki and Akita). The latter specimens were not specifically derived from a prior survey; however the identification is indisputable, and its relevance to the study is guaranteed. Additionally, for each species, the top ~10 search results from NCBI (retrieved on March 10, 2023) were downloaded and saved in the FASTA format. To select appropriate sequences from the public database, the combined FASTA-formatted sequences were initially aligned using MAFFT version 7 (Katoh and Standley 2013) on the MAFFT web server (https://mafft.cbrc.jp/alignment/software/) with default parameters. A maximum likelihood (ML) tree was then constructed from the aligned sequences using the IQ-TREE web server (Trifinopoulos et al. 2016) with automatic selection of substitution model. Node support was estimated using bootstrapping (Felsentein 1985) with 1000 pseudoreplicates. The phylogenetic tree was visually inspected, and one or two sequences forming a group with many other sequences were selected for use as a reference database for taxonomic assignment (Table 1). Preference was given to nucleotide sequence data from Japanese and Korean species whenever possible.

Processed reads were assigned to each taxonomic assignments with > 98% sequence identity to the reference sequences (query coverage ≥ 90%, allowing two or three nucleotide mismatches) using the “usearch_global” command. Reads with sequence identities between 80–98% were designated as “U98 OTU” and clustered at 98% similarity using the “cluster_smallmem” command. After taxonomic assignment, all rare molecular taxonomic units (MOTUs) with read counts representing less than 0.01% of the total sequences per sample (<2 reads) were excluded from the taxonomic table to ensure conservative estimates of MOTU diversity (Miya et al. 2022).

Quantitative real-time PCR (qPCR) was performed using the ITS-PHLichenF and ITS-PHLichenR primer pair and eDNA samples with a StepOnePlus real-time PCR system (Thermo Fisher Scientific). The alternative primer pair (PALichenF and ITS-PALichenR) was excluded from the qPCR assay due to low amplification efficiency observed in preliminary experiments. Each 20 µl reaction contained 2 µl of eDNA template and 900 nM of each primer in PowerUp SYBR Green Master Mix (Thermo Fisher Scientific). The thermal cycling profile consisted of an initial 10-minute denaturation at 95 °C, followed by 40 cycles of 95 °C for 15 sec and 60 °C for 1 min. Three replicates were amplified for each eDNA sample, as well as for 100, 10, and 1 pg of standard DNA (a mixture of tissue-derived DNA from Kashiwadia orientalis, Hyperphyscia crocata, and Dirinaria applanata), and a negative PCR control. Among the species used in the standard DNA mixture, K. orientalis and D. applanata, which were primary target taxa, were also used for sequencing library preparation. The primers used in this study target a partial sequence of the internal transcribed spacer 1 (ITS1) and a partial sequence of the 5.8S ribosomal RNA gene.

To examine whether there were statistically significant differences in the number of detected species and the concentration of detected DNA among the three DNA extraction methods, Dunn’s multiple comparison test was performed. The test was conducted in R v4.2.2 (R CoreTeam 2022) using the dunn.test package (Dinno 2017). The regression analysis was performed using the mdl Rissanen method (Rissanen 1987) in R v4.2.2 (R Core Team 2022). To investigate the variation in species composition dissimilarity among the three eDNA extraction methods, the Jaccard dissimilarity index was calculated. All dissimilarity calculations were performed using the vegdist function in the vegan package (Oksanen et al. 2022) in R v4.2.2 (R Core Team 2022).

The pooled 102 samples (including three negative controls and 84 samples from other projects) were sequenced, and the MiSeq run yielded a total of 2,747,887 reads, with an average of 94.2% base calls having Phred quality scores of ≥ 30.0 (Q30; error rate = 0.1% or base call accuracy = 99.9%). This run was highly successful, considering that the manufacturer’s guidelines (Illumina Publication no. 7702011-001 as of May 27, 2014) recommend > 80% bases ≥ Q30 at 2 × 150 bp.

A total of 452,297 reads were assigned to the 15 samples, and the number of raw reads for each sample ranged from 16,486 to 45,390 with an average of 30,153 reads. Merging the two overlapping paired-end fastq files yielded 444,249 reads (98.2%). The sequences from which the primer sequences were removed were subjected to quality filtering to eliminate low-quality reads, resulting in 438,472 reads (a 96.9% retention). The remaining reads were dereplicated for subsequent analysis, and singletons to tripletons were removed from the unique sequences (Miya et al. 2022). Thereafter, the reads were denoised to remove putatively erroneous and chimeric sequences. The remaining 404,170 reads (a 92.2% retention) were subjected to taxonomic assignment, with 403,891 reads (99.9% of the denoised reads) putatively considered lichen sequences.

We analyzed these 403,891 reads (average 26,926 reads per sample) from the 15 samples, with five replicates for each of the three eDNA extraction methods. The average read counts for the three methods without bead-beating (NB), with bead-beating (BB), and with frozen bead-beating (FB) were similar, being 29,535, 28,608, and 22,635, respectively. Negative controls yielded no denoised reads across all methods.

Results from the automatic taxon assignments, with manual adjustments based on the genus-level ML trees, are summarized in Table 2. In total, 11 species from four families and eight genera were detected across the three methods, with a total detection frequency of 100.

The number of detected species per sample was the lowest in the eDNA extraction method without bead-beating (NB; Fig. 1), ranging from 2 to 6 species across the five samples, with a total of 7 species detected (Table 2). In contrast, the number was higher in the two methods with bead-beating (BB and FB), ranging from 7 to 9 species across the 10 samples, with a total of 9 species detected for each method (Fig. 1, Table 2). Statistically significant differences were found between the methods with and without bead-beating (NB vs. BB and NB vs. FB; p < 0.05, Dunn’s multiple comparison test), while no significant difference was found between the two bead-beating methods (BB vs. FB; p = 1.000, Dunn’s multiple comparison test) (Fig. 1).

Figure 1.

Box plots showing the variation in the number of detected species per sample among the three eDNA extraction methods. Numerals beside the boxes indicate the total number of detected species across the five samples for each method. Different letters denote significant differences. The ‘x’ inside the box represents the mean value (NB: 4.2, BB: 8.0, and FB: 7.8). The box edges represent the interquartile range (IQR), where the upper edge indicates the 75^th percentile and the lower edge indicates the 25^th percentile. The upper whisker represents the maximum value (NB: 6, BB: 9, and FB: 9), while the lower whisker represents the minimum value (NB: 2, BB: 7, and FB: 7). Different letters indicate significant differences between the DNA extraction methods p < 0.05).

As with the number of detected species (Fig. 1), the DNA concentration was the lowest in the eDNA extraction method without bead-beating (NB), with a mean of 7.2 ± 4.2 pg/2 µL (mean ± 1SD; Fig. 2). In the two methods with bead-beating (BB and FB), it was significantly higher, with means of 53.5 ± 15.3 and 53.1 ± 11.0 pg/2 µL, respectively (Fig. 2). Statistically significant differences were found between the methods with and without bead-beating (NB vs. BB and NB vs. FB; p < 0.01, Dunn’s multiple comparison test), while no significant difference was found between the two bead-beating methods (BB vs. FB; p = 1.000, Dunn’s multiple comparison test) (Fig. 2).

Figure 2.

Box plots showing the variation in DNA concentration, as determined by qPCR, among the three eDNA extraction methods. Different letters denote significant differences. The ‘x’ inside the box represents the mean value (NB: 7.2 pg/2 μL, BB: 53.5 pg/2 μL, and FB: 53.1 pg/2 μL). The box edges represent the interquartile range (IQR), where the upper edge indicates the 75^th percentile and the lower edge indicates the 25^th percentile. The upper whisker represents the maximum value (NB: 12.9 pg/2 μL, BB: 75.1 pg/2 μL, and FB: 64.7 pg/2 μL), while the lower whisker represents the minimum value (NB: 2.1 pg/2 μL, BB: 37.7 pg/2 μL, and FB: 39.4 pg/2 μL). Different letters indicate significant differences between the DNA extraction methods (p < 0.05).

As expected from the similar patterns of variation in the number of detected species (Fig. 1) and DNA concentration (Fig. 2), a significant positive correlation was found between the two variables, with a coefficient of determination (R²) of 0.707 (p < 0.01) (Fig. 3). These observations indicate that bead-beating residues from stemflow not only increase the amount of extracted DNA but also consequently lead to an increase in the number of detected species.

Figure 3.

Scatter plot showing the relationship between DNA concentration and the number of detected species.

Of the 11 detected species, 6 were consistently found in all bead-beating samples (BB and FB; detection frequency = 5; Table 2). In contrast, the detection frequencies of these 6 species in the 5 samples without bead-beating (NB) varied from 2 to 5, with a mean of 3.33 (Table 2). Additionally, two species, viz. Phaeophyscia rubropulchra and Physciella melanchra, which were undetected in only a few of the 10 samples in BB and FB, were not detected at all in NB (Table 2). Furthermore, the remaining three species were detected only once across the three methods. These patterns of detection frequencies were consistent with the overlap of detected species among the three extraction methods, as shown in Fig. 4.

Figure 4.

Venn diagram showing the overlap of detected species among the three eDNA extraction methods.

The average read number per sample for the 11 species was plotted against detection frequency (Fig. 5). The read number per sample for the two species with the highest detection frequency of 15 (Dirinaria applanata and Kashiwadia orientalis) was exceptionally high, both exceeding 10⁵. It ranged between 10³ and 10⁴ for the four species with intermediate detection frequencies (12 and 13), while it was below 10³ for the remaining five species with lower detection frequencies (1–9). These results indicate that detection frequency is reflected in the read number from MiSeq sequencing and suggest that sequencing depth should be increased for comprehensive detection of species composition in the extracted DNA (Gweon et al. 2019; Singer et al. 2019).

Figure 5.

Relationship between the detection frequencies of 11 species and the average read number per sample from MiSeq sequencing.

Fig. 6 shows the variation in pairwise dissimilarities (Jaccard distance) among the 15 samples (105 pairs), partitioned into three within-pre-treatment groups (NB, BB, and FB) and two between-treatment groups (BB vs. FB and NB vs. BB + FB). In the two pre-treatments involving bead-beating (BB and FB), which resulted in higher amounts of extracted DNA (Fig. 2), the dissimilarities were relatively low both within and between samples, with averages ranging from 0.1 to 0.2. In contrast, the average dissimilarity within samples that were not pre-treated with bead-beating (NB), which had significantly lower amounts of extracted DNA, was approximately 0.5, indicating a higher level of dissimilarity. Similar levels of dissimilarity were observed between samples with and without bead-beating pre-treatment (NB vs. BB and NB vs. FB).

Figure 6.

Box plots showing the variation in pairwise dissimilarities of species composition among the three eDNA extraction methods. The dissimilarities are partitioned into three within-pre-treatment groups (NB, BB, and FB) and two between-pre-treatment groups (BB vs. FB and NB vs. BB + FB). Different letters denote significant differences. The ‘x’ inside the box represents the mean value (NB: 0.493, BB: 0.161, FB: 0.117, BB vs FB: 0.142, and NB vs BB + FB: 0.500). The box edges represent the interquartile range (IQR), where the upper edge indicates the 75^th percentile and the lower edge indicates the 25^th percentile. The upper whisker represents the maximum value (NB: 0.667, BB: 0.250, FB: 0.222, BB vs FB: 0.250, and NB vs BB + FB: 0.778), while the lower whisker represents the minimum value (NB: 0.167, BB: 0, FB: 0, BB vs FB: 0, and NB vs BB + FB: 0.143). Different letters indicate significant differences between the DNA extraction methods (p < 0.05).

When environmental DNA is present at very low concentrations (e.g., due to dilution or degradation), PCR amplification can be less efficient and may fail to consistently amplify target sequences within the extracted DNA, leading to inconsistency in the detection of specific species (Deiner et al. 2017). This reduced detection reliability can result in variability in the number and types of species detected across samples, potentially increasing the perceived dissimilarity between them (Kelly et al. 2014). Consequently, the detected species may fluctuate randomly due to low template concentration, reducing the similarity between samples (Ficetola et al. 2016).

Of the 6 foliose lichen species visually observed on the tree, 4 species were detected in the BB and NB samples, while 5 species were detected in the FB samples via eDNA metabarcoding.

Based on their frequency and coverage as determined by visual observation, we classified them into three groups: “dominant” (Dirinaria applanata and Kashiwadia orientalis), “common” (Parmotrema austrosinense, P. tinctorum, and Flavoparmelia caperata), and “rare” (P. clavuliferum).

The two “dominant” species were consistently detected across all samples from the three methods. Among common species, Parmotrema austrosinense and P. tinctorum were detected more frequently in BB and FB samples than in NB samples, suggesting that the choice of sampling method affects detection efficiency. The absence of Flavoparmelia caperata across all methods may be attributed to DNA degradation, primer mismatches, or differences in cell wall composition affecting DNA extraction efficiency. Our results demonstrate that eDNA metabarcoding successfully identified most of the visually observed foliose lichen species, although detection rates varied across different sample types.

Additionally, eDNA metabarcoding detected by 4 additional foliose lichens (Hyperphyscia cf. adglutinata, H. crocata, Phaeophyscia rubropulchra, and Physciella melanchra) and two crustose lichens (Amandinea cf. punctata and Lecanora cf. fulvastra).

Notably, excluding species detected only once or twice, non-target lichen species in NB samples were either undetected or detected at lower frequencies compared to BB and FB samples.

Our reference database may appear small compared to standard databases often used in eDNA metabarcoding studies (Table 1). We acknowledge the limitations of using a small reference database for taxonomic identification. However, in the case of lichens, many publicly available sequences have uncertain taxonomic assignments. Expanding the database indiscriminately would increase the risk of erroneous taxonomic classifications. As a matter of fact, a previous study supports the relevance of our curated reference set (Sakata et al. 2023). Since our study focused on a tree where only 6 foliose lichens target species that were already known and present in a region with relatively low species diversity were visually observed, the use of a large database was not deemed necessary.

Furthermore, the number of unmatched sequences was minor in terms of both read count and species diversity, suggesting that this does not significantly impact the validity of our findings. These unmatched sequences were further examined using BLAST searches. The results showed that sequences with “no significant similarity found” accounted for a total of 228 reads (0.06% of the total reads) across 13 sequences from 6 samples. Additionally, one algal sequence was detected in one sample [12 reads (0.003% of the total reads)], while two non-lichenized fungal sequences were identified in 2 samples, totaling 23 reads (0.006% of the total reads).

This study highlights several methodological considerations for lichen eDNA metabarcoding. (1) Environmental variability: In rivers, it is suggested that the detectability of eDNA across varying densities is influenced by factors such as species, stream size, discharge rate, and season (Goldberg et al. 2011). Similarly, in stemflow, factors such as rainfall intensity may influence eDNA collection, potentially impacting comparability across different sampling periods and locations. Standardizing sampling conditions remains a key challenge for future studies; (2) Sampling scope: Our study focused on a single tree, which limits the applicability of the results. However, since the primary aim was to compare extraction methods rather than characterize lichen diversity, this limitation does not necessarily compromise our methodological conclusions. Future studies incorporating multiple trees will help assess potential site-specific effects; (3) Potential methodological biases: Differences in bead-beating intensity are known to have a significant impact microbial community composition (Albertsen et al. 2015). However, no method-specific biases have been identified so far in the analysis of lichen species composition. Additionally, the application of bead-beating for detecting metazoans without cell walls may negatively affect DNA recovery and detection efficiency. Further investigation is required to evaluate how these methodological choices influence ecological interpretations; (4) Broader applicability: Although bead-beating is widely used for eDNA extraction (Kuske et al. 1998; Guo and Zhang 2013; Albertsen et al. 2015; Lever et al. 2015; Ushio 2019), differences in homogenization protocols across studies make direct comparisons challenging. Future research should explore the consistency of these methods across diverse environmental contexts; and (5) Low template DNA concentration: Our template DNA concentrations were relatively low (<100 pg/µL), which makes PCR amplification highly sensitive to random fluctuations, primer binding efficiency, and potential contamination. However, due to constraints in our experimental system and sample availability, working with higher DNA concentrations was not feasible.

This study demonstrated that bead-beating as a pre-treatment for eDNA extraction efficiently increases both the yield of lichen DNA and the number of detected species from environmental samples. In contrast, the use of liquid nitrogen freezing before eDNA extraction, which is effective for DNA extraction from lichen tissue fragments (Cubero et al. 1999), did not show any significant effect on the results of eDNA metabarcoding. Considering the time, effort, and cost associated with freezing, bead-beating alone appears to be sufficient as a pre-treatment for extracting lichen eDNA. However, further research is needed to determine whether these findings are consistent across different environmental samples and other arboreal organisms. The application of this optimized method could enhance the accuracy and efficiency of biodiversity monitoring and conservation efforts involving lichens.

The authors have declared that no competing interests exist.

No ethical statement was reported.

This study was supported by JSPS KAKENHI grant 23K11516.

AS and MM conceived and designed the study. AS selected and provided appropriate materials and performed the field survey. AS, TM, TS, and MM conducted the laboratory experiments and data analysis. ROG reanalyzed datasets. AS, TM, and MM wrote and edited the first draft of the manuscript. All authors discussed the results and contributed to the development of the manuscript.

Ayumi Sakata https://orcid.org/0009-0005-2153-7119

Toshifumi Minamoto https://orcid.org/0000-0002-5379-1622

Tetsuya Sado https://orcid.org/0000-0003-4149-9824

Ryo O. Gotoh https://orcid.org/0009-0006-3816-4774

Masaki Miya https://orcid.org/0000-0002-9791-9886

All of the data that support the findings of this study are available in the main text or Supplementary Information.