Global Diversity and Biogeography of the Zostera marina Mycobiome

Abstract

Seagrasses are marine flowering plants that provide critical ecosystem services in coastal environments worldwide. Marine fungi are often overlooked in microbiome and seagrass studies, despite terrestrial fungi having critical functional roles as decomposers, pathogens, or endophytes in global ecosystems. Here, we characterize the distribution of fungi associated with the seagrass Zostera marina, using leaves, roots, and rhizosphere sediment from 16 locations across its full biogeographic range. Using high-throughput sequencing of the ribosomal internal transcribed spacer (ITS) region and 18S rRNA gene, we first measured fungal community composition and diversity. We then tested hypotheses of neutral community assembly theory and the degree to which deviations suggested that amplicon sequence variants (ASVs) were plant selected or dispersal limited. Finally, we identified a core mycobiome and investigated the global distribution of differentially abundant ASVs. We found that the fungal community is significantly different between sites and that the leaf mycobiome follows a weak but significant pattern of distance decay in the Pacific Ocean. Generally, there was evidence for both deterministic and stochastic factors contributing to community assembly of the mycobiome, with most taxa assembling through stochastic processes. The Z. marina core leaf and root mycobiomes were dominated by unclassified Sordariomycetes spp., unclassified Chytridiomycota lineages (including Lobulomycetaceae spp.), unclassified Capnodiales spp., and Saccharomyces sp. It is clear from the many unclassified fungal ASVs and fungal functional guilds that knowledge of marine fungi is still rudimentary. Further studies characterizing seagrass-associated fungi are needed to understand the roles of these microorganisms generally and when associated with seagrasses. IMPORTANCE Fungi have important functional roles when associated with land plants, yet very little is known about the roles of fungi associated with marine plants, like seagrasses. In this study, we report the results of a global effort to characterize the fungi associated with the seagrass Zostera marina across its full biogeographic range. Although we defined a putative global core fungal community, it is apparent from the many fungal sequences and predicted functional guilds that had no matches to existing databases that general knowledge of seagrass-associated fungi and marine fungi is lacking. This work serves as an important foundational step toward future work investigating the functional ramifications of fungi in the marine ecosystem.

Keywords: 18S rRNA; ITS2; Zostera marina; abundance-occupancy; core; dispersal limited; eelgrass; global distribution; marine fungi; microbial eukaryotes; mycobiome; plant selected; seagrasses.

FIG 1

Within-sample diversity varies across tissues and sites. Here we depict sample Shannon diversities for the ITS2 region for each sample type (leaf, root, sediment) using box plots (A) and within each sample type at each collection site (see Table 1 for site codes) using bar charts (B). Shannon diversities for the 18S rRNA gene amplicon data are also plotted for each sample type (leaf, root, sediment) using box plots (C) and within each sample type at each collection site (see Table 1 for site codes) using bar charts (D). For panels B and D, the standard error of the mean Shannon diversity at each site for each sample type is represented by an error bar, and bars are colored by sample type. Post hoc Dunn results for pairwise site comparisons for each sample type can be found in Tables S1 to S6 in the supplemental material. For panels A to D, zero values indicate samples that are dominated by a single taxon.

FIG 2

Community structure varies between tissues and ocean basins. Principal-coordinate analysis (PCoA) visualization of Hellinger distances of fungal communities associated with leaves, roots, and sediment based on ITS2 region amplicon data (A and C) and 18S rRNA gene amplicon data (B and D). (A and B) Points in the ordination are colored and represented by shapes based on sample type: leaf, yellow circles; root, green triangles; sediment, blue squares. (C and D) Points in the ordination are colored by ocean (Pacific, blue; Atlantic, orange) and represented by shapes based on sample type (circles), root (triangles), or sediment (squares).

FIG 3

Mantel tests suggest a distance decay relationship. Scatterplots depict the weak but significant positive distance-decay relationship between leaf fungal community beta diversity (Hellinger distance) using the ITS2 region amplicon data and geographical distance (km) between sites from the Pacific Ocean (r = 0.1767, P = 0.0001) (A) and Atlantic Ocean (r = 0.1057, P = 0.0001) (B).

FIG 4

Overlap between predicted core mycobiomes of individual Z. marina tissues. Venn diagrams representing shared core ASVs as defined by abundance-occupancy distributions for each sample type (leaf, root, sediment) for ITS2 region amplicon data (A) and 18S rRNA gene amplicon data (B). For comparison, overlap of the entire mycobiome of individual Z. marina tissues is shown in Fig. S5 in the supplemental material.

FIG 5

Abundance-occupancy distributions reveal core mycobiomes. Abundance-occupancy distributions were used to define core members of the leaf (A), root (B), and sediment (C) mycobiomes for the ITS2 region amplicon data. Each point represents an ASV, with predicted core members indicated by a color (leaf = yellow, root = green, sediment = blue) and noncore ASVs in white. Ranked ASVs were predicted to be in the core based on a final percent increase equal to or greater than 10%. A solid line represents the fit of the neutral model, and dashed lines represent 95% confidence around the model prediction. ASVs above the neutral model are predicted to be selected for by the environment (e.g., by the host plant, Z. marina), and those below the model are predicted to be selected against or dispersal limited.