Phylogenomics of fescue grass-derived fungal endophytes based on selected nuclear genes and the mitochondrial gene complement

Endophyte isolates and DNA extraction

Paired-end library preparation and sequencing

Genomic DNA was fragmented in a Covaris instrument (Woburn, MA, USA) to an average size of 100–900 bp. For each endophyte DNA sample, paired-end libraries with inserts c. 400 bp in size were prepared using the standard protocol (TruSeq DNA Sample Prep V2 Low Throughput: Illumina Inc., San Diego, USA) with paired-end adaptors. Library quantification was performed using the KAPA library quantification kit (KAPA Biosystems, Boston, USA). Paired-end libraries were pooled according to the attached adaptors and sequence analysed using the HiSeq2000 platform (Illumina) following the standard manufacturer’s protocol.

Processing and assembly of sequence data

All generated sequence reads were quality controlled by filtering and trimming of reads based on quality using a custom Python script, which calculates quality statistics, and stores trimmed reads in several fastq files. Data assembly was performed using the Linux-based de novo assembler Velvet ver.1.1.06 [29]. For Velvet assembly, different hash lengths (K-mer sizes) ranging from 39 to 51 were tested as appropriate for different sequence read sets, and the minimum contig length was always defined as 200 bp. Values for estimated coverage and coverage cut-off were set to auto.

Assembly of nuclear gene sequences

Presence and copy number of tub2, tefA, and perA genes (using Genbank accession numbers:tub2: AY722412, tefA: FJ660614, perA: AB205145) in each endophyte genome were initially determined through nucleotide BLAST (Basic Local Alignment Search Tool) [30] analysis using contigs from the optimised Velvet assembly of total reads as the database. In order to assemble each gene copy, matching reads for each reference gene were identified from a database of all trimmed reads from each endophyte genome, through a similarity search using BLAST, defining the E value threshold as 0.1. From this BLAST output, all corresponding paired reads (both forward and reverse) were extracted from the database, and the second reads were reverse-complemented. Each of the first and second reads was concatenated and used as BLASTN queries against a database consisting solely of each reference gene sequence for assembly. From this search, reads that matched in an anti-sense orientation were reverse-complemented in order to orientate the concatenated reads in 5′- to 3′-orientation against the reference gene sequence. Subsequently, reads were separated into two distinct sets (designated read1 and read2) and the two groups were used individually as BLASTN queries against the gene sequence that was to be assembled. The aim of these individual BLAST searches was to estimate positions for each read against the reference database sequence. After assigning positions, each read was padded to the appropriate position along the reference gene sequence using a customised PERL script. The padded read pairs (when both reads were selected from the BLAST output) were then concatenated and saved in FASTA format. Using the graphical multiple sequence editor SeaView [31], the padded reads were manually assembled into the defining gene sequences. Functional sequences were predicted for multiple perA gene copies through translation of each gene sequence using ExPASy online DNA sequence translation tool [32].

Phylogenetic analysis of nuclear genes

Multiple alignments of the complete gene sequences for each selected gene were performed individually using ClustalW [33] with default parameters. To reconstruct tree topology, parsimony, maximum likelihood (ML) and neighbour-joining (NJ) methods were used as implemented in MEGA 5 [34] with default parameters and 1,000 bootstrap replicates. Gene sequences available in Genbank from related endophyte species were used appropriately, and corresponding accession numbers are provided in each tree diagram. After identification of proposed origin from an individual genome for each gene copy based on the individual gene trees, the three nuclear genes were concatenated in the order tub2-tefA-perA and aligned using ClustalW with the default parameters. For the concatenated multiple sequence alignment, phylogenetic topology was reconstructed using MEGA with default parameters and 1,000 bootstrap replicates. Phylogenetic networks were constructed for aligned concatenated gene sequences using the NeighborNet algorithm [35] on the Nei–Li pairwise distance matrix, and network diagrams were produced using the program SplitsTree4 [36].

Assembly of mitochondrial genomes

Contigs of mitochondrial (mt) genome origin were initially identified using nucleotide BLAST at an E value threshold of 0.001 through alignment of a database containing all contigs from the optimised Velvet assembly (as described above) against the mt genome of the N. lolii standard endophyte (SE) strain as a reference [37] (Genbank accession number KF906135). For each candidate endophyte, a set of contigs with higher read depth coverage were shown to have a significant match to the reference mt genome sequence. A cut-off value for read depth coverage of each mt genome was identified based on the BLAST output, and a second Velvet assembly was performed with assignment of this value as the coverage cut-off value in order to filter contigs derived from the mt genome. A range of k-mer values were tested, and a final assembly was accepted on the basis of features such as total number of assembled contigs, N50 value and cumulative contig length. For those mt genomes with few (2–5) contigs, ordering was performed with BLASTN (E value threshold of 0.001) based on the pre-existing SE mt genome sequence, and overlapping regions were manually linked. In order to confirm gaps observed in comparison to the SE mt genome, alignment was performed for trimmed reads using the Burrows-Wheeler Alignment (BWA) tool [38] with the maximum number of gap openings set to five. Mapped reads were viewed using Tablet 1.12.02.06 [39], a graphical viewer for sequence alignment. Observed gaps were further confirmed through grouping of observed gap positions within each endophyte species or taxon.

Identification of protein-coding gene sequences was performed using each mt genome sequence as the query database against the individual mitochondrial protein gene sequences from the clavicipitacean entomopathogenic fungus Metarhizium anisopliae (Genbank accession number NC008068). Identified protein-coding genes were concatenated according to the order observed in the M. anisopliae mt genome.

Phylogenetic analysis of mt genome protein coding gene complement

Multiple alignment of concatenated mitochondrial protein-coding gene complements from the 19 endophytes and counterparts in M. anisopliae was performed using the M-LAGAN program within the mVISTA on-line suite of computational tools [40], with default parameters. Alignments were manually edited for mis-alignments that may have accumulated due to overlapping gene fragments. To reconstruct the tree topology, parsimony, ML and NJ methods were used as implemented in MEGA 5 with default parameters and 1,000 bootstrap replicates. Furthermore, to study the level of identity of each mt genome protein-coding complement relative to that of M. anisopliae, aligned sequences were visualized through use of a VISTA plot [41].

Identification of individual nuclear gene copies

Presence of the three nuclear genes (tub2, tefA and perA) was determined for 13 fescue-derived endophyte genomes and all other reference taxa (Epichloë spp. and N. lolii) that were used for this study. The observed number of copies for each gene within each taxon was then determined (Table 2). For each gene, 3 copies were observed within individual N. coenophialum genomes, while 2 copies were observed from all other taxa apart from N. uncinatum, which contained a single copy of the tub2 gene. Moreover, one of the assembled perA gene copies from FaTG-2, FaTG-3 and UNS genomes revealed a common variant structure based on large- to moderate-sized deletions (coordinates 1251-1878 bp and 4590-4918 bp; Figure 1).

	Gene	Endophyte taxon
Observed number of gene copies for		N. coenophialum	FaTG-2	FaTG-3	FaTG-3 like	UNS	N. uncinatum
	perA	3	2	2	2	2	2
	tub2	3	2	2	2	2	1
	tefA	3	2	2	2	2	2

Phylogenetic relationships based on nuclear genes

Phylogenetic relationships were reconstructed based on individual gene sequences of perA, tub2, and tefA genes, as well as the concatenated sequences of all three genes, using parsimony ML and NJ methods (related alignments can be found at URL: http://purl.org/phylo/treebase/phylows/study/TB2:S14923). All positions containing gaps and missing data were eliminated during each analysis. Corresponding gene-specific DNA sequences from other taxa that were previously deposited in GenBank were included. Similar tree topologies were observed for all three methods, and the most parsimonious phylograms were selected (Figures 2, 3 and 4 and Additional files 1, 2, 3, 4 and 5). The GenBank accession numbers of the nuclear genes derived from fescue-derived endophytes are provided in Additional file 6. Maximum parsimony analysis based on the individual gene sequences of the three genes resolved more than one most parsimonious tree, and so bootstrap consensus trees inferred from 1000 replicates are displayed for each gene. In each phylogram, individuals from the same endophyte taxa were clustered together, and were separated from the Epichloë isolates that were used in the study. Apart from the placement of E. baconii and E. amarillans within the perA-and tefA-specific trees, phylogenetic analysis of individual genes resolved similar genomic relationships between endophyte taxa. All major clades that were defined by single gene phylogeny were also strongly supported by the phylogeny of the concatenated perA, tub2, and tefA genes, due to observation of similar tree topology. Maximum parsimony analysis of concatenated genes yielded a single most parsimonious tree with high level of bootstrap support for the majority of the individual branches.

In all instances, the phylogenetic trees were predominantly separated into two major groups representing different Epichloë species. For example, in the tefA gene-specific phylogram (Figure 3), Group 1 contained the taxa E. festucae, E. baconii, E. amarillans and E. bromicola, while Group 2 contained E. sylvatica, E. typhina, and E. clarkii. Furthermore, individual gene copies from the fescue-derived endophytes were located closely adjacent to the reference E. festucae-, E. bromicola- and E typhina-derived sequences, suggesting affinities to putative sub-genome components of the heteroploid taxa. However, in all instances members of the fescue-derived endophyte (N. coenophialum, FaTG-2, FaTG-3, FaTG-3-like and UNS) gene copy 1 (FGC1) clade were not so closely related to sequences from any of the currently included Epichloë endophytes (Figures 2, 3 and 4 and Additional file 2). In the perA-based phylogeny, E. amarillans formed a sister group to FGC1, while addition of a partial gene sequence from E. amarillans to the tefA-based phylogeny generated a separate clade also containing E. baconii. Nevertheless, in the absence of E. amarillans, the tub2-and concatenated gene-based analysis placed FGC1 as a sister group to E. baconii.

The genomes of the known heteroploid endophyte taxa N. coenophialum, N. uncinatum and FaTG-3 contributed gene copies to both major phylogram groups, consistent with hybrid origins either from one of the relevant Epichloë species, or a closely related taxon. In contrast, the multiple gene copies from FaTG-2 and UNS were located only within Group 1. High levels of similarity between sequences from different endophyte taxa suggested the presence of common sub-genomic components. For instance, N. coenophialum and N. uncinatum genotypes always contributed gene copies that were identical or very closely related to one another, and similar to those from E. typhina. Although N. coenophialum gene copies showed close affinities to those from other fescue endophyte taxa within Group 1, distinct clusters were generated in all instances. A further point of interest was the close relationships between the FaTG-2-derived perA and tefA gene copies and those obtained from both E. festucae and its asexual anamorph, N. lolii.

As an alternative way to explore the reticulated evolutionary relationships between sexual Epichloë progenitor species and heteroploid Neotyphodium endophytes, a phylogenetic network diagram was constructed based on the concatenated nuclear gene sequence (Figure 5). Distinct clusters for asexual endophyte-derived gene copies were formed around the reference E. festucae and E. typhina sequences, corresponding to respective clades in Group 1 and Group 2 of the phylograms, and further supporting the hypothesis of progenitor relationships. Separation of E.baconii from FGC1 gene copies was also consistent with the phylogram structure. The network analysis also served to further demonstrate the differentiation of FaTG-2 and UNS-derived gene copies within the E. festucae-containing clade of the phylograms, despite similarity in both instances to E. festucae. Similar results were obtained for FaTG-3- and FaTG-3-like-derived gene copies corresponding to members of FGC1 within the phylograms, and, to a lesser extent, Group 2.

Mitochondrial genome sequence structure

Phylogenetic relationships based on mitochondrial genome comparisons

Phylogenetic relationships were reconstructed based on the concatenated sequences of 13 mitochondrial protein-coding genes from the fescue endophytes, while the equivalent for M. anisopliae was used as the out-group (Figure 6). Similar tree topology was observed for parsimony, ML and NJ methods (related alignments can be found at URL: http://purl.org/phylo/treebase/phylows/study/TB2:S14923). As expected, individual sequences from the same taxon clustered together. A number of putative progenitor relationships, such as that between E. typhina and N uncinatum, were more readily apparent from the phylogram. Close relationships were revealed between the N. lolii and LpTG-2 mt genomes and that of their putative sexual progenitor, E. festucae, and similar, but less close, relationships were apparent for N. coenophialum and FaTG-2. Commonalities of mitochondrial genome structure were evident between the FaTG-3 and FaTG-3-like mt genomes, albeit with lower bootstrap support, but affinity to potential Epichloë genomes was not so obvious, although an E. festucae mt genome provides the most obvious candidate. This result was inconsistent with the data from the nuclear gene analyses in the present study, that revealed closer relationships of gene copies from the FaTG-3 and FaTG-3-like isolates to the putative FGC1 progenitor, and to E. typhina. A clear differentiation of the UNS mt genome from that of the preceding groups, with higher levels of sequence similarity to E. baconii and E. typhina rather than E. festucae, was evident from this analysis.

Copy number variation of selected nuclear genes

The presence of 3 copies for each of the tefA, tub2 and perA genes in the N. coenophialum genome suggests a tri-parental hybrid origin, consistent with previous studies [17, 23]. Similarly, the observation of 2 copies for each gene in the genomes of other heteroploid endophyte taxa (FaTG-2, FaTG-3, FaTG-3-like and UNS) is compatible with a series of bi-parental hybrid origins. Although N. uncinatum has also previously been inferred to have arisen as a bi-parental hybrid, the presence of a sole tub2 gene copy suggests selective gene loss, as previously proposed to account for the current heteroploid constitution of this and other taxa. These results, apart from concordance with earlier sequence-based studies [14, 19], are also consistent with the complexity of SSR profiles from the same accessions [9], which typically contained up to 3 distinct amplicons from N. coenophialum genotypes, and up to 2 amplicons from the other taxa in this study.

Phylogeny of previously described fescue-derived endophytes

The present study has also permitted identification of those Epichloë species that are likely to be most closely related to the taxa that participated in hybrid origins. Previous phylogenetic studies based on two of the nuclear genes used in this study (tub2 and tefA), as well as the act1 actin gene, have provided evidence for progenitor identity [15]. Three tall fescue-derived endophyte taxa have previously been included in such studies, N. coenophialum was proposed to have originated from E. festucae-, E. baconii- and E. typhina-like ancestors [18, 19, 42], while FaTG-2 and FaTG-3 were suggested to be derived from E. festucae- and E. baconii-like, and E. baconii- and E. typhina-like ancestors, respectively [17]. As summarised in Figure 7, the present study was consistent with these predictions in terms of affinities to contemporary E. festucae and E. typhina genotypes, but more distant relationships were observed for E. baconii. The group designated FGC1 in this study, which cannot be unequivocally attributed to a E. baconii-like progenitor, was also identified in a previous study and termed the ‘Lolium-associated endophyte’ (LAE) clade [23, 43]. Furthermore, two distinct E. typhina lineages appear to have contributed to formation of the N. coenophialum/N. uncinatum and FaTG-3/FaTG-3-like heteroploid genomes, respectively (Figure 7), based on interpretation of the tree and network diagrams.

Phylogenetic reconstruction based on the perA gene sequence revealed a closer relationship between E. amarillans and the FGC1 than for the E. baconii genotype that was used in this study. However, the E. amarillans-derived tefA gene sequence demonstrated a close genetic relationship to E. baconii, and E. amarillans formed a sub-clade with E. baconii, E. festucae, and N. lolii as well the FGC1 clade, consistent with previous studies [43]. Observed anomalies between the gene-specific phylogenies in the present study may be due to different rates of molecular evolution between endophyte-specific (perA) and housekeeping (tefA and tub2) genes. Further to this, addition of the entire tefA gene sequences of E. amarillans to the phylogenetic analysis may provide a higher level of resolution.

Phylogeny of mt genomes

For N. coenophialum, FaTG-2 and FaTG-3, the mitochondrial gene complement-based analysis revealed closest relationships to E. festucae, suggesting that this or a closely related sexual species donated the cytoplasmic genome to the known heteroploid taxa. This conclusion is again consistent with previous studies, apart from the status of FaTG-3, which does not show strong similarity to the mt genomes of either E. baconii or E. typhina. This anomaly may be due to the effects of recombination between progenitor mitochondrial genomes following generation of a heteroplasmon by parasexual processes [24]. Such mechanisms have been demonstrated to operate in sexual crosses between E. typhina endophytes, although in general, uniparental inheritance is observed [44]. Alternatively, accelerated evolutionary rates of mt DNA relative to nuclear DNA, which have been observed in animals, fungi and in certain protist species [27], may contribute to lower phylogenetic affinities. In support of this explanation, substantial size differences were observed between the mt genomes of the two E. typhina isolates used for this study, suggesting that extensive surveys of intraspecific diversity may be required to identify suitable candidates for progenitor status. However, mt genomes within a given heteroploid taxon were relatively uniform in size, suggesting that limited opportunities for evolutionary divergence have arisen.

In contrast, the N. uncinatum mt genome protein coding gene complement is very closely related to that of one E. typhina lineage, consistent with an origin from this species during heteroploid formation. Furthermore, the mt genome of LpTG-2, which was inferred to have arisen as an N. lolii x E. typhina hybrid [45], demonstrates a high genetic similarity to the N. lolii mt genome, as confirmed by previous studies [37, 46]. This latter observation suggests a relatively recent origin in evolutionary time, and the absence of complicating factors such as recombination between mt genomes.

Phylogeny of novel fescue-derived endophytes

A close relationship was apparent between the previously described taxon FaTG-2 and the novel UNS endophyte group, both of which show affinities to E. festucae and to the putative progenitor of the FGC1/LAE lineages. This result was consistent with the SSR-based phenetic analysis, in which FaTG-2 and UNS accessions were located in sister groups within the same super-cluster of the phenogram [9]. Despite this close affinity, the present study was able to confirm that these two taxa are distinct, based on the formation of separate sub-clusters in both the FGC1 and Group 1 E. festucae-containing clades of the phylograms, as well as in the network diagram. The mt genome phylogram reinforced this distinction, suggesting that the UNS mt genome may have been contributed by the ancestor of the FGC1/LAE lineages, while the FaTG-2 mt genome, as previously described, is most closely related to E. festucae. In combination, the data suggests that these two heteroploid taxa may be derived from hybridisation events in reciprocal mode between two pairs of closely related haploid species, or that divergence from common origins has occurred within each lineage.

Similarly, close relationships between the FaTG-3 and FaTG-3-like endophytes was revealed through the initial endophyte-specific SSR analysis [9]. In the present study, both FaTG-3 and FaTG-3-like endophytes display phylogenetic affinities to both E. typhina and the FGC1/LAE ancestor, with particularly strong similarity for the E. typhina-like gene copies. Similar relationships have been obtained for FaTG-3-like endophytes (later designated as FaTG-4) in a previous study of tub2 phylogeny [14]. Furthermore, both groups also showed E. festucae-like mt genome structure. However, the two groups were identified in differing host grass taxa, FaTG-3 genotypes being detected in Mediterranean hexaploid tall fescue accessions, while FaTG-3-like endophytes were obtained from decaploid tall fescue [9]. At the genomic level, commonly observed deletions of the FGC1/LAE perA gene copy in FaTG-3 were not present in the FaTG-3-like endophytes. In general, the FGC1 copies of each nuclear gene were distinct between FaTG-3 and FaTG-3-like endophytes, suggesting origin of this genomic sub-component from related but distinct taxa.

Variation of nuclear gene structure

Previous phylogenetic studies of tub2, tefA, act1 and several alkaloid biosynthesis genes, including perA were performed by PCR amplification and subsequent sequencing of amplified PCR products [17, 21, 23]. More recently, a study of tub2 phylogeny made partial use of whole genome sequence data [13]. In contrast, the present study was solely based on whole genome sequencing and subsequent independent assembly of entire tub2,tefA, and perA genes, allowing comprehensive identification of insertion-deletion events. Two common deletions were identified in the FGC1-specific perA gene copies from FaTG-2, FaTG-3 and UNS endophytes. A previous study of perA gene phylogeny reported the presence of a 328 bp deletion within the coding region of the FaTG-2 genome [23], similar to observations of the FaTG-2, FaTG-3 and UNS genomes in this study, and an additional adjacent deletion of 627 bp within the same gene copy was here identified for all three taxa. As FaTG-2 endophytes have previously been demonstrated to effectively produce the alkaloid, perA gene function must be due not to the FGC-1 gene copy, but the alternate copy that is putatively derived from an E. festucae-like progenitor.

In addition to these structural changes, the perA sequence of E. baconii endophyte 9707 exhibited an identical deletion to that reported in the E. festucae endophyte E2368 [47], through loss of the reductase domain-encoding sequence at the 3′-terminus. Although this deletion is common among endophytes [3] none of the sequenced novel fescue endophytes were found to contain this deletion.