A novel bioinformatics pipeline to discover genes related to arbuscular mycorrhizal symbiosis based on their evolutionary conservation pattern among higher plants

Most land plants engage in symbiotic associations with fungi (Glomeromycota) that colonize their roots and provide them with phosphate and other mineral nutrients [1]. This association, referred to as arbuscular mycorrhiza (AM), is found in most major taxa of land plants [2], and is thought to have emerged monophyletically in an early progenitor of the vascular plants [3]. The strongest argument for this assumption is the fact that mycorrhizal development requires a conserved signalling pathway that consists of approximately 10 genes that encode receptor components such as SYMRK, and signalling intermediates such as CCaMK [4]. The genes involved in this pathway are conserved between monocots and dicots, and occur also in lycopods and mosses [5],[6], suggesting that the origin of AM dates back to the early vascular plants at the time when land became colonized [3]. This assumption is consistent with the fossil record, which provides evidence for AM-like associations in the sediments of the Rhynie chert that is estimated to originate from the Ordovician period around 450 My ago [7]. More than 350 My after the evolution of AM, a subsequent event in a small subset of the dicots (Fabales, Fagales, Cucurbitales, Rosales), allowed for the emergence of a new form of symbiosis, root nodule symbiosis (RNS) with rhizobacteria [8]-[10]. Interestingly, RNS, as well as the actinorrhizal symbiosis with cyanobacterial endosymbionts [11],[12], involve the same signalling pathway as AM, which therefore is referred to as common symbiosis signalling pathway (common SYM pathway) [4],[13]. A central element of the common SYM pathway is calcium spiking, a rythmic change in perinuclear calcium concentration, which is perceived and transmitted by calcium and calmodulin-dependent protein kinase (CCaMK) to induce symbiotic gene expression [14],[15].

AM is formed by more than 80% of the vascular plants [1], indicating that this association provides a significant selective advantage over non-mycorrhizal plants. However, some plant taxa do not form AM, among them the Brassicaceae with the best-characterized model plant species Arabidopsis thaliana, and the Chenopodiaceae with the economically important crop species Beta vulgaris (sugar beet). AM-related genes are often conserved among AM-competent plant species, while they are less conserved or even missing in non-mycorrhizal species. This phenomenon has been described for VAPYRIN (VPY), which is essential for infection and development of the fungal feeding structures, the arbuscules, in Petunia hybrida as well as in Medicago truncatula [16],[17]. VPY is entirely missing from non-mycorrhizal plant species [16]-[18], and the same is true for numerous genes that are expressed specifically in AM [19],[20]. Such a pattern of conservation was also observed in AM-related genes that are members of large ubiquitous gene families such as the ABC transporters STUNTED ARBUSCULE (STR) and STR2, or the GRAS-type transcription factor REQUIRED FOR ARBUSCULAR MYCORRHIZA1 (RAM1), which both belong to subfamilies which are restricted to AM-competent plant species [21],[22]. Similarly, several components of the common SYM pathway are missing from the non-mycorrhizal model species Arabidopsis thaliana [6], whereas they are conserved among AM-competent dicots and monocots.

Traditionally, AM-related genes have been identified either by mutant screenings followed by characterization of the mutated gene (forward genetics), or by transcript profiling, followed by mutational analysis of AM-induced genes (reverse genetics). Considering the increasing number of sequenced plant genomes, the loss of AM-related genes from the genomes of non-mycorrhizal species could serve as a criterion to detect new AM-related genes by comparative genomics. This third way of gene discovery could potentially identify AM-related genes that have escaped characterization via traditional genetic approaches because of functional redundancy, lethal phenotypes, or constitutive gene expression. Conceptually, such an approach represents a substractive procedure where the proteomes of non-mycorrhizal plants such as A. thaliana are substracted from a panel of proteomes of AM-competent species to result in a set of proteins that are consistently conserved among AM-competent plants and absent from non-mycorrhizal reference species.

Here, we describe a novel approach to identify new AM-related genes based on genome substraction. The approach consists of a multistep procedure that uses protein sequence conservation and gene expression as criteria for the enrichment of potential symbiosis-related genes. We have compared a set of six AM-competent angiosperm species (Solanum lycopersicum, Solanum tuberosum, Vitis vinifera, Medicago truncatula, Glycine max, Populus trichocarpa), and three non-AM species (Arabidopsis thaliana, Arabidopsis lyrata, and Brassica rapa) with an initial clustering approach that resulted in a set of potential candidate proteins comprising approximately 10% of the entire proteome. At a second step, a clustering of these genes based on gene expression patterns in Medicago truncatula provided a set of conserved genes that are induced in AM. Finally, to focus on conserved constitutively expressed genes (such as the common SYM genes), a proteome blast of the conserved genes to an extended panel of proteomes (including monocots) allowed to perform quantitative statistics on the E-values, hence providing a set of proteins that are significantly more conserved among AM-competent plant species than towards the non-mycorrhizal Brassicaceae. This strategy was validated with a number of known AM-related genes that passed our selection scheme. The resulting list of predicted AM-related proteins will be functionally tested by reverse genetic approaches.

During the past decade, many screens for AM- or RNS-defective mutants have been carried out and a significant number of genes involved in signalling and symbiotic functioning have been identified [4]. In addition, transcript profiling has identified many AM- and RNS-related transcripts, however, only few of them have been functionally characterized by reverse genetic approaches. A central finding of the majority of these genetic studies was, that the genes involved in symbiosis are well conserved among AM-competent species, whereas they are conserved to a lesser extent, or entirely missing, in non-mycorrhizal species such as A. thaliana. Based on this finding, we designed a strategy for a third avenue to identify yet unexplored AM-related genes using the conservation of their coding region as the main criterion. We compared the proteomes of AM-competent and non-mycorrhizal plant species, thereby isolating the proteins that are conserved only, or to a higher degree, in AM-competent species. Essentially, this approach represents an in silico substraction procedure, where the proteomes of the non-mycorrhizal species are substracted from the proteomes of the AM-competent species to yield a fraction of the proteome that is expected to be enriched in AM-related proteins.

Since function is selected for at the level of the protein, and since sequence conservation at the DNA level is affected by the degeneracy of the genetic code, we compared the species at the proteome level. Ideally, an in silico substractive approach would involve an equal number of proteomes of the AM-competent and the non-mycorrhizal group. However, due to the wide success of AM in nature, non-mycorrhizal species are the exception (estimated 10-20% of the angiosperm species), and only few of them have been sequenced to a degree that they could be used for such an approach. As non-mycorrhizal species, only A. thaliana, A. lyrata, and B. rapa were available as proteomes in the ENSEMBL database. In contrast, a wide range of AM-competent species is available, among them various dicots (legumes, Solanaceae, etc.) and cereals. Hence, in a first approach, we selected two AM-competent monocots (O. sativa, S. bicolor), four AM-competent dicots (S. lycopersicum, V. vinifera, M. truncatula, and G. max), and the two non-mycorrhizal species A. thaliana and B. rapa.

The first step in this approach was to establish the orthologous relationships between the proteomes of the eight selected species. This task was performed with Hieranoid, which establishes clusters of orthologous/paralogous proteins in a hierarchical step-wise procedure that proceeds according to a previously defined phylogenetic tree of the involved species. The subsequent step is to identify protein clusters that have members in all or most of the AM-competent species, but none in the non-mycorrhizal species. However, after a first run with Hieranoid, we noticed that in many cases, the orthologous trees of known single proteins fell into a monocot and a dicot tree, thus preventing our strategy. Therefore, we decided to perform a clustering just using dicot species, namely G. max, M. truncatula, P. trichocarpa, V. vinifera, S. lycopersicum, and S. tuberosum, together with the three non-mycorrhizal species A. thaliana, A. lyrata, and B. rapa (see Additional file 2: Figure S1). In order to assess the efficiency of the procedure, we used a number of known AM-related genes and a few housekeeping genes for reference (Additional file 5: Table S1).

Clustering of all nine proteomes resulted in a total of 28’528 protein clusters. Subsequent selection of clusters that contain orthologues only in AM-competent species established a first set of potential AM-related proteins (Figure 2, Task6). In order to account for cases were individual orthologues may be missing from one or two proteomes, we relaxed the stringency of the filtering by allowing also for clusters were 1 or 2 members of AM-competent species were missing (Task5 and Task4, respectively) (Table 1, Figure 2). An even more permissive filter (Task3) was applied to obtain a large population of proteins that were subsequently subjected to quantitative assessment of amino acid sequence conservation based on the comparison of E-values (Figure 2, Additional file 12: Table S4, Additional file 13: S5 and Additional file 14: S6). The fact that Task3 and Task4 yielded rather long lists of candidate proteins (>1000 clusters) was not considered a problem, since subsequent filtering and sorting narrowed the list further down (see below).

A next step was to use gene expression as a criterion for selection, hence, we defined expression ratios relevant for AM (Table 2). Based on these criteria we sorted the list of proteins identified by Task4. By applying a stringent threshold induction ratio of 3-fold, we restricted the analysis to genes with a robust transcriptional induction in mycorrhizal roots (Figure 3, Additional file 7: Table S2). Besides known AM-related proteins such as PT4, AM3, VAPYRIN, RAM2, and IPD3, this list included several interesting proteins like the LysM receptor kinase LYR1, a close homologue of the Lotus japonicus nod factor receptor NFR5 [48], and several ABC transporters and peptide transporters, that may play a role in AM development. Highly induced candidates with yet unknown functions in AM include a triacylglycerol lipase (MTR_7g081050.1), neutral ceramidase (MTR_3g079190.1), an enzyme with a role in sphingolipid metabolism [49], and reticuline oxidase (MTR_7g034740.1) an enzyme involved in the biosynthesis of the alkaloid intermediate reticuline [50]. The fact that many known AM-induced genes were recovered by our selection procedure (Additional file 7: Table S2 and S3) documents the potential of gene discovery by comparative phylogenetics.

The main goal of this study was to identify new AM-related proteins based on conservation of the coding sequence, i.e. independently of AM-related gene induction. In order to obtain a quantitative measure for conservation of the open reading frame (ORF) in AM-competent species, we performed a secondary filtering step by blasting the conserved consensus sequences derived from the 6 AM-competent species (MSA) against a wider panel of plant proteomes including dicots, monocots, the liverwort Selaginella moellendorfii, the moss Physcomitrella patens, and the unicellular alga Chlamydomonas reinhardtii. The interest to include these three lower plants is based on the finding that several of the AM-related genes were found in liverworths and/or mosses [5],[6], potentially revealing the evolutionary origin of the common SYM pathway in lower land plants [8]. The psiblast was performed with the results of the relatively permissive Task3 (protein clusters containing at least 3 of the 6 AM-competent species). With this approach 1329 proteins were identified for which the E-values between the three relevant groups of plant species were significantly different (Wilcoxon test, p < 0.05) between the non-mycorrhizal species, and the AM-competent dicots and monocots, respectively (Additional file 12: Table S4, Figure 6), indicating that they are conserved particularly well in AM-competent species. Interestingly, the 50 proteins with the lowest p-values (i.e. highest AM-related conservation) included 13 proteins with annotations as receptors, protein kinases, phosphatases, disease resistant proteins, or G-protein regulating protein, implying roles in signalling.

In view of the essential role of AM in plant nutrition, a highly conserved ammonium transporter (Additional file 12: Table S4, rank 572), and the oligomineral nutrient transporter Nramp5 (rank 495) are of particular interest. Ammonium is the form in which nitrogen is thought to be transferred from the AM fungus to the plant [51],[52], and therefore, ammonium transporters are of particular importance in AM [41]. Transporters for other mineral nutrients such as zinc and copper remain to be identified. Many of the conserved proteins were also found to be well conserved in the liverwort S. moelledorfii, and some of them even in the moss P. patens, whereas C. reinhardtii proved to be phylogenetically much more distant (Additional file 12: Table S4).

In order to obtain a direct quantitative measure for the degree of AM-related conservation, we defined the conservation ratio, which represents the ratio of the averaged E-values from psiblast (as log10) between the non-mycorrhizal species, and the AM-competent species (dicots or monocots, respectively) (Additional file 14: Table S6a). For comparison, the results from Task9 (proteins found in all species) were included as representatives of generally well conserved housekeeping proteins. The larger the conservation ratio, the more a protein sequence is divergent between the non-AM plants and the AM-competent species. Therefore, the proteins identified by Task3 have mostly positive values, whereas the housekeeping proteins (Task9) have neutral or negative values (Figures 6 and 7). Of particular interest were proteins with a conservation ratio >100, which comprised for example PT4 and RAM1 (Figure 7a), which are known to show a highly AM-related conservation pattern (Figure 1). A general overview over the conservation ratios (Additional file 14: Table S6a) revealed that for 220 proteins, the conservation ratio could not be calculated because they were entirely missing from the psiblast results in the non-mycorrhizal species (NaN in Additional file 14: Table S6a). In addition, a further 272 proteins exhibited conservation ratios of >100, indicative of a high degree of AM-related conservation. In order to define a range above which the conservation ratios can be considered as being skewed towards the AM-competent species, 150 proteins identified by Task9 (houskeeping genes) were included for reference in the results from Task3 (see Additional file 14: Table S6a; proteins marked in yellow).

Among the 220 proteins that were missing from the non-mycorrhizal species (Additional file 14: Table S6a), 123 have unknown functions. Interestingly, three proteins each are predicted chitinases (TASK3_1488, TASK3_1548, TASK3_1649), and glucanases (TASK3_566, TASK3_581, TASK3_4113), indicating potential roles in cell wall modification or signaling. Three predicted CLE peptides (TASK3_3616, TASK3_3634, TASK3_4020) might be involved in AM-related signaling as it was shown for RNS [53]. A known component in symbiotic signaling comprised INTERACTING PROTEIN of DMI3 (IPD3)/ CYCLOPS (TASK3_1003), an interactor of the central symbiosis regulator CCaMK encoded by DMI3 [54]-[56]. High conservation ratios were also found for VAPYRIN (TASK3_795), which is required both for AM and RNS [16],[17],[26] and VAPYRIN-like (TASK3_1469).

16 proteins exhibited a value >100 for both conservation ratios (Additional file 14: Table S6b), i.e. they are highly conserved among both, AM-competent dicots and monocots, relative to the non-mycorrhizal species. Of particular interest is TASK3_196 which encodes a predicted α-glucosidase (Additional file 12: Table S4 and Additional file 14: S6b) which is in some species classified as α-xylosidase (e.g. XP_006358190.1 in S. tuberosum). While the conserved protein is missing from A. thaliana (Figure 8a), there are related glycosidases that occur in all examined plant species including A. thaliana, indicating that the non-mycorrhizal model species misses only the AM-related isoform (Figure 8b). The potential role for α-glucosidase or α-xylosidase in AM is unknown, but given the important contribution of xylose in the hemicelluloses of the cell wall [57], it may be involved in the modifications of cell walls during AM infection and functioning [20]. Interestingly, a predicted α-xylosidase of petunia is among the most strongly repressed genes under conditions of high phosphate, which inhibits AM [19].

Five of the 16 proteins with high conservation ratios (>100) represent receptors and protein kinases, of which three (Task3_354, Task3_238, and Task3_4055) are closely related receptor-like protein kinases (RLKs) (Additional file 14: Table S6b). They are similar to S-locus receptor kinase (SRK) [58], which mediates sporophytic self-incompatibility in the Brassicaceae [59], and to the maize receptor kinase ZmPK1 [60]. Phylogenetic analysis showed that they indeed fall into a large clade of proteins found in all AM-competent species, but with considerably lower conservation in the non-mycorrhizal species (Additional file 11: Figure S5). A conserved C-rich stretch in the amino acid sequence is shared among the three AM-related RLKs, as well as with SRK, and ZmPK1 (Additional file 18: Figure S6), indicating that they may have similar tertiary structures mediated by cysteine bridges. Interestingly, ZmPK1 is expressed in roots and young seedlings, rather then floral tissues [60] suggesting that this family of receptor kinases can have other functions than self-incompatibility. Taken together, the examples of α-glucosidae/xylosidase and the RLKs show that our approach can identify proteins with AM-related conservation patterns even if they are members of large families. Furthermore, the case of the three related RLKs highlights the potential of our approach to identify proteins that could potentially escape genetic mutant screens because of functional redundancy.

Transcription is controlled primarily by cis-acting regulatory sequences in the promoters of the genes. Hence, co-regulated genes can be expected to share common regulatory elements in their promoters. Based on this logic, we screened the promoters of 190 AM-induced genes for overrepresented sequences that may represent cis-acting regulatory promoter elements. The prediction program MEME identified two new palindromic elements (elements Myc1 and Myc2) and found the GCCGGC element to be overrepresented in AM-induced promoters as well as in the promoters of genes identified by Task4. Statistical analysis confirmed the overrepresentation of these elements in AM-related promoters (Additional file 16: Table S7). The statistical significance of element Myc1 and Myc2 is emphasized by their rather unusual length (14 nt and 13 nt, respectively) compared to other cis-acting elements that are typically less then 10 nt in length [61]. In addition, element Myc1 and Myc2 are often coupled (see Additional file 17: File S5), underscoring their relevance for AM-related gene induction, and suggesting that they may function in conjunction.

The predominant position of the GCCGGC element in a range of 500 bp upstream of the ATG start codon is indicative of functional relevance as well. The GCCGGC element can be bound by the transcription factor RAP.2.11 in ethylene-responsive promoters involved in the potassium starvation response of A.thaliana [62], and it was classified as wound- and pathogen-inducible [63]. GCCGGC represents a palindromic variant of the common GCC element (sequence GCCGCC) that is widely found in the promoters of pathogenesis-related genes and in ethylene-responsive promoters of various plants [64]-[70]. In fact, GCC elements occur in different variants, which are recognized by specific members of the ETHYLENE RESPONSE FACTOR family of transcription factors [69]. The overrepresentation of the GCCGGC element in AM-related promoters may explain the frequency of predicted stress- and pathogen-related transcripts in the transcriptome of mycorrhizal roots in plant species such as rice, potato, petunia, parsley, M. truncatula, and L. japonicus [19],[32]-[35],[71]-[73].

Surprisingly, the AM-related CTTC/MYCS element (sequence: CTTGTTC) [46],[74] was not found by MEME in the 190 AM-induced promoters. Furthermore, it was not overrepresented in the genes resulting from Task4, and only moderately in the AM-induced genes (Additional file 16: Table S7). However, an extended form (GACTTGTTC), was found to be significantly overrepresented, indicating that this may represent the complete AM-inducible element in M. truncatula. An overview over the occurrence of predicted promoter elements in the analyzed promoters is provided in Additional file 19: Table S8.

Taken together, we have identified two new potentially AM-associated cis-regulatory promoter elements. In addition, we identified the GCCGGC element as a potential AM-related element that may indicate an involvement of ethylene in AM development, and we present a modified (elongated) form of the CTTC element. Future work will address the functional relevance of these elements for AM development and functioning,

We describe a non-targeted approach to identify new AM-related genes by genome-wide comparative sequence analysis of AM-competent and non-mycorrhizal plant species. The validity of this approach has been confirmed by the finding that it identified many known AM-related genes from the genomes of AM-competent plant species. Following this strategy, we identified many new potentially AM-related genes based on the conservation pattern of the coding sequence, on gene expression pattern, and on predicted regulatory promoter elements. These genes will be further evaluated for their predicted function (pathway), the degree of their conservation, and the copy number, or the number of similar genes in the case of gene families. This process will result in a short list of genes that will be further functionally analyzed by insertional mutagenesis in P. hybrida, M. truncatula or L. japonicus using established transposon-related protocols [74]-[76].

Microarray data are available in the Medicago truncatula Gene Expression Atlas (MtGEA version 3; http://mtgea.noble.org/v3). Plant proteomes were downloaded from the Ensembl Plants website (www.ensemblgenomes.org; database release 20). Promoter sequences of M. truncatula were downloaded from Ensembl Plants BioMart (http://plants.ensembl.org/biomart/martview). Phylogenetic trees were prepared with the online tool www.phylogeny.fr using the standard mode (“One-Click”) with the sequences provided in Additional file 1: File S1. All Supporting Information has been deposited on LabArchives (http://www.labarchives.com).