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.