The Long chain Diol Index: A marine palaeotemperature proxy based on eustigmatophyte lipids that records the warmest seasons

In order to gain insight into the seasonal occurrence of LCDs and their source organisms, Mediterranean Sea surface water samples were taken from mid-spring to early autumn in 2019 (Table 1). The coldest in situ temperatures were measured on 30 April, when sea surface temperatures (SSTs) ranged between 14.3 and 16.0 °C, and the warmest SSTs of 26 °C were measured on 27 August. On the other sampling days, SSTs ranged between 20 and 22 °C. Measured particulate organic carbon (POC) values gradually decreased from 248 µg L⁻¹ in April to 167 µg L⁻¹ in August, with a minor increase to 198 µg L⁻¹ in September (SI Appendix, Table S1).

LCDs were detected as free lipids in all samples from the Mediterranean Sea, and concentrations were not affected by acid or base hydrolysis (see SI Appendix, Supplementary Note 1). LDI-LCD concentrations increased substantially between 2 and 27 August, when a strong increase in SSTs was observed (Fig. 1 A and G, Table 1, and SI Appendix, Table S1). In accordance with the SSTs, LDI-LCD concentrations on 17 September were decreased again, suggesting that the sources of LDI-LCDs thrive under warm water conditions. In the Mediterranean Sea water, LDI-LCD concentrations, as well as C₃₀ 1,15-keto-ol concentrations, followed a distinct pattern over time, reaching maximum values for 27 August (Fig. 1 A and B and SI Appendix, Supplementary Note 2). These concentrations were in the same range as sterols (Fig. 1 D–F), suggesting a relatively high abundance of LDI-LCD sources at that time. Comparison of LDI-LCD concentration patterns with p23S rRNA read numbers for individual amplicon sequence variants (ASVs) was performed to ascertain the LDI-LCD sources. The ASV numbering is based on the order of total read counts. Read numbers of ASV_003 (i.e., the third most abundant ASV) showed the strongest resemblance to the LDI-LCD concentrations (Fig. 1 and SI Appendix, Tables S1 and S2) and was assigned to the cyanobacterium Synechococcus sp. ASV_001 and ASV_002 were also assigned to Synechococcus cyanobacteria, but their reads showed high numbers on 25 June and 2 August, when LCD concentrations were low. An explanation is that only a select group of Synechococcus species produces LCDs, something similar as observed for the diatoms, where only Proboscia spp. produce 1,14-LCDs (7).

ASV_005, which was not directly assigned to a specific organism by the basic local alignment search tool (BLAST) procedure, also showed a pattern similar to LDI-LCD concentrations (Figs. 1 and 2 and SI Appendix, Tables S1 and S2). These genes were absent in samples from 30 April and 25 September, whereas 64 reads were detected on 2 August, and 868 reads were recorded on 27 August. The highest number of ASV_005 reads was identified on 17 September, both for absolute numbers (1,117 reads) and relative to the total reads per sample (13%). Comparing this gene sequence with data from GenBank (https://www.ncbi.nlm.nih.gov/genbank/), the closest available organism, at 91.6% identity was the eustigmatophyte Nannochloropsis oceanica strain LAMB2011 (accession number CP038136). To get a more accurate identification of ASV_005, we performed phylogenetic analyses on our p23S rRNA dataset combined with relevant corresponding sequences obtained from GenBank. In the obtained phylogenetic trees (Fig. 2 and SI Appendix, Figs. S2–S5), ASV_005 formed a cluster with ASV_026 and ASV_031 which branched off from the known eustigmatophytes before their species diversified. Separating eukaryotes from the prokaryote reads, ASV_005 contributed 76% of all eukaryote reads for 27 August and 59% of the reads for 17 September, a pattern similar to the LDI-LCD concentrations (Fig. 2). ASV_005 was the only ASV tentatively linked to eustigmatophytes for 27 August; ASV_026 contributed an additional 27 and 3% of all eukaryote reads on 2 August and 17 September, respectively; and ASV_031 contributed 1, 9, and 5% of all eukaryote reads for 30 April, 25 June, and 2 August, respectively. These results suggest a shift in LDI-LCD producing species over time, which could be related to SST, nutrient availability, or other environmental conditions. Given the similarity in sequence reads and LDI-LCD concentrations over time and the occurrence of LDI-LCDs in eustigmatophytes, these algae are a more likely source for the LDI-LCDs than cyanobacteria.

In the calculated phylogenetic trees based on p23S sequences, the distinctive ASVs clustered together with six gene fragments extracted from Holocene and Late Pleistocene sediments from the salt lake Qinghai (39) and with one operational taxonomic unit (OTU) from water column samples from the Mariana Trench (40) (Fig. 2 and SI Appendix, Figs. S2–S5, cluster b). Indeed, LDI-LCDs have previously been reported in sediments from both environments (41, 42). In the study of the Mariana Trench, eustigmatophyte genes were reported from the bathypelagic and hadal zones (40). We followed the procedures described in that study to cluster their 23S rRNA data into OTUs (40) but used the improved USEARCH version 11.0 instead of version 7.0 (43). With that approach, we obtained only one eustigmatophyte OTU (OTU_18) clustering with Nannochloropsis species instead of two OTUs that were reported to have high read numbers at 8,320- and 4,000-m depth. In addition to this Nannochloropsis-related sequence, we obtained OTU_56, clustering within the eustigmatophyte sister group identified in our Mediterranean Sea study. More than 85% of this OTU was recorded in the upper epipelagic zone at 4-m depth, and the remaining reads were distributed over all other depths.

Whereas a link between eustigmatophytes and LDI-LCDs seems highly plausible in the current work, the study by Balzano and coworkers (27) aimed at identifying the sources of LCDs in the western tropical North Atlantic Ocean and detected eustigmatophyte reads in only 2 out of 22 surface water samples, and the read numbers were in no proportion to the measured LDI-concentrations. Acknowledging the limitations of the automated BLAST results, we performed phylogenetic analyses on their DNA data for improved identification (Fig. 3 and SI Appendix, Figs. S6–S9). The dataset from Balzano et al. was combined with 18S DNA data from 19 sampling sites across the coastal areas of the Mediterranean Sea (44) and with 18S V4 region sequences from eukaryote species obtained from GenBank. Our results confirm the identification of denovo2075 as a species closely related to Nannochloropsis. This species was exclusively found in the surface waters of station 12, one of the shallowest stations located on the continental shelf near the coast (27). Our results do not support the identification of denovo229 as an unspecified eustigmatophyte. On the other hand, our phylogenetic analyses suggest that denovo162, denovo1571, and denovo2431 belonged to an early-branching clade of eustigmatophytes. Denovo162 was identified in almost half of the surface and mixed-layer (15-m depth) water samples from locations farther off the coast (27), but nevertheless, the number of reads is far too low to link the detected LCD concentrations directly to these species. Based on the concentrations of LCDs detected in N. oceanica (45), Balzano et al. indicated that even if all of the nanoplankton organisms detected in their study produced LCDs, the total cell number would still have been too low to explain the measured LCD concentrations (27). It was therefore hypothesized that the LCDs originated predominantly from degraded material whose DNA was not preserved.

In the dataset from Penna et al., eustigmatophyte sequences were detected in 18 out of 20 Mediterranean Sea samples (44). An additional BLAST search in GenBank targeting the early-branching eustigmatophyte sequences revealed several other related sequences from various marine settings. These results support the existence of a globally present, yet still unidentified, group of marine eustigmatophytes that originated before the known clade of predominantly freshwater eustigmatophytes diversified. Phylogenetic analyses on an 18S rRNA dataset consisting of sequences with >1,000 base pairs further substantiate the existence of this early-branching group of marine eustigmatophytes (SI Appendix, Figs. S10 and S11). In all 18S rRNA phylogenetic trees, the potential LDI-LCD producers cluster as a sister group to the currently known eustigmatophytes (Fig. 3 and SI Appendix, Figs. S6–S11), and the same is observed for the Maximum Likelihood phylogenetic tree based on a subset of p23S rRNA sequences (SI Appendix, Fig. S5). Bootstrap or posterior probability (PP) values for the p23S rRNA phylogenetic trees in Fig. 2 and SI Appendix, Figs. S2–S4 do not reject this possibility. However, there is still the possibility that the clade of potential LDI-LCD producers represents another lineage of stramenopile algae independent of Eustigmatophyceae. Additional sequence, morphological, and physiological data are still required to test if these LDI-LCD producers meet the criteria needed for inclusion in the Eustigmatophyceae.

At our sampling locations in the Mediterranean Sea, LDI-LCDs were predominant at the time with the highest SSTs (Fig. 1 A and G), which is consistent with earlier observations that the LDI temperature proxy shows the strongest correlation with summer SSTs (34, 35, 46). To explore this in more depth, we combined the surface sediment LCD datasets from a number of studies (34, 47–49) to compare LCD distributions with mean annual, seasonal, and monthly values of temperature, nutrients, salinity, and mixed-layer depths (Fig. 4) at a 1° and 0.25° grid resolution, obtained from the World Ocean Atlas 2018 (50). We observed the strongest correlations between the LDI and surface water temperatures when using 1° grid datasets with data from all available years combined; if no data directly corresponding to the sample coordinates were available, data from the nearest adjacent grid cell were selected. Neither the use of data from 0.25° grid datasets nor the use of datasets from specific time spans resulted in stronger correlations.

We observed the strongest correlations for the LDI with SSTs from the warmest month at each of the sampling locations (Figs. 4A and 5 A and B). The water conditions during these months were characterized by, on average, low surface water salinities, low nutrient concentrations, and shallow mixing-layer depths (Fig. 4). These findings further corroborate that LDI-LCD producers mainly occur in the upper water layers during periods of enhanced stratification. The results also support the rationale of the 1,14-LCD seasonal upwelling indices; these are based on the assumption that 1,14-LCDs are mainly produced during periods of seasonal upwelling and high nutrient concentrations, whereas the LDI-LCDs are linked to nonupwelling periods (15, 16). The predominance of LDI-LCD producers in surface waters, as proposed in previous LDI calibration studies (34, 35), is also reflected by the predominance of OTU 56 from the Mariana Trench in samples at 4-m depth (40) and the highest LCD concentrations and numbers of eustigmatophyte reads in the upper water layers of the western tropical North Atlantic Ocean (27). Furthermore, the predominance of LDI-LCD producers at times with the shallowest mixed-layer depths matches the observation that in Lake Geneva, LCDs were most abundant during the period of thermal stratification of the water column (51). The low nutrient conditions prevailing during these months are consistent with the finding that eustigmatophytes were predominant in Qinghai Lake during periods with oligotrophic conditions (39).

The LDI–temperature correlation greatly improved when selecting SSTs from the warmest months instead of annual mean SSTs (Figs. 4A and 5 A and B), but the significant scatter that remained suggests there are additional factors besides temperature affecting the LDI. It has been suggested that contributions of 1,13-LCDs from specific Proboscia diatom species may affect the LDI (52), and the contribution of these species can be quantified by the fractional abundance of C₂₈ 1,12-diol (FC₂₈ 1,12-diol) (34):

So far, few studies have reported FC₂₈ 1,12-diol values. For the 525 samples reported in the study by de Bar et al. (34), high FC₂₈ 1,12-diol values were predominantly identified in areas with mean SSTs from the warmest months ranging between 9 and 14 °C (e.g., Fig. 5 A and B), and their abundances showed no clear trend with temperature, salinity, or nutrients. The correlation of the LDI with mean SSTs from the warmest months increased from R² 0.88 to 0.90 by excluding 18 samples with FC₂₈ 1,12-diol values >0.15, while excluding additional samples with lower FC₂₈ 1,12-diol values had little effect (SI Appendix, Fig. S12A).

It has also been indicated that freshwater input or low salinity may affect the LDI (3, 32, 33, 35), and it was shown that the LDI–temperature correlation can be improved by excluding samples from areas with annual mean salinities <32 ppt (34). For the combined LDI dataset, excluding samples with annual mean or warmest month salinities <32 ppt resulted in an R² increase from 0.91 to 0.92 (SI Appendix, Fig. S13). This minor improvement of the correlation may not justify the considerable number of excluded samples (79 and 119, respectively, from a dataset of 688 samples in total). Moreover, exclusion of samples deposited under low salinity conditions is often not feasible for paleoclimate studies where salinity data are typically unavailable. When dividing the dataset into sample sets with annual mean salinities above and below 32 and correlating LDI values of each group individually to annual mean SSTs, it is remarkable that the slopes of the two correlations are almost identical, whereas the intercept for the lower salinity samples is 0.17 higher (SI Appendix, Fig. S13). Inferring LDI-derived SSTs using the >32 salinity correlation, thus, results in temperatures that are more than 5 °C higher than when using the <32 salinity correlation. This can be interpreted as an indication that also in low-salinity areas, the LDI is affected by temperature but reflecting colder-than-average periods of the year. When using the SSTs from the warmest months instead of annual mean SSTs, the offset between the two LDI calibrations is no longer constant, and the sample sets are more comparable to each other, which results in the smaller effect on the correlation when excluding low-salinity samples (SI Appendix, Fig. S13).

It has been noted that LCD distributions from lacustrine environments and near rivers are often characterized by relatively high C₃₂ 1,15-diol abundances (2, 3, 30, 32, 33, 53, 54). Whereas FC₃₂ 1,15-diol values were <0.03 in our Mediterranean Sea samples, values >0.1 have been reported in central North and South Atlantic Ocean sediments from areas that are unlikely to be affected by freshwater input (34). This suggests that at least some marine LDI-LCD producers also biosynthesize C₃₂ 1,15-diol. In a previous study, it was reported that the removal of samples with high C₃₂ 1,15-diol abundances had little effect on the LDI–temperature calibration (34). However, samples with high FC₃₂ 1,15-diol values seem to have a different correlation with SST compared to samples with lower FC₃₂ 1,15-diol values (Fig. 5 A and B). We have arranged the data of the combined LCD dataset in the order of their FC₃₂ 1,15-diol values and calculated moving slopes and intercepts for sets of 21 samples along the list. Data from sample sets with temperature ranges <15 °C were excluded, and the remaining values were plotted against their average FC₃₂ 1,15-diol values (Fig. 5 C and D). For samples with FC₃₂ 1,15-diol values <0.3, the calculated slope values showed a weak negative trend and the corresponding intercept values a weak positive trend with increasing FC₃₂ 1,15-diol, but those trends were of no proportion to the nonlinear variation observed for those samples. For samples with FC₃₂ 1,15-diol values >0.3, slope values showed a distinct decrease with higher FC₃₂ 1,15-diol, whereas the intercept showed a distinct increase.

Based on these observations, we propose to use the average slope of 0.033 and the average intercept of −0.008 when reconstructing SSTs for samples with FC₃₂ 1,15-diol values <0.3. For samples with FC₃₂ 1,15-diol values >0.3, the slope and intercept for each sample can be calculated individually, using the following equations:

Using these values, SSTs from the warmest months can then be estimated by

For the resulting dataset, estimated versus measured SSTs for the warmest months have an R² of 0.92 (Fig. 5E). For samples with FC₃₂ 1,15 diol <0.3, the SD for estimated − measured temperatures is 2.4 °C, and for samples with FC₃₂ 1,15-diol >0.3, this SD is 4.6 °C. Sample MT20 from the Gulf of Lion, showing an anomalously high LDI-inferred temperature of 34.3 °C, stood out from the rest of the dataset (estimated − measured temperature >3 × SD). In the original paper, this sample was left out when displaying the spatial distribution of C₃₂ 1,15-diol (33), possibly for showing an anomalous FC₃₂ 1,15-diol value. By excluding this sample from our dataset, the SD for estimated − measured temperatures is reduced to 3.3 °C. The LDI temperature calibration for samples with FC₃₂ 1,15 diol values >0.3 is currently based on only 33 samples, with a relatively large contribution of Baltic Sea samples (8) and Kara Sea samples (7); therefore, one should remain critical when including the results from samples with high FC₃₂ 1,15-diol values.

We also used Bayesian statistics following MATLAB code that has been used for previous proxy-calibration studies (55) to calibrate the relationship between the LDI and annual mean SSTs, mean summer SSTs, or mean SSTs from the warmest months for samples with FC₃₂ 1,15-diol values <0.3. The obtained slopes, intercepts, and root-mean-squared errors for the obtained Bayesian calibration models were similar to those obtained with the linear regression models. For samples with FC₃₂ 1,15-diol values <0.3, LDI temperatures that had previously been inferred with the original calibration (35) can simply be adjusted to warmest month temperatures by adding 2.64 °C to the initial value.

As evident from their identical slope values, LDI calibrations to annual mean SST and mean SST from the warmest month both have a 30.3 °C temperature range, but the warmest month temperatures are more evenly distributed along the calibration line (Fig. 5 A and B). Perhaps the most obvious improvement is the plateau with LDI values 0.95 and corresponding annual mean SST values ranging from 15.5 to 30.3 °C, whereas for the new calibration, >0.95 values were observed in areas where the maximum month temperatures ranged between 23.3 and 30.3 °C. An alternative calibration has been proposed for South China Sea samples with annual mean SST values >27 °C, where a negative correlation was observed between the LDI and annual mean SST (49). When correlating the LDI data to summer temperatures, a weak positive correlation over a very narrow range of measured temperatures had been observed. For our new calibration, the residual errors of those samples range from −2.5 to 1.8 °C with no obvious trend, thereby disputing the added value of the alternative calibration.

Temperatures calculated from the LDI values of our Mediterranean Sea water samples exceeded the measured SST (SI Appendix, Fig. S14A). Although the highest LDI values were indeed obtained at the warmest sampling day, the corresponding temperature variations were in no proportion to the measured data. Possibly, organisms living deeper in the water column have different LDI-LCD distributions, and the LDI preserved in sediments may be the result of mixed signals, but LDI-derived temperatures from surface sediment samples from the Gulf of Lion—near our sampling locations—also generally overestimated measured temperatures (SI Appendix, Fig. S14B). Some sediment trap studies also reported minimal seasonal LDI variation and LDI-derived temperatures that did not reflect the temperatures of the overlying waters, and similar observations are known for other organic temperature proxies (56). For LDI values obtained in Mozambique Channel sediment trap samples, it was hypothesized that resuspended material affected the sediment trap signal (56), but for our current study, the similar patterns observed for LDI-LCDs and eustigmatophyte sequence reads suggest that LDI-LCDs in our Mediterranean Sea samples were predominantly derived from fresh biomass.

For other well-established paleotemperature proxies, it often remains uncertain how exactly the sediment distributions originate. Sediment trap studies on the ${\text{U}}_{37}^{\text{K'}}$, an index based on the proportion of di-unsaturated to triunsaturated C₃₇ alkenones (57, 58), show that there is no consistent seasonal pattern in alkenone fluxes (59), and yet, the global ${\text{U}}_{37}^{\text{K'}}$ values strongly correlate with mean annual SSTs (60). Moreover, the global ${\text{U}}_{37}^{\text{K'}}$ calibration is statistically similar to the first published calibration derived from culturing experiments with the haptophyte alga Emiliania huxleyi (61). It has been observed that lipid distributions in different alkenone-producing species are differently affected by temperature, and correlations of the ${\text{U}}_{37}^{\text{K'}}$ for freshwater species do not correspond to the global marine ${\text{U}}_{37}^{\text{K'}}$ calibration (62–64). Eustigmatophyte culturing experiments have also demonstrated that these algae do adjust their LCD composition to the growth temperature (3). However, different eustigmatophyte families show different responses to temperature, and so far, none of the obtained LDI–temperature correlations corresponded to the marine LDI calibration (3). This demonstrates the relevance of isolating species from the eustigmatophyte clade for culturing experiments to further test the LDI calibration.

Surface water samples (0-m depth) were taken at four locations near Toulon, at a maximum distance of 5.6 km from the coast (Table 1 and SI Appendix, Fig. S1). Med. 2 and 3 were located farthest from the coast into the Mediterranean Sea, and Med. 4 was located at a shallower location at the edge of the Gulf of Lion. Surface water was collected using a bucket and, after measuring the temperature, stored for transport in 20-L Nalgene polycarbonate carboys. SPM matter was filtered within 2 d using precombusted 0.7 µm GF/F filters (Whatman). Filters used for lipid analyses were stored at −20 °C, and samples for DNA analyses were stored at −80 °C.

For the determination of POC, aliquots of freeze-dried filters were analyzed with a LECO RC 612 multiphase carbon analyzer. Filters were placed in silver cups and heated from 120 to 1,000 °C in a stream of O₂ (750 mL/min). The time- and temperature-resolved CO₂ release was detected with an infrared absorption cell.

Filters containing SPM from 8.7 to 29.8 L water were freeze-dried and ultrasonically extracted with 3× methanol (MeOH), 3 mL, 3× MeOH: dichloromethane (DCM) 1:1, 3 mL, and 3× DCM, 3 mL. The extracts were combined and dried under N₂. Milli-Q water was added to remove the salts, lipids were extracted with DCM, and samples were split into two equal aliquots and dried under N₂. Milli-Q water and one drop of 6% KOH in MeOH were added to one aliquot of each sample, and neutral lipids were extracted with n-hexane (3×). Thereafter, the water was acidified to pH 1 with 10% HCl, and polar lipids were extracted with n-hexane (5×). We added 3 mL 6% KOH in MeOH in a 7-mL vial to the remaining half of the sample extracts and to the extracted filters for base hydrolysis. The vials were closed and heated for 3 h at 80 °C in a heating block. After cooling down, neutral lipid fractions were extracted with n-hexane (3×), the remaining samples were acidified to pH 1 with 10% HCl, and polar lipids were extracted with n-hexane (5×). A second set of filters containing SPM collected 27 August 2019 and 17 September 2019 were extracted as described above, after which 3 mL 2n HCl in MeOH were added to half of the dried extracts and to the extracted filters for acid hydrolysis. Sample vials were closed and heated for 3 h at 80 °C in a heating block, and after cooling down, lipids were extracted with n-hexane (5×). C₂₂ 5,17-diol was added as an internal standard (52) to each of the obtained lipid fractions, and the solvent was evaporated under a stream of N₂. Fatty acid groups in the polar fractions were methylated according to Poerschmann et al. (65); 1 mL trimethylchlorosilane (TMCS) MeOH (1:9) was added to the dried extracts, and the sample vials were closed and heated for 90 min at 80 °C in a heating block. After cooling down, lipids were extracted with n-hexane (3×) and dried under N₂. Prior to gas chromatography–mass spectrometry (GC-MS) analyses, alcohol groups of all samples were silylated by adding 15 µl BSTFA [N,O-bis(trimethylsilyl)trifluoroacetamide] containing 5% TMCS and 15 µl pyridine and heating for 30 min at 60 °C. After cooling down, 100 µl n-hexane were added to all samples.

GC-MS analyses were performed using a Thermo Fisher Trace 1310 GC coupled to a Thermo Fisher Quantum XLS Ultra MS. The GC was equipped with a capillary column (Phenomenex Zebron ZB-5, 30-m length, 0.25-μm film thickness, inner diameter 0.25 mm), using helium as a carrier gas (flow rate of 1.5 mL/min). Derivatized extracts were injected into a splitless injector and transferred to the GC column at 270 °C. The GC oven temperature was kept at 80 °C for 1 min and then increased to 310 °C at 5 °C/min, at which it was held for 20 min. Electron ionization MS were recorded in full scan mode at 70 eV, a mass range of m/z 50 to 600 and a scan time of 0.42 s, or in selected ion monitoring mode of the characteristic fragmentation ions of the environmental LCDs and the LCD standard (i.e., m/z 187, 299, 313, 327, and 341; refs. 35 and 52) for the quantification of the LCDs. Correction factors of 0.124, 0.068, and 0.230, obtained from full scan MS, were applied for calculating the concentrations of the saturated and unsaturated environmental LCDs and the LCD standard, respectively. As small amounts of the neutral lipids were identified in the polar fractions—in particular, after base hydrolysis—lipids were quantified in both fractions, and the results of this study consist of the combined data.

With the methods used for the current study, C₂₉ 12-hydroxy methyl alkanoate coeluted with C₂₈ diols, complicating the identification of C₂₈ 1,12-diol (see SI Appendix, Supplementary Text).

On ice, 1/4 sections of the GF/F filters from location Med. 3, containing SPM from ca. 10 L seawater, were cut into small pieces using sterile scalpels and tweezers. Filter pieces were transferred into 2 mL microtubes, and DNA was extracted using a DNEasy Powerlyzer PowerSoil Kit (Qiagen) following the manufacturer’s instructions. The universal primers NGS-p23SrV_f1 (5′-TCG TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG GGA CAG AAA GAC CCT ATG AA-3′) and NGS-p23SrV_r1 (5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GTC AGC CTG TTA TCC CTA GAG-3′) (66) were used to amplify the plastid 23S rRNA as described by Sherwood and Presting (67). Sequencing was conducted at the G2L: Next-Generation Sequencing at the Göttingen Genomics Laboratory, Department of Microbiology, University of Göttingen (Germany). All the sequence data generated in this study have been deposited in NCBI BioSample database under the accession numbers SAMN24256204, SAMN24256205, SAMN24256206, SAMN24256207, and SAMN24256208. Eustigmatophyte sequences ASV_005, ASV_026, ASV_031, and ASV_157 have also been deposited individually under the accession numbers OL958649, OL958650, OL958651, and OL958652, respectively.

We obtained 1,531,336 raw sequences, which were further processed using appropriate packages in the R software environment (68). Sequences containing an “N” and sequences without primer detection were discarded, using the DADA2 (69), Biostrings (70), and ShortReads (71) packages. Thereafter, primers were removed from the reads using the Cutadapt software (72). DADA2 quality filtering was performed using a maximum expected error threshold of 2 for the forward reads and 4 for the reverse reads, and the dereplication of identical reads, sample inference, merging of paired-end reads, and chimera removal were executed using default values (69). The taxonomic assignment of the ASV sequences was performed using the nucleotide–nucleotide basic local alignment search tool (BLASTN) software version 2.10.1+ (73) with the BLASTN algorithm and the June 2021 NCBI GenBank nucleotide database version (74), using the Tidyverse package (75).

In order to verify sequence identifications, phylogenetic analyses were performed by Bayesian Inference algorithms, using the program MrBayes 3.2.7-WIN (76). The generalized time-reversible (GTR) model was used, with gamma-distribution rate variation across sites and a proportion of invariable sites, whereby five hot chains were run in addition to the standard (cold) chain. The Bayesian search was run for 2,500,000 generations, saving every thousandth tree, and the first 700 trees were discarded. In addition, phylogenetic analyses were performed by Maximum Likelihood algorithms, using the program MEGA-X (77) with standard settings and 1,000 bootstrap replications.