Total genomic DNA was extracted from axenic fungal cultures following Arnold and Lutzoni (2007) or with the Extract-N-Amp Plant Mini Kit (Sigma-Aldrich, St. Louis, MO, USA) following the manufacturer's instructions. The former is a phenol:chloroform based extraction and the latter a more rapid two-step kit that was adopted for convenience. Both methods provide comparable results (Sandberg et al., 2014). We used the forward primers ITS1F or ITS5 and reverse primers LR3 or ITS4 (10 μM) to amplify the nuclear ribosomal internal transcribed spacers and the 5.8S gene (ITS ribosomal DNA [rDNA]), and when possible, the first 600 base pairs (bp) of the large subunit (partial LSU rDNA). Fragment sizes ranged from ca. 600 (ITS rDNA) to 1200 bp (ITS rDNA–partial LSU rDNA). PCR methods followed Arnold and Lutzoni (2007) unless genomic DNA was obtained using the Extract-N-Amp kit, in which case methods followed the manufacturer's instructions.

Reactions were run on PTC-200 thermal cyclers (Bio-Rad, Hercules, CA, USA) with the following cycling parameters: 94°C for 3 min; 35 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 2 min; and 72°C for 10 min. PCR products were evaluated by visualizing with SYBR Green I (Molecular Probes, Invitrogen, Carlsbad CA, USA) after electrophoresis on 1% agarose gels. Positive PCR products producing single bands were sequenced directly (see below). Negative controls included water instead of template and were always blank, indicating no contamination during the DNA extraction or amplification process.

Products producing multiple bands or displaying weak amplification were cloned with StrataClone cloning kits (Agilent Technologies, Santa Clara, CA, USA) following the manufacturer's instructions, except that we used one-half the recommended reagent volumes for each reaction. Following blue/white screening, successful transformants were transferred to fresh plates and allowed to grow an additional 24 h to increase colony size. Eight positive clones per isolate were amplified using primers M13F(−20) and M13R(−27; numbers indicate specific primer variants) (10 μM) (Hoffman and Arnold, 2010). The cycling parameters were 94°C for 10 min; 35 cycles of 94°C for 1 min, 54°C for 1 min, and 72°C for 2 min; and 72°C for 10 min. Up to eight positive PCR products per isolate were selected for sequencing.

Positive PCR products showing one band of the appropriate length were cleaned with ExoSAP-IT (Affymetrix, Santa Clara, CA, USA) following the manufacturer's instructions. Cleaned products were diluted 1:2 with sterile water prior to sequencing.

Positive products were sequenced bidirectionally with the original PCR primers (5 μM) at the University of Arizona Genetics Core (Applied Biosystems BigDye Terminator v3.1 cycle sequencing kits; Applied Biosystems AB3730XL DNA Analyzer; Foster City, CA, USA). An assembly pipeline consisting of phred and phrap (Ewing and Green, 1998; Ewing et al., 1998) driven by Chromaseq (Maddison and Maddison, 2005) in Mesquite v.2.75 (Maddison and Maddison, 2009) was used to call bases and assemble reads into contigs. Base calls were verified by manual inspection of chromatograms in Sequencher v.5.1 (Gene Codes Corp., Ann Arbor, MI, USA). ITS rDNA–partial LSU rDNA sequences were submitted to GenBank (accessions KU977534-KU978121, KU978403-KU978448, and KX530755-KX530763).

We used BLASTn comparisons with GenBank (nucleotide collection [nr/nt]) (Altschul et al., 1990) and the Ribosomal Database Project (RDP) Fungal ITS Classifier (Deshpande et al., 2015) to select 223 strains for further study. These were isolates that were tentatively identified as Fusarium, Xylaria, or close relatives of those genera. Together these fungi represented isolates from 29 species of vascular plants from BCI and LS (three species of ferns, foliage only; 26 woody angiosperms, foliage and/or seeds) (Supplementary Table 1).

Top matches in database searches are not always closest phylogenetic neighbors and consequently may be inappropriate for identification (see U'Ren et al., 2009; Koski and Golding, 2012; Porter and Golding, 2012; U'Ren et al., 2016). Therefore, we used phylogenetic analyses to clarify the taxonomic placement of these strains and the resulting phylogenies to inform the distribution of seed-associated vs. foliar endophytic habits in each focal group.

BLASTn- and RDP classification with parallel examination of preliminary sequence alignments led us to assemble 19 data sets that encompassed Fusarium and close relatives in the Hypocreales, and Xylaria and close relatives in the Xylariales (Supplementary Table 1). Care was taken to match the phylogenetic information provided by ITS rDNA or ITS rDNA–partial LSU rDNA with the appropriate breadth of taxon sampling in each focal data set. We delimited our taxon sampling in part by reviewing existing phylogenetic hypotheses for each group, and by evaluating the limitations imposed by ambiguous alignments of ITS rDNA at broad taxonomic levels (see Schroers, 2001; Zhang et al., 2006; O'Donnell et al., 2008; Hsieh et al., 2010; Lombard et al., 2010; Chaverri et al., 2011; Short et al., 2013; U'Ren et al., 2016).

For each dataset, we first used Muscle v3.7 (Edgar, 2004) to align sequences with well-curated references from GenBank (i.e., those from type specimens and/or published species descriptions, and with information on host, substrate, and geographic origin, if available; Supplementary Table 2). We then edited each alignment by eye in Mesquite v2.75 (Maddison and Maddison, 2009). Ambiguously aligned characters were excluded. Matrix details are listed in Supplementary Table 3 and alignments submitted were to TreeBASE².

Maximum likelihood analyses were conducted in RAxML v.8.2.4 (Stamatakis, 2014) with bootstrap support determined on the basis of 1000 pseudoreplicates. Bayesian analyses were conducted in MrBayes v.3.2.6 (Huelsenbeck and Ronquist, 2001; Ronquist and Huelsenbeck, 2003) with six million generations, initiated with random trees, four chains, sampling every 1000th tree, and a burn-in of all trees with a standard deviation of split frequencies ≥ 0.01. Bayesian analyses did not converge after six million generations for the Fusarium solani species complex, so the run was extended to 30 million generations. For both ML and Bayesian analyses, we used jModelTest2 (Posada, 2008; Darriba et al., 2012) to select the appropriate model of sequence evolution. Models for each alignment are listed in Supplementary Table 3. Muscle, RAxML, Mr. Bayes, and jModelTest2 were implemented using the CIPRES Science Gateway (Miller et al., 2010).

Relationships remained unresolved in the analyses outlined above for two taxonomic groups: the Fusarium solani species complex (FSSC) and Gliocladiopsis. We therefore used multi-locus data to improve our inferences. For FSSC we used ITS rDNA–partial LSU rDNA, the RNA polymerase II second largest subunit gene (RPB2; 1500–1800 bp in length), and the translation elongation factor 1-alpha gene (TEF; 600–900 bp in length) following Short et al. (2014). For Gliocladiopsis we used ITS rDNA, the RNA polymerase II largest subunit gene (RPB1; 400 to 600 bp in length), and TEF following Castlebury et al. (2004) and Hirooka et al. (2013). PCR primers and cycling conditions are detailed in the Supplementary Materials and Methods. These sequences were submitted to GenBank (accessions KU978449-KU978453 [RPB1], KU978454-KU978547 [RPB2], KU978548-KU978605 [TEF]).

For each taxon, we aligned sequences with references, separately assembled and manually edited alignments, and conducted preliminary ML analyses for each gene as above. We explored topological congruence among the best ML gene trees for each taxon with the Wilcoxon signed-rank Templeton test (WT), implemented in PAUP^*v.4.0a150 (Swofford, 2002). Following the exclusion of 13 isolates due to missing data (see Supplementary Table 1), we detected conflict among single-locus datasets for the FSSC (WT, p ≤ 0.05 for 2/3 comparisons). However, because the multilocus framework adopted here has been used to classify phylogenetically diverse and clinically relevant fusaria and members of the FSSC in recent studies, our dataset includes many of these previously described isolates, and ITS data did not yield well-resolved or well-supported topologies in this group, we proceeded with using all of the data in a combined analysis (see Chang et al., 2006; Zhang et al., 2006; O'Donnell et al., 2008; Short et al., 2013). For Gliocladiopsis, we included all isolates and detected no conflict among the three single-locus datasets (WT, p > 0.1 for all comparisons). Single-locus alignments and trees for the FSSC and Gliocladiopsis were submitted to TreeBASE2.

Following assessment of congruence of individual gene trees, we concatenated single-locus alignments for each taxon into a supermatrix, and carried out phylogenetic reconstructions as above while allowing for independent evolution of individual loci by specifying partitions in each analysis. Bayesian analyses did not converge after six million generations for the FSSC, so the run was extended to 30 million generations. Matrix details are listed in Supplementary Table 3, and final concatenated alignments were submitted to TreeBASE².

Together these analyses clarified the phylogenetic placement of seed-associated fungi and foliar endophytic fungi, and showed that these functional groups have a shared evolutionary history (see below). Using this information we then structured our screening of EHB to evaluate whether bacterial endosymbionts differed in prevalence, diversity, or composition in seed-associated vs. foliar endophytic strains.

DNA extracted directly from apparently axenic fungal cultures (above) was screened for bacterial 16S rRNA through PCR using the primers 27F and 1492R following Hoffman and Arnold (2010). PCR products were visualized using gel electrophoresis as above. Positive PCR products producing single bands were prepared and sequenced directly as above. Positive products producing multiple bands or displaying weak or no amplification were cloned, and positive transformants were amplified and sequenced as above. All reads were processed as above. 16S rRNA gene sequences were submitted to GenBank (accessions KU978122-KU978402).

Fungal cultures that were positive for bacterial 16S rRNA were examined by microscopy with the LIVE/DEAD BacLight™ Bacterial Viability Kit (Invitrogen, Carlsbad, CA, USA) following Hoffman and Arnold (2010). We examined exemplar fungi representing at least three unique bacterial 16S sequences for each of the major taxonomic groups of EHB recovered here (Supplementary Table 1).

We prepared fungal samples for visualization by removing a small piece of mycelium (≤ 2-mm²) from the growing edge of single colonies growing on 2% MEA. Fragments were aseptically transferred to glass slides containing 20 μL of 1:1:18 LIVE/DEAD stain (component A: component B: diH₂O), teased apart using sterile insect mounting needles (size 00; BioQuip, Rancho Dominguez, CA, USA), covered with a coverslip, and incubated in darkness for 15 min. After incubation, we washed the mycelium by pulling sterile distilled water through the slide mounts with bibulous paper, and sealed the slides with two coats of nail polish. We used a Leica DM400B compound microscope with a 100-W mercury arc lamp for fluorescent imaging. Samples were viewed at room temperature with a Chroma Technology 35,002 filter set (480-nm excitation/520-nm emission) and 100 × APO oil objective.

Visible fluorescence of nucleic acids distinct from fungal mitochondrial or nuclear DNA, combined with the absence of extrahyphal bacteria and successful amplification of the bacterial 16S rRNA gene from fungal genomic DNA, provided evidence for the presence of viable EHB (Hoffman and Arnold, 2010; Arendt et al., 2016). We next evaluated the richness, phylogenetic placement, and composition of these EHB.

We used Sequencher v.5.1 (Gene Codes Corp., Ann Arbor, MI) to cluster EHB sequences at ≥97% sequence similarity in order to estimate EHB operational taxonomic units (OTUs) and calculate richness (Stackebrandt and Göbel, 1994; Kembel et al., 2014). We treated each fungal isolate as a distinct sample unit (i.e., containing an EHB community), and used the R package vegan to calculate diversity of EHB OTUs as Fisher's alpha (R Core Team, 2015; Oksanen et al., 2016). OTUs were classified tentatively by taking into account BLASTn comparisons with GenBank (16S rRNA collection; Altschul et al., 1990), the RDP Bacterial 16S Classifier (Wang et al., 2007), and the SILVA Incremental Aligner (“Search and Classify”; Pruesse et al., 2012; Quast et al., 2013). We looked for agreement among the three methods and carefully considered cases of disagreement. Results were used to frame taxon sampling for subsequent phylogenetic analyses.

We optimized (1) taxon sampling and (2) the taxonomic breadth appropriate for obtaining maximum resolution by assessing preliminary reconstructions at the level of bacterial phylum (or at the class level for Proteobacteria) (data not shown), and by considering published phylogenetic hypotheses for focal groups (Ludwig and Klenk, 2001; Vaneechoutte et al., 2004; Kämpfer et al., 2006; Tomitani et al., 2006; Nhung et al., 2007; Williams et al., 2007, 2010; Berrendero et al., 2011; Brady et al., 2013; Shih et al., 2013; Verma et al., 2013; Spring et al., 2015; Tank and Bryant, 2015). Overall, we generated 22 alignments that together contained 284 sequences representing EHB from our surveys. Our data sets also included 61 EHB sequences identified previously from endophytic Ascomycota (Hoffman and Arnold, 2010), 28 EHB sequences from root-associated Mucoromycotina, Mortierellomycotina, and Glomeromycota (Bianciotto et al., 2003; Partida-Martinez and Hertweck, 2005; Partida-Martinez et al., 2007a; Naumann et al., 2010; Sato et al., 2010; Lackner et al., 2011; Desirò et al., 2014, 2015), and 470 non-EHB reference sequences from GenBank representing the type or otherwise named strains with greatest homology to our unclassified EHB (Supplementary Table 4).

For each of the 22 datasets, we aligned 16S rRNA sequences with selected references from GenBank and inferred phylogenetic relationships as described above. Bayesian analyses did not converge after six million generations for Enterobacteriaceae, Pseudomonadales, and Pasteurellaceae (Gammaproteobacteria) such that these runs were extended to 20, 10, and 10 million generations, respectively. Matrix details are listed in Supplementary Table 3, and alignments were submitted to TreeBASE².

Treating each fungal isolate screened as a distinct sample, we evaluated whether communities of EHB differed as a function of the habit of the fungal strains from which they were isolated (i.e., seed-associated vs. foliar endophytic). We used Bray-Curtis dissimilarity, which was calculated using the relative abundances of all non-singleton EHB OTUs, and visualized results using non-metric multidimensional scaling (NMS, 999 runs). After finding a stable solution, we determined the proportion of variance explained by fungal functional group (seed-associated vs. foliar endophytic) with a post-hoc goodness-of-fit test and an independent analysis of similarity (ANOSIM, 999 permutations). Distance matrices, NMS, and ANOSIM were implemented using the R package vegan (R Core Team, 2015; Oksanen et al., 2016).

We used indicator species analysis (ISA) to explore whether certain taxa of EHB were strongly associated with the seed- or leaf-associated habit of host fungi. We first coded all non-singleton EHB OTUs by bacterial taxonomy (order or below except for the Bacteroidetes [phylum]) and the fungi in which they occur by habit (seed-associated or foliar endophytic), and then estimated the indicator value for each EHB taxon-host fungal habit combination. We assessed the significance of each indicator value by comparing with mean values obtained from a randomization test (999 permutations). Estimation of ISA parameters and indicator value significance were conducted using the R package indicspecies (De Cáceres and Legendre, 2009; R Core Team, 2015).

We annotated tips on bacterial phylogenies with data from original collections, GenBank, and/or published work to visualize patterns regarding fungal habit (i.e., seed-associated vs. foliar endophytic), fungal taxonomy (genus level or below), host plant species, and geography (Supplementary Tables 1, 2, 4, 5). We separately implemented multinomial logistic regression in JMP v.12.2 (SAS Institute, Cary, NC), combining those variables as predictors in an additive model to predict bacterial taxonomy at the level of phylum.

We measured congruence among the phylogenies of EHB, their fungal hosts, and the host plants from which the EHB-fungal associations were isolated. To obtain the host plant phylogeny, we used the Phylomatic tool implemented in Phylocom v.4.2 (Webb et al., 2008) to trim the Angiosperm Phylogeny Group megatree (APG III., 2009) to include only those plant lineages examined here. Cophylogeny between host plants and fungi or bacteria was assessed manually for all pairs of trees.

We used metadata for EHB in each bacterial phylogeny to compare phylogenetic diversity (PD) based on the habits (occurring in seeds vs. leaves) and geographic origins (temperate vs. tropical) of host fungi. We further compared EHB detected in this study to those known from root-associated Mucoromycotina, Mortierellomycotina, and Glomeromycota (MMG), as the latter represent ecologically and evolutionarily distinct groups in which EHB have been well-documented (Bianciotto et al., 2004; Partida-Martinez et al., 2007a; Naumann et al., 2010; Sato et al., 2010; Desirò et al., 2015). We calculated PD for each analysis as the sum of branch lengths spanning the minimum path connecting all taxa belonging to defined groups (Vane-Wright et al., 1991; Faith, 1992; Humphries et al., 1995; Faith and Baker, 2006; Arnold et al., 2007). Because PD typically increases with sampling effort and comparisons can be misinterpreted if sample sizes of groups are not standardized, we reduced the number of sequences in each focal group to that of the group with the smallest sample size prior to PD calculations (i.e., rarefied PD; Nipperess and Matsen, 2013). We randomly selected sequences for removal during rarefaction, and performed this process three times for each comparison to obtain mean PD and standard deviation for each group. We used the R package picante for trimming trees and inferring PD (Kembel et al., 2010; R Core Team, 2015). Rarefied PD was compared between groups using Welch two-sample t-tests.

To quantify the degree to which EHB among fungi of a particular habit form phylogenetically distinct groups, we grouped EHB in each bacterial phylogeny as for PD above and calculated phylogenetic signal (Blomberg et al., 2003). Although phylogenetic signal in continuous traits can be quantified in many ways (see Gilbert and Webb, 2007; Hardy and Pavoine, 2012; Münkemüller et al., 2012), few methods are available for estimating phylogenetic signal in binary traits. We quantified phylogenetic signal using the sum of sister-clade differences (character dispersion, D, where D ≥ 1 indicates weak signal or convergent evolution and D ≤ 0 indicates strong signal; Fritz and Purvis, 2010). Strong signal implies relative phylogenetic clumping, or monophyly, and would indicate phylogenetically distinct groups. For each bacterial phylogeny, we estimated D for each group, and tested values for significant departure from a model in which traits have a phylogenetically random distribution (D = 1; 1000 permutations), and a model in which traits are clumped as if evolved by Brownian motion (D = 0; 1000 simulated walks). Analyses were conducted using the R package caper (Orm et al., 2013; R Core Team, 2015).