RESULTS
Comparisons of sequencing and bioinformatic processing The ITS1F dataset included 146 samples (105 CA, 42 FL), had the most reads (14.2 million), and resulted in 7,609 OTUs. The LR22F dataset included 142 samples (106 CA, 36 FL), had the fewest reads (7.4 million), and resulted in 10,028 OTUs. The LROR dataset included 144 samples (107 CA, 37 FL), had an intermediate number of reads (9.1 million) and resulted in 5,786 OTUs (FIG 1). The ITS1F dataset had the highest percentage of reads discarded due to primer incompatibility (6.37%) while the LR22F set had the lowest percentage of primer incompatibility (0.38%), but the highest percentage of reads discarded due to short length (5.76%). Data from the LR22F primer set were mostly assembled into contigs that used both forward and reverse reads (FIG 1). In contrast, pairing between forward and reverse reads for ITS1F data was widely variable, whereas reads from the LROR primer set generally could not be compiled into contigs. This resulted in the majority of LROR OTUs being comprised of only forward reads (283 bp). Rarefaction plots indicate sufficient sampling from each state for each primer set (SUPP FIG 3).
The proportion of fungal OTUs was different for each of the three primer pairs and the two taxonomy assignment methods (FIG 2, SUPP FIG 4). The ITS1F dataset recovered the most fungal OTUs followed by LROR and LR22F according to the hybrid (for ITS1F) and UTAX (for LSU) methods of taxonomy assignment. The RDP classifier assigned 100% of ITS1F OTUs to Fungi but assigned a greater proportion of LROR and LR22F OTUs to non-fungal groups than the UTAX method (47.2-54% RDP vs 38.8-33.2% UTAX) (FIG 2 C). The majority of fungal reads from all three datasets were assigned to Ascomycota followed by Basidiomycota but the three datasets differed in the proportion of sequences that were assigned to EDF groups (FIGS 2, 4). The maximum OTU length was similar for all primer sets (ITS1F 534 bp, LR22F 548 bp, LROR 546 bp). Correlation plots of OTU length versus read number indicate that these variables are significantly negatively related for ITS1F, but not for the LSU datasets (SUPP FIG 5). Finally, the ITS1F dataset had the lowest percentage of unidentified OTUs with 1.9%, followed by LROR with 23.8%, and LR22F with 35.5% using the hybrid and UTAX methods, respectively. The RDP classifier recovered fewer unidentified OTUs with 0% for ITS1F, 18.3% for LR22F, and 24.1% for LROR.
Within Fungi, several orders of Dikarya were dominant across all primer sets and samples (SUPP FIG 6; SUPP TABLE 3). The greatest differences between primer sets were among less OTU-rich groups that were recovered by only one or two of the three primer sets (FIG 3; SUPP FIG 6). The LROR and LR22F datasets recovered more EDF OTUs from Blastocladiales, Calcarisporiellales, Chytridiales, Cladochytriales, Endogonales, Entomophthorales, Gromochytriales, Monoblepharidales, Neocallimastigales, Zoopagales, and Microsporidia than ITS1F (FIG 3). The LROR and LR22F datasets also recovered more OTUs from fungal-like organisms in class Oomycetes (kingdom SAR) and slime molds in Physariida (kingdom Amoebozoa). In other cases, the ITS1F dataset had several times the number of OTUs than either of the LSU datasets. For example, Rozellomycota was the dominant taxon among EDF for ITS1F, but Chytridiomycota was dominant for both LSU markers (FIG 3). Likewise, the ITS1F dataset had 36 OTUs assigned to Kickxellales compared to five or fewer for each of the LSU primer sets. This mirrors the inflation of Kickxellales in the mock community in the ITS1F dataset (see below).
Comparison of mock communities Inspection of the ITS1F mock community shows some likely errors in the OTU taxonomic assignment (TABLE 1). The AMPtk hybrid taxonomy assignment method identified two ITS1F OTUs as Stramenopiles (SAR), but each had BLAST matches to Acaulopage (Zoopagales), a mock community member (73% coverage, 91.5% identity and 72% coverage, 78.08% identity). The OTU identified as Dimargaris cristilligena had a BLAST match of 100% coverage and 97.4% identity to Cokeromyces recurvatus,the host of this mycoparasite. Five other OTUs found in the ITS1F Zoopagomycota mock community were identified as Basidiomycota, Metazoa, or Haptista. These Dikarya and non-fungal OTUs were not found in the negative controls and they were not returned in the mock communities sequenced with the LSU primers. These OTUs could have been amplified from the mixed genomic DNA present in the non-axenic mock isolates. The ITS1F hybrid dataset also returned 12 OTUs assigned to Kickxellales in the mock community whereas only nine taxa were actually included. Similarly, the ITS1F dataset had only three Zoopagales OTUs even though 21 isolates were originally added. The RDP classifier method identified 5 putative Mucoromycota host fungi OTUs and one Kickxellales OTU from the mock community. The remaining 21 OTUs were classified only to kingdom Fungi by RDP.
Mock community recovery was more accurate with the LSU datasets using the UTAX taxonomy assignment and RDP+SILVA database than the ITS1F dataset (TABLE 1). Both LSU primer sets recovered almost all the members of the mock community except that the LR22F dataset identified one lessCoemansia OTU. Both LSU primer sets also recovered OTUs assigned to putative Mucoromycota host fungi of the mycoparasites included in the mock. Comparison of UTAX against the RDP+SILVA LSU database to the RDP classifier shows that few of the mock community members were identified by RDP (TABLE 1). Many mock isolates were classified as Metazoa or remained unclassified with RDP but were accurately identified as fungi by UTAX and the RDP+SILVA LSU database. Across all primer sets, the Zoopagales mock members were underrepresented with several isolates remaining undetected by all three primer sets. There are multiple possible reasons these taxa did not amplify or sequence well: 1) we found that Acaulopage dichotoma , A. tetraceros , andStylopage species lack part of the ITS2 priming site (see alignment files at https://osf.io/cz3mh/), 2) amplification competition from shorter host DNA fragments present in the community, and/or 3) these fungi may have uncharacterized sequence features that reduce amplification, such as high G:C content (Dutton et al., 1993).
Ordination plots and statistical analyses – Community analyses are based on the subset of OTUs identified as Fungi or the subset of OTUs identified as EDF. The ITS1F primer set recovered 6,140 fungal OTUs, 1,757 of which were EDF, LROR recovered 2,095 fungi (369 EDF), and LR22F recovered 3,126 fungi (447 EDF). There were no observable differences in the fungal community recovered between primer sets based on the grouping of samples in the NMDS plots (FIG 4) (individual plots for each primer set separately are given in SUPP FIG 7). However, fungal communities from CA and FL were distinct with additional partitioning based on sampling site. The Evey Canyon and San Jacinto sites clustered separately from the Mojave Desert Sweeney sites in CA, whereas the two FL sites overlapped (FIG 4B, D). For the EDF, CA site clusters were separated in the ITS1F dataset, but those communities overlapped in the LSU datasets (FIG 4 A). PERMANOVA analyses found significant effects of site and sample type on the fungal communities among each of the primer sets (TABLE 2), but the effect size (R2) was small. Generally, sample type had a greater effect on EDF and fungal communities from FL, whereas site had a stronger impact on communities in CA. The ANOVA results of theB sim matrices for CA were all significant except for ITS1F site (fungi) and sample (EDF). Mantel test comparisons between primer sets were all significant for both fungi and EDF, indicating that the recovered communities were significantly correlated between the three datasets (TABLE 3). The correlation (R) was less than 0.50 for fungal ITS1F vs. LROR and LROR vs. LR22F and less than 0.60 for all EDF comparisons. The highest Mantel R statistic was for fungal ITS1F vs. LR22F, with a value of 0.71.
Plots of fungal richness between each primer set, site, and sample type showed that overall fungal richness was similar across all three primer sets, but soil samples generally had the highest diversity (SUPP FIG 8). Evey Canyon and San Jacinto had among the highest OTU richness for all sample types across all three primers. Fungal richness was lower for invertebrate samples across all primer sets but was higher for ITS1F than either LSU primer set. However, EDF diversity patterns differed. For ITS1F, many water, mud, and invertebrate samples had greater EDF diversity than soil (FIG 5). Conversely, soil and mud from Sweeney sites in the Mojave Desert had among the highest EDF diversity for both LSU markers (FIG 5). Tukey’s test results varied for EDF and all fungi and by primer (SUPP TABLE 4) and separation of EDF communities by sample type and site are also observable in the NMDS plots (SUPP FIG 7).
Phylogenetic reconstruction of LSU OTUs The 50 LROR and 50 LR22F OTUs identified only as “Fungi” were added to the sequence alignment that included 436 taxa. After exclusion of ambiguous sites, the LSU alignment contained 1,104 characters and OTUs had a final length of 216-239 bp for LR22F and 186-230 bp for LROR. Figure 7 shows the phylogeny and SUPP TABLE 5 lists the LSU OTUs used in the phylogeny and BLAST results for each OTU. Backbone nodes of the phylogeny were mostly unsupported whereas nodes near the tips had higher bootstrap support (>70). Both primer sets recovered monophyletic clades of OTUs that had BLAST matches to protist sequences (SUPP TABLE 5), but the LR22F dataset had fewer than LROR (FIG 6). The majority of these LROR OTUs had higher BLAST identity scores to the protist sequences than the LR22F OTUs, a maximum of 81% coverage and 100% identity for LROR, but only 20% coverage and 98.86% identity for LR22F. The LROR OTUs that had BLAST matches to protists were resolved in two different clades, one with three OTUs that had matches to Stramenopiles, and a larger clade with matches to Rhizaria. The LR22F protist clade was nested within the Chytridiomycetes, but all the matches (except the one mentioned above) had identity scores in the 75-78% range. Most fungal OTUs from both primer sets were resolved in the Chytridiomycetes and Orbiliomycetes. The OTUs placed in the Orbiliomycetes had close BLAST matches to Orbiliomycetes. In contrast, many OTUs placed in the Chytridiomycetes had matches to other fungal orders. Only LROR OTUs were placed within Aphelidiomycetes, and clades in the Basidiomycota, Eurotiomyctes, Umbelopsidomycetes, and Zoopagomycetes. The Archaeorhizomycetes, Endogonomycetes 1, and Sordariomycetes contained LR22F OTUs but none from LROR. The Glomeromycetes contained four LROR OTUs and one LR22F OTU. In the Zoopagomycota-only phylogeny, the two subphyla are recovered as polyphyletic, contrary to other studies (Davis et al., 2019), with the Dimargaritales and Ramicandelaber nested within the Zoopagomycotina (FIG 7). Additionally, the Harpellales are sister to the two subphyla rather than nested within Kickxellomycotina as found by other studies (Wang et al., 2019). Only OTU 6654 (LR22F) did not place in a clade matching its taxonomic classification from the taxonomy pipeline.