Utility of large subunit for environmental sequencing of arbuscular mycorrhizal fungi: a new reference database and pipeline

Abstract

Arbuscular mycorrhizal fungi (AMF – phylum Glomeromycota) form symbioses with most plant species worldwide and play critical roles in plant nutrient and water uptake, pathogen resistance and soil aggregation (Smith & Read, 2008; Delavaux et al., 2017; Brundrett & Tedersoo, 2018). Because AMF community composition influences ecological function (van der Heijden et al., 1998; Vogelsang et al., 2006; Koziol et al., 2018), understanding patterns of AMF composition is a research priority. Hyphae of AMF species are not morphologically distinguishable, and therefore quantification of AMF species diversity and community composition has increasingly relied on metabarcoding of ribosomal RNA (rRNA) gene sequences from field samples (Öpik et al., 2014). However, to date, no single region of the rRNA gene has been universally accepted as optimal for AMF environmental sequencing. The internal transcribed spacer (ITS) region of the rRNA gene has been suggested as the universal fungal marker (Schoch et al., 2012; Lindahl et al., 2013) and has been used for AMF biogeographical studies (Tedersoo et al., 2014) and environmental sequencing (Öpik et al., 2014). However, this region is suboptimal as a marker gene for AMF (Stockinger et al., 2010; Schoch et al., 2012). The sequence matching approach used for ITS sequences with other fungi is of limited utility for AMF because of the poor representation and poor curation of AMF sequences in ITS sequence databases (Bidartondo, 2008; Stockinger et al., 2010). This database problem cannot be easily rectified because a high proportion of AMF encountered in environmental samples are undescribed. While phylogenetic approaches can be used to identify new sequences as AMF, this approach cannot be used for ITS amplicons because its rapid sequence evolution (Nilsson et al., 2008) does not generate reliable trees. The most commonly used region of the rRNA gene for environmental sequencing of AMF is the small subunit, or SSU (Öpik et al., 2014). The utility of this region is enhanced by a well-developed and curated database for AMF (Öpik et al., 2010; Davison et al., 2015). However, the SSU region has the disadvantage of being slow-evolving and therefore not sufficiently variable to adequately resolve AMF species (Krüger et al., 2009; Bruns & Taylor, 2016; Schlaeppi et al., 2016). By contrast, the large subunit (LSU) region consistently shows greater utility for taxonomic resolution for AMF (Krüger et al., 2012; Hart et al., 2015; House et al., 2016), making it potentially more useful in environmental AMF sequencing. Thus far, the LSU region has rarely been used in environmental sequencing of AMF (Gollotte et al., 2004; Lekberg et al., 2013; House & Bever, 2018; Vieira et al., 2018; Schütte et al., 2019), perhaps because of bioinformatical challenges in implementation. Here, we aim to expand the utility and ease the adoption of the LSU for amplicon sequencing of AMF by providing a well-curated LSU reference database, a reference backbone tree for phylogenetic placement and a computational pipeline easily implemented using current bioinformatical tools. We present a current curated database and reference tree using AMF sequences from several sources: a subset of sequences published by Krüger et al. (2012) available on the National Center for Biotechnology Information (NCBI), a database of unpublished spore-derived sequences from the International Culture Collection of (Vesicular) and Arbuscular Mycorrhizal Fungi (INVAM, Morgantown, WV, USA) and additional recently described sequences from NCBI. Sequences amplified by Krüger et al. (2012) generally used primers SSUmAf and LSUmAr, with a second amplification round with SSUmCf and LSUmBr or SSU-Glom1-NDL22, while sequences amplified by INVAM used primers 1TS1 and NDL22, followed by a second amplification round using primers LR1 and NDL22 (Morton & Msiska, 2010); additional primers used in recently described sequences from NCBI can be found in each respective publication associated with the accession numbers detailed in Fig. 1. To build the reference backbone tree, representative sequences were chosen to maintain the tree structure, conserving clear clades within the tree (using the Interactive Tree of Life to view the tree; Letunic & Bork, 2019), but were kept at a minimum, to make the use of the tree as a reference computationally feasible. These sequences were aligned using Mafft (Katoh & Standley, 2013) and a tree constructed using Raxml v.8 (Stamatakis, 2014) with 1000 bootstrap replicates and the evolutionary model Gtrgamma in Qiime2 (https://qiime2.org). Outgroups were Mortierella elongata (MH047197, Mucoromycota), Exophiala spinifera (MH876260; Basidiomycota) and Rhodotorula hordea (AY631901; Ascomycota). In addition, we included LSU sequences of a plant, Citrus limon (X05910, Rutaceae), and an animal, Rutilus rutilus (EF417167, Cyprinidae). We use the reference database as the reference for open (closed and then de novo) operational taxonomic unit (OTU) clustering and use a phylogenetic tree generated from these same sequences (Fig. 1; Supporting Information Fig. S1) as a backbone constraint in constructing phylogenies to place study sequences. All previously unreported reference sequences have been uploaded to Genbank (MT832155–MT832238). We present a pipeline starting from raw Illumina LSU sequences (MiSeq V3, 2x300 bp) and ending with an OTU table of phylogenetically defined putative AMF OTUs for downstream analyses (Fig. 2). This pipeline is built for use with a high-performance computing (HPC) cluster using a Simple Linux Utility for Resource Management (SLURM) workload manager. Our pipeline covers key bioinformatical steps, including primer removal, quality control, and OTU clustering (Vsearch algorithm (Rognes et al., 2016); for a discussion of why OTUs may be preferable over amplicon sequence variants (ASVs) for this particular application, see Methods S1. Importantly, our pipeline maintains nonoverlapping forward and reverse reads for each sequence, retaining all possible data in long, non-paired reads (700–900 bp; see Table S1 for a forward and reverse read concatenation test). In addition, the pipeline places representative OTU sequences within a tree using our backbone phylogeny and subsequently extracts those that fall within the AMF clade. Operationally, our software handles many batches of OTUs in parallel, thereby greatly improving processing speed. Finally, the output of all batches is joined into a single OTU table containing counts of putative AMF OTUs in each sample, along with a Fasta file containing the representative sequences of all phylogenetically determined putative AMF OTUs. The analogous data files for the OTUs that fall outside of the AMF clade are also provided by this final step. The full bioinformatical pipeline and description can be found in Methods S1; all required files are also supplied (Methods S2; https://github.com/c383d893/AMF-LSU-Database-and-Pipeline). In the age of metagenomics, it is attractive to declare one marker suitable for all fungi (Lekberg et al., 2018). Universal fungal primers would facilitate efforts to compare relative abundances between taxonomic and functional groups and to identify global scale biogeographic patterns. Nonetheless, here we confirm that the general primers targeting the ITS region are not adequate for detecting a majority of undescribed AMF species. We estimate that nearly 90% of phylogenetically defined putative AMF OTUs in our test dataset derived from a Midwestern US grassland are undescribed (i.e. have no Glomeromycota Blast match in NCBI) and c. 30% do not group with described families in our phylogenetic tree (Table S2). Given this severe limitation in building a database and current database constraints (Stockinger et al., 2010; Schoch et al., 2012; Hart et al., 2015), the sequence matching algorithms are not adequate for environmental sequencing of AMF regardless of rRNA gene region. The phylogenetic approach we use here can accommodate undescribed AMF taxa in environmental samples. While the ITS region evolves too quickly to allow reliable tree construction, the LSU region can be used to build a phylogenetic tree and place all study sequences inside or outside of the conserved AMF clade. This allows identification of undescribed putative AMF through clade placement of any environmental sequence. We illustrate this benefit in two studies in which analyses of LSU amplicons reveals significant environmental patterns in AMF composition in US grasslands and boreal forests that were not evident in analyses of ITS amplicons (Methods S3; Table S3). The SSU region, like the LSU region, can be used to build phylogenetic trees and place environmental sequences in the AMF clade (Öpik et al., 2014; Stefani et al., 2020). In addition, the SSU amplicons from commonly used primers have been easier to handle bioinformatically because the small length of c. 500 bp allows for merging of short, paired-end Illumina reads (Lee et al., 2008; Dumbrell et al., 2011). Second, environmental sequences can be directly assigned to virtual taxa in a large publicly available database (Öpik et al., 2010, 2014). However, the SSU region has the disadvantage of being slowly evolving, thereby limiting inferences to taxonomically coarse designations (Krüger et al., 2009; Stockinger et al., 2010; Schoch et al., 2012). Individual virtual taxa assigned using the MaarjAM database include many distinct species and have been suggested to be analogous to genera (Bruns & Taylor, 2016). Analogous problems are present in the LSU in that sequence variation within isolates can be attributed to different OTUs, and individual OTUs can include sequences of different species, particularly for the Claroideoglomeraceae (Stockinger et al., 2010; House et al., 2016). Nonetheless, the LSU does a much better job of capturing other AMF families (Krüger et al., 2012; House et al., 2016) and has been suggested as the most suitable region within short read length restrictions (Stockinger et al., 2010). How these differences translate into inferences on issues such as environmental dependence of AMF distributions or frequency of endemism remains to be evaluated. Here we provide a well-curated LSU reference database, a backbone phylogeny and a computational pipeline that uses these resources to process environmentally derived amplicon sequences. We are optimistic that this set of tools will facilitate molecular work with AMF within the LSU region, leading to finer scale assessments of ecological inferences from AMF community structure. The authors thank the Center for Resource Computing at the University of Kansas as well as RAxML and QIIME2 communities for their support. The authors acknowledge fruitful conversations with Benjamin Sikes about the manuscript pipeline, Bill Wheeler for assistance with INVAM culture meta-data, Rob Ramos for help with initial code development, friendly reviews by Ylva Lekberg and Maarja Öpik, and financial support from the National Science Foundation (DEB-1556664, DEB-1738041, OIA-1656006, IOS-2016351, DEB-0076066). SLS would like to thank the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) for a Research Assistantship (Process 307.995/2019-4) and Universidade Regional de Blumenau (FURB) for supporting a sabbatical leave at the University of Kansas. CSD and JDB designed the study. JBM generated the LSU sequences from INVAM isolates. CSD, JDB and SLS led the phylogenetic tree analysis. MRW and CSD developed the bioinformatical pipeline. CSD and JDB led interpretation and manuscript writing, with substantial contributions from all authors (CSD, SLS, MRW, US, JBM, JDB). The unpublished sequences that support the findings of this study are openly available at NCBI at https://www.ncbi.nlm.nih.gov/, accession numbers MT832155–MT832238. Data used to test the pipeline are uploaded under NCBI Project no. PRJNA648993. All other code is presented in Methods S1 and S2 and can be found at https://github.com/c383d893/AMF-LSU-Database-and-Pipeline. Fig. S1 Backbone tree with bootstrap support values. Methods S1 Bioinformatical pipeline. Methods S2 Pipeline scripts and instructions. Methods S3 Comparing ecological inference between ITS database and LSU phylogenetic results. Table S1 Forward and reverse read concatenation test. Table S2 Family classification of OTUs generated from our pipeline with test dataset. Table S3 Comparing ecological inference between ITS database and LSU phylogenetic results. Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office. Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.

Topics & Concepts