Genetics of nodulation in Aeschynomene evenia uncovers mechanisms of the rhizobium–legume symbiosis

Among legumes (Fabaceae) capable of nitrogen-fixing nodulation, several Aeschynomene spp. use a unique symbiotic process that is independent of Nod factors and infection threads. They are also distinctive in developing root and stem nodules with photosynthetic bradyrhizobia. Despite the significance of these symbiotic features, their understanding remains limited. To overcome such limitations, we conduct genetic studies of nodulation in Aeschynomene evenia, supported by the development of a genome sequence for A. evenia and transcriptomic resources for 10 additional Aeschynomene spp. Comparative analysis of symbiotic genes substantiates singular mechanisms in the early and late nodulation steps. A forward genetic screen also shows that AeCRK, coding a receptor-like kinase, and the symbiotic signaling genes AePOLLUX, AeCCamK, AeCYCLOPS, AeNSP2, and AeNIN are required to trigger both root and stem nodulation. This work demonstrates the utility of the A. evenia model and provides a cornerstone to unravel mechanisms underlying the rhizobium–legume symbiosis.

L egumes (Fabaceae) account for~27% of the world's primary crop production and are an important protein source for human and animal diets. This agronomic success of legumes relies on the capacity of many species to establish a nitrogen-fixing symbiosis with soil bacteria, collectively known as rhizobia, forming root nodules 1 . Promoting cultivation of legumes and engineering nitrogen fixation in other crops will decrease the input of chemical nitrogen fertilizers and to will help to achieve short-and long-term goals aimed at a more sustainable agriculture 2 .
Intensive research mainly performed on two temperate model legumes, Medicago truncatula and Lotus japonicus, has yielded significant information on the mechanisms controlling the establishment and functioning of the legume-rhizobium symbiosis 1 . In the general scheme, plant recognition of key rhizobial signal molecules, referred to as Nod factors, triggers a symbiotic signaling pathway leading to the development of an infection thread that guides bacteria inside the root and to the distant formation of a nodule meristem where bacteria are delivered and accommodated to fix nitrogen 1 .
To further advance our understanding of the rhizobial symbiosis, there is a great interest in tracking the origin of nodulation 3,4 and in uncovering the whole range of symbiotic mechanisms 5,6 . In this quest, some semi-aquatic tropical Aeschynomene species constitute a unique symbiotic system because of their ability to be nodulated by photosynthetic bradyrhizobia that lack the canonical nodABC genes, necessary for Nod factor synthesis 7,8 . In this case, nodulation is not triggered by a hijacking Type-3 secretion system present in some non-photosynthetic bradyrhizobia 9,10 . Therefore, the interaction between photosynthetic bradyrhizobia and Aeschynomene represents a distinct symbiotic process in which nitrogen-fixing nodules are formed without the need of Nod factors. To unravel the molecular mechanisms behind the so-called Nod factorindependent symbiosis, Aeschynomene evenia (400 Mb, 2n = 2x = 20) has emerged as a genetic model [11][12][13] .
A. evenia is also a valuable legume species because: (i) it uses an alternative infection process mediated by intercellular penetration as is the case in 25% of legume species 14,15 ; (ii) it is endowed with stem nodulation, a property shared with very few hydrophytic legume species 16,17 ; and (iii) it groups with Arachis spp., including cultivated peanut (Arachis hypogaea) in the Dalbergioid clade, which is distantly related to L. japonicus and M. truncatula 11 . Previous transcriptomic analysis from root and nodule tissues did not detect expression of several known genes involved in bacterial recognition (e.g., LYK3 and EPR3), infection (e.g., RPG and FLOT), and nodule functioning (e.g., SUNERGOS1 and VAG1) 12,18 . Such data support the presence of distinct or divergent symbiotic mechanisms in A. evenia in comparison with other well-studied model legumes. In addition, they comfort A. evenia as a system of interest to study the evolution and diversity of the rhizobial symbiosis.
In this work, to efficiently conduct genetic studies of nodulation in A. evenia, we produce a genome sequence for this species along with de novo RNA-seq assemblies for 10 additional Nod factor-independent Aeschynomene spp. These genomic and transcriptomic datasets allow us to perform a comparative analysis of known symbiotic genes, leading to the evidence of singular symbiotic mechanisms in Aeschynomene spp. Finally, we use the available genome sequence in a forward genetic approach to conduct the genetic dissection of nodulation in A. evenia and we identify a receptor-like kinase that is not present in model legumes. This finding uncovers an important molecular step in the establishment of the Nod fcator-independent symbiosis.

Results
A reference genome for the Nod factor-independent legume Aeschynomene evenia. As a support to forward genetic and comparative genetic studies of nodulation, a reference genome assembly was produced for A. evenia using the inbred CIAT22838 line 12 . To the single-molecule real-time (SMRT) sequencing technology from PacBio RSII platform was used to obtain a 78× genome coverage (Supplementary Tables 1-4). The resulting assembly was 376 Mb, representing 94-100% of the A. evenia genome, considering the estimated size of 400 Mb obtained by flow cytometry 12,16 or of 372 Mb derived from k-mer frequencies ( Supplementary Fig. 1). PacBio scaffolds were integrated in the 10 linkage groups of A. evenia using an existing genetic map 12 , an ultra-dense genetic map generated by genotyping-by-sequencing (GBS), and scaffold mapping was subsequently refined on the basis of synteny with Arachis spp. 19 (Supplementary Figs. 2 and 3). The final 10 chromosomal pseudomolecules anchored 302 Mb (80%) of the genome (Supplementary Fig. 4 and Supplementary Table 4). Protein-coding genes were annotated using a combination of ab initio prediction and transcript evidence gathered from RNA sequenced from nine tissues/developmental stages of nodulation using both RNA-sequencing (RNA-seq) and PacBio isoform sequencing (Iso-Seq) (Supplementary Tables 5 and 6). The current annotation contains 32,667 gene models ( Fig. 1a and Supplementary Table 7). Their expression pattern was also determined by developing a Gene Atlas from the RNA-seq data obtained here (Supplementary Table 8) and from an earlier nodulation kinetics 18 . The identification of 94.4% of the 1440 genes in the Plantae BUSCO dataset (Supplementary Table 9) confirmed the high quality of the genome assembly and annotation. Approximately 72% of the genes were assigned functional annotations using Swissprot, InterPro, Gene Ontology (GO), and KEGG (Supplementary Table 10). Additional annotation of the genome included the prediction of 6558 non-coding RNAs (ncRNAs), the identification of repetitive elements accounting for 53.5% of the assembled genome and mainly represented by LTRs, the effective capture of 16 out of the 20 telomeric arrays, and the distribution of sequence variation along chromosomes based on the resequencing of 12 additional A. evenia accessions 20 (Fig. 1a, Supplementary Fig. 4, and Supplementary Tables 11-15). Finally, all the resources were incorporated in the AeschynomeneBase (http://aeschynomenebase.fr), which includes a genome browser and user-friendly tools for molecular analyses.
To trace back the history of the A. evenia genome, it was compared to the genomes of Arachis duranensis and Arachis ipaiensis, which belong to the same Dalbergioid clade. Aeschynomene and Arachis lineages diverged~49 Ma (million years ago) but are assumed to share an ancient whole-genome duplication (WGD) event that occurred~58 Ma at the basis of the Papilionoid legume subfamily 19,[21][22][23] . The shared WGD event, the Aeschynomene-Arachis divergence, and the A. duranensis-A. ipaiensis speciation were apparent in the synonymous substitutions in coding sequence (Ks) analysis between and within the A. evenia-A. duranensis-A. ipaiensis genomes (Fig. 1b). Modal Ks values are~0.65 for A. evenia, i.e., more similar to those reported for L. japonicus and G. max (both~0.65) than to those of A. duranensis (~0.85) and A. ipaiensis (~0.80) that were already reported to have evolved relatively rapidly 19 . In the case of A. evenia, it is worth noting that no more recent peak of Ks is visible, indicating it did not undergo any further WGD event. We identified paralogous A. evenia genes and orthologous A. evenia-Arachis spp. genes using synteny and Ks value criteria. This revealed the blocks of conserved collinear genes resulting from the WGD event~58 Ma in the A. evenia genome (Fig. 1c). A comparison of A. evenia with A. duranensis and A. ipaiensis shows that extensive synteny remains prominent along chromosome arms despite multiple rearrangements (Fig. 1d). To be able to compare A. evenia to other Aeschynomene spp. that also use a Nod factor-independent process, we performed de novo RNA-seq assemblies from root and nodule tissues for 10 additional diploid Aeschynomene spp. (Supplementary Tables 16 and 17). Groups of orthologous genes for A. evenia, related Aeschynomene spp., and several species belonging to different legume clades were then generated using OrthoFinder (Supplementary Table 18). A consensus species tree inferred from single-copy orthogroups perfectly reflected the legume phylogeny and, for the Aeschynomene clade, the previously observed speciation with the early diverging species Aeschynomene filosa, Aeschynomene tambacoundensis, and Aeschynomene deamii, and a large group containing A. evenia 16,20 (Fig. 2a).
Symbiotic perception, signaling, and infection. In addition to their ability to nodulate in the absence of Nod factors 8,11 , A. evenia and related Aeschynomene spp. use an infection process that is not mediated by the formation of infection threads 14 . This prompted us to perform a phylogenetic analysis of known symbiotic genes 1 based on the orthogroups containing Aeschynomene spp. and to exploit the Gene Atlas developed for A. evenia. This comparative investigation revealed that the two genes encoding the Nod factor receptors, NFP and LYK3, are present but that LYK3 is barely expressed in A. evenia (Supplementary Data 1). What is more, transcripts of both genes were not detected in the transcriptome of all other Aeschynomene species (Fig. 2a). In line with this observation, the gene coding for NFH1 (Nod factor hydrolase 1), which mediates Nod factor degradation in M. truncatula, was not found in any Aeschynomene data ( Fig. 2a and Supplementary Data 2). Interestingly, EPR3, which inhibits infection of rhizobia with incompatible exopolysaccharides in L. japonicus, was not found in the A. evenia genome (Fig. 2a)   Considering that AeLYK3 is probably not involved in Nod factorindependent symbiosis, it remains to be investigated how the presence of two copies of AeSYMRK and AePUB1 in A. evenia might contribute to the diversification of the signaling mechanisms 25 . Downstream of SYMRK, the symbiotic signaling pathway leads to the triggering of the plant-mediated rhizobial infection. Determinants such as VPY, LIN, and EXO70H4 26 , which are required both for polar growth of infection threads and subsequent intracellular accommodation of symbionts in M. truncatula 1 , have symbiotic expression in A. evenia (Supplementary Data 2). This expression pattern is probably linked to the later symbiotic process since rhizobial invasion occurs in an intercellular manner in A. evenia 14 . In contrast, other key infection genes 1 are expressed at very low levels, as is the case of NPL and CBS1, or absent: RPG was undetectable in Dalbergioid legume species and FLOT genes were completely missing in Aeschynomene spp., suggesting mechanistic differences in the infection process (Fig. 2a, Supplementary Fig. 12, and Supplementary Data 2).
Nodule development and bacterial accommodation. During nodule development, differentiating plant cells undergo endoreplication leading to an increase in ploidy levels and cell size. The mitotic inhibitor CCS52A, a key mediator of this nodule development process 27,28 , is conserved in all Aeschynomene spp. (Fig. 3a). However, earlier transcriptomic studies 12,18 failed to detect two genes coding for components of the DNA topoisomerase VI complex, subunit A (SUNERGOS1) and an interactor (VAG1). In L. japonicus, these two genes are required for cell endoreplication during nodule formation 29,30 . From previous Arabidopsis studies, the DNA topoisomerase VI is known to contain two other components, the subunit B (BIN3) 31 and a second interactor (BIN4) 32 , which were both successfully identified in legumes but not in A. evenia (Fig. 3a). Synteny analysis based on genomic sequence comparison with Arachis spp. substantiated the specific and complete loss of SUNERGOS1, BIN3, and BIN4, and the partial deletion of VAG1 in A. evenia (Supplementary Fig. 13 and Supplementary Data 2). A similar pattern could be observed for most Aeschynomene spp. However, SUNERGOS1, BIN3, BIN4, and VAG1 with a distinct truncation were detected in A. deamii and the full gene set was present in A. filosa and A. tambacoundensis as is the case for peanut (A. hypogaea), indicating that these gene losses are disconnected from the Nod factor-independent character ( Fig. 3a and Supplementary Fig. 14). To link these different gene patterns with variations in nodule cell endoreplication, roots and nodules of several species were analyzed by flow cytometry. Contrary to our expectations, whereas no difference in ploidy levels was observed in A. filosa, A. tambacoundensis or peanut, we discovered higher ploidy levels in nodule cells than in root cells of A. deamii, A. evenia, A. scabra, A. selloi, and A. sensitiva (Fig. 3b). Taken together, these data reveal a case of gene co-elimination affecting the Topoisomerase VI complex, but the functional relevance of this loss of genes on the nodule cell endoreplication process needs to be investigated. Nodule formation is also accompanied by the differentiation of nodule cell-endocyted rhizobia into nitrogen-fixing bacteroids. In M. truncatula, this differentiation is mediated by the expression of a wide set of plant genes coding for nodule-specific cysteinerich peptides (NCRs) 1 . Although NCRs were long thought to be restricted to the IRLC clade to which M. truncatula belongs 33 , A. evenia and other Aeschynomene spp. were recently shown to express NCR-like genes 34 . We identified 58 such genes in the A. evenia genome ( Fig. 4a and Supplementary Data 3). The AeNCR genes are mainly organized in clusters (Fig. 4b) and they are typically composed of two exons encoding the signal peptide and the mature NCR (Fig. 4c). Most NCR genes display prominent nodule-induced expression in A. evenia that correlates with the onset of bacteroid differentiation ( Fig. 4d and Supplementary Data 3). All predicted NCRs contain one of the two previously described cysteine-rich motifs 34,35 (Fig. 4e). Thus, 26 NCRs of A. evenia harbor the cysteine-rich motif 1 similar to M. truncatula NCRs while 32 NCRs of A. evenia have the defensin-like motif 2 ( Fig. 4a, Supplementary Figs. 15 and 16). In A. duranensis and A. ipaiensis, no NCRs with the cysteine-rich motif 1 could be found, whereas 10 and 5 NCR-like genes, respectively, with the defensinlike motif 2 were identified (Fig. 4a, Supplementary Figs. 15 and 16) and the expression of most of them is induced in the nodule (LegumeMines database). These features of Dalbergioid NCRs raise questions as to how they emerged and whether they evolved for symbiotic or defense functions.
Nodule functioning involves leghemoglobins derived from class 1 phytoglobins. In nitrogen-fixing nodules, maintaining a low but stable oxygen concentration is crucial to protect the nitrogenase complex. To ensure this function, legumes have recruited leghemoglobins (Lbs) that evolutionary derive from non-symbiotic hemoglobins (now termed phytoglobins; Glbs), and that occur at high concentrations in nodules 36 . We found six globin genes in the A. evenia genome (Supplementary Data 4). Two of them are homologous to class 3 Glb genes and were not studied further. Two genes show moderate expression, have homology to class 1 and class 2 Glb genes, and were accordingly designated AeGlb1 and AeGlb2 (Fig. 5a). The two other globin genes were found to be highly and almost exclusively expressed in nodules (Fig. 5a). This observation suggested that they encode Lbs and were termed AeLb1 and AeLb2. To unequivocally classify the four proteins, they were purified and characterized for heme coordination (Fig. 5b). Both AeGlb1 and AeGlb2 showed hexacoordination in the ferric and ferrous forms, confirming that they correspond to class 1 and class 2 Glbs, respectively. AeLb1 shows complete pentacoordination in both ferric and ferrous form, whereas AeLb2 is hexacoordinate in the ferric form and almost fully pentacoordinate in the ferrous form. AeLb1 is therefore a typical Lb but AeLb2 appears to be an unusual one. All four globins were found to bind the physiologically relevant ligands, O 2 and nitric oxide ( Supplementary Fig. 17). However, the unexpected discovery was that both AeLb1 and AeLb2 cluster with class 1 Glbs and not with class 2 globins, as observed for other legumes 36 (Fig. 5c, d). In the globin phylogeny, AeLb1 and AeLb2 cluster tightly with certain class 1 Glb genes of Arachis and also of the more distantly related legume Chamaecrista fasciculata. The Arachis genes are also highly expressed in nodules (LegumeMines database) and probably encode Lbs. Among the C. fasciculata genes, one was previously evidenced to be highly expressed in root nodules and to code for a putative ancestral Lb named ppHB 37 (corresponds to the Chafa1921S17684 gene in Fig. 5c). Sequence and synteny analysis further indicated that A. evenia Lbs and class 1 Glb genes are similar and located in a single locus that is conserved in legumes . This supports the hypothesis that A.evenia Lbs arose from class 1 Glbs by local gene duplication, and the presence of probably such Lbs in Arachis and Chamaecrista legumes further suggests this evolution to be ancient. The finding of Lbs originating from a class 1 Glb challenges our view on the evolution of Lbs in legumes and is only comparable to panHBL1, the symbiotic Glb1 of the non-legume Parasponia 3 . However, panHBL1 appears to be different from A. evenia Lbs (Fig. 5c). These Lbs thus offer a valuable case to study the convergent evolution of O 2transporting Lbs.
Genetic dissection of root and stem nodulation. To uncover genes underpinning the singular symbiotic traits evidenced in A.

Dataset
Signaling The collection of mutants was subjected to targeted sequence capture on a set of selected genes with a potential symbiotic role. Analysis of EMS-induced SNPs allowed the filtering of siblings originating from the same screening bulks and led to the identification of candidate mutations. We decided to focus our genetic work on the Nod − mutants since they are most probably altered in genes controlling the early steps of nodulation. Moreover, they provide an opportunity to test the role of these genes in stem nodulation whose genetic control is completely unknown so far. For this, Nod − mutants were backcrossed to the WT line and segregating F 2 progenies were phenotyped for root and stem nodulation after sequential inoculation. These analyses always pointed to a single recessive gene controlling both root and caulinar nodulation ( Fig. 6a and Supplementary Table 20). For Nod − mutants associated with a candidate mutation, the mutations were validated as being causative by genotyping F 2 backcrossed mutant plants and by performing targeted allelism tests. This produced allelic mutant series for five genes of the symbiotic signaling pathway 1 : AePOLLUX (6 alleles), AeCCaMK (4 alleles), AeCYCLOPS (2 alleles), AeNSP2 (4 alleles), and AeNIN (6 alleles) (Fig. 6b and Supplementary Tables 20-22). Among these signaling genes, AePOLLUX was found to be consistently expressed in all plant organs, whereas the other genes are expressed only in symbiotic organs. AeCCaMK is constantly expressed in roots and in all stages of nodule development, AeCYCLOPS and AeNIN are induced during nodulation, and AeNSP2 is down-regulated during nodulation (Fig. 6c and Supplementary Data 2). Thus, mutant analysis revealed that the signaling pathway, described in M. truncatula and L. japonicus, is partially conserved in A. evenia and is necessary for stem nodulation. However, not all known signaling genes were evidenced with the mutant approach (Fig. 6d). In particular, no consistent mutation was found in any member of the LysM-RLK family. Although it cannot be excluded that our mutagenesis was not saturating, this observation again supports the lack of a key role for LysM-RLKs in the Ae05 Ae06 Ae07 Ae08 Ae09 Ae10

Signal pepƟde Mature NCR
Aeschynomene NCR motif 1 SP-X n -C-X n -C-X 3 -C-X n -C-X n -C-X 1 -C-X n Medicago NCR motif a SP-X n -C-X 5 -C-X n -C-X n -C-X 4 -C-X 1 -C-X n Aeschynomene NCR motif 2 SP-X n -C-X n -C-X n -C-X 3 -C-X n -C-X n -C-X 1 -C-X 3 -C-X n Defensin-like motif SP-X n -C-X n -C-X n -C-X 3 -C-X n -C-X n -C-X 1 -C-X 3 -C-X n e Fig early steps of the symbiotic interaction in A. evenia. Neither was a causative mutation found for the two paralogs of SYMRK in A. evenia. In an earlier study, we used RNAi to target AeSYMRK (actually AeSYMRK2), which reduced the number of nodules 13 . Because AeSYMRK1 and AeSYMRK2 are 82% identical in the 296-pb RNAi target region, they were probably both targeted. The functioning of the two receptors during nodulation remains to be investigated.
A receptor-like kinase mediates the symbiotic interaction. Two Nod − mutants, defective in both root and stem nodulation, were not associated with any known genes and were consequently good candidates to uncover novel symbiotic functions (Fig. 7a). To identify the underlying symbiotic gene, we used a mapping-bysequencing approach on bulks of F 2 mutant backcrossed plants.
Linkage mapping for each mutant population identified the same locus on chromosome Ae05, where mutant allele frequencies reached 100% (Fig. 7b). Analysis of the region containing the symbiotic locus identified mutations in a gene that encodes a cysteine-rich receptor-like kinase (CRK) 38 , henceforth named AeCRK (Supplementary Table 23). The predicted 658-aa-long protein harbors a signal peptide, two extracellular DUF26 (domains of unknown function) domains, a transmembrane domain (TM), and an intracellular serine/threonine kinase domain ( Fig. 7c and Supplementary Fig. 22). In the mutated forms, the G2228A SNP alters a canonical intron/exon splice boundary probably generating a truncated protein while the G1062A SNP leads to the replacement of G354E in the highly conserved glycine-rich loop of the kinase domain (Supplementary  Table 20). Allelism tests performed with the two Nod − mutant lines (I10 and J42) indicated that they belong to the same complementation group (Supplementary Table 22). Hairy root transformation of the I10 mutant with the coding sequence of AeCRK, fused to its native promoter, resulted in the development of nodules upon inoculation with the Bradyrhizobium ORS278 strain, while no nodules were produced in control plants transformed with the empty vector ( Supplementary Fig. 23 and Supplementary Table 24). The identification of genetic lesions in the two independent Aecrk alleles together with the transgenic complementation of the mutant phenotype provide unequivocal evidence that AeCRK is required for the establishment of symbiosis A. evenia.
AeCRK was found to be expressed in roots with significant upregulation in nodules, in agreement with its symbiotic function ( Fig. 7d and Supplementary Data 5). Notably, AeCRK is part of a cluster of five CRK genes in A. evenia, but genes of this cluster are interspersed within the CRK phylogeny ( Fig. 7e and Supplementary Data 5). Although similar CRK clusters are located in syntenic regions in other legumes, no putative ortholog to AeCRK could be found in M. truncatula or L. japonicus, and actually in no Papilionoid legume using a root hair-and infection threadmediated infection process (Fig. 7f, Supplementary Fig. 24, and Supplementary Data 5). To gain further insights into the molecular evolution of AeCRK, we ran branch model by estimating different substitution rates (ω) using the phylogenetic tree topology. These analyses, performed on the entire gene sequence and on the four functional domains of AeCRK orthologs separately (signal peptide, extracellular, transmembrane, and kinase domains), revealed a higher purifying/negative selection acting on the extracellular domain part in the Aeschynomene clade (ω BG = 0.480 and ω FG = 0.187, p = 0.017214) (Fig. 7f and Supplementary Table 25). This purifying selection suggests that AeCRK could have evolved to adapt nodulation with Nod genelacking photosynthetic bradyrhizobia. These data support that AeCRK is a key component of the pathway used by A. evenia to trigger symbiosis in the absence of Nod factors and infection threads.

Discussion
A. evenia and a handful of other Aeschynomene spp. have gained renown for triggering efficient nodulation without recognition of rhizobial Nod factors nor infection thread formation 8,9,14 . To accelerate the deciphering of this original symbiosis, we conducted in A. evenia forward genetics based on an EMS mutagenesis and developed a reference genome sequence to enable resequencing strategies of nodulation mutants. This work leads to the demonstration that the triggering of nodulation in A. evenia is mediated by several components of the Nod signaling pathway described in model legumes, AePOLLUX, AeCCaMK, AeCY-CLOPS, AeNSP2, and AeNIN, thus significantly extending a previous report of the involvement of AeSYMRK, AeCCaMK, and AeLHK1 genes in root nodulation 13 . The present study also reveals that this symbiotic signaling pathway controls not only root but also stem nodulation in A. evenia. This dual nodulation is present in few half-aquatic legume species 16 nodulation has remained unknown so far. With the forward genetic screen, not all known genes of the Nod signaling pathway were recovered. Indeed, no causative mutation could be found in AeCASTOR or in AeNSP1, whereas CASTOR and NSP1 are known to act in concert with POLLUX and NSP2, respectively 1 .
In addition, there are no obvious paralogs reported that may function redundantly, as it is probably the case for SYMRK in A. evenia. Therefore, either both these genes were unfortunately not targeted by the EMS mutagenesis or a special evolution of AePOLLUX and AeNSP2 rendered them sufficient for symbiosis as evidenced for DMI1/POLLUX in M. truncatula 39 . Also striking is the failure of the mutant approach to demonstrate the involvement of any LysM-RLK member, most notably the Nod factor receptors. In agreement with this observation, LYK3 is not expressed in A. evenia. Conversely, NFP remains expressed in A. evenia, putatively because of a function in the arbuscular mycorrhizal (AM) symbiosis, which is likely ancestral 40 . Therefore, a comparative genetic analysis of NFP and LYK3 between A. evenia and Aeschynomene patula, which displays a Nod factordependent nodulation and which was recently selected as a suitable complementary model 16 , should illuminate their recent evolution and clarify if NFP has any role in A. evenia.
Finding that the core Nod signaling pathway, but not the upstream Nod factor receptors, is conserved in A. evenia suggests that one main difference with other legumes comes from the symbiotic receptor plugged-in in the pathway 41 . In line with this idea, a receptor-like-kinase belonging to the large CRK family 38 was discovered as being required to trigger nodulation in A. evenia. In the legume phylogeny, this gene is present only in Papilionoid lineages using an intercellular infection process and also in the Caesalpinioid legumes, C. fasciculata and Mimosa pudica. Such distribution of the AeCRK orthologs suggests that their presence is ancestral in legumes. Molecular evolutionary analysis further evidenced the extracellular domain of AeCRK orthologs to be under purifying selection in the Aeschynomene clade, arguing for a particular evolution with the Nod factorindependent symbiosis. CRKs are repeatedly pointed as important actors of plant early signaling during immunity and abiotic stress 42,43 . They are supposed to be mediators of reactive oxygen species (ROS)/redox sensing through their DUF26 extracellular domains and to transduce the signal intracellularly via their cytoplasmic kinase 43 . Another putative function of their DUF26 domains was recently proposed, based on strong similarity to fungal lectins, as mediating carbohydrate recognition 44 . Therefore, characterization of AeCRK will be crucial to provide information on pending questions: Has AeCRK retained the ancestral function or has it been neofunctionalized? Is AeCRK involved in the direct perception of photosynthetic bradyrhizobia or does it mediate ROS/redox sensing during early signaling/infection? Is the Nod factor-independent activation inherently linked to intercellular infection? This could be probably the case since genetic studies in L. japonicus evidenced that double-mutant lines were occasionally able to develop nitrogen-fixing nodules in a Nod factor-and infection thread-independent fashion 45 . Additionally, the ability of L. japonicus to be infected intracellularly or intercellularly, depending on the rhizobial partner, was recently used to provide insights into the genetic requirements of intercellular infection 46 . It was showed that some determinants required for the infection thread-mediated infection are dispensable for intercellular infection, among which RPG is found. This finding echoes the observed absence of RPG in A. evenia and other Dalbergioid legumes for which intercellular infection is the rule. However, other infection determinants (LIN, VPY, EXO70H4, and SYN) that are also involved in intracellular accommodation of symbionts are present in A. evenia, suggesting that both the core symbiotic signaling pathway and the machinery mediating intracellular accommodation are conserved, as a general feature of endosymbioses 47 . Continuing the mutantbased gene identification in A. evenia will increase our knowledge on the mechanisms of the as yet under-explored intercellular infection process.
In addition to the intercellular infection process, several symbiotic features present in A. evenia are shared with other legumes, including peanut, for which the molecular basis of nodulation is subject of recent investigations 48 . As evidenced previously 34 and in the present work, Aeschynomene and Arachis spp. express NCR-like genes during bacterial accommodation, in a similar fashion to IRLC legumes, but their symbiotic involvement remains to be clarified. Most remarkable is the discovery that Aeschynomene and Arachis spp. have recruited some class 1 Glbs as Lbs transporting O 2 in nodule infected cells. Indeed, it is well established that in legumes some class 2 Glbs have evolved to Lbs to ensure such a crucial function 36 , but the Dalbergioid lineage appears to be an exception to this pattern of Lb utilization. Comparative genomic analysis in Papilionoid legumes revealed a striking parallel with the presence of two conserved loci where both Glb and Lb genes belonging to class 1 and class 2, respectively, can be found across species. It is therefore tempting to hypothetize that Lbs arose from Glbs by gene tandem duplication and divergent evolution in these two loci, and that they were differentially lost depending on the legume lineages. In Caesalpinioid C. fasciculata, the presence of a hemoglobin that has some characteristics of Lb 37 and is closely related to Dalbergioid Lbs supports that this feature is ancient in legumes. In addition, the presence also in nodulating non-legume species of class 1derived Lbs (e.g., Parasponia) or class 2-derived Lbs (e.g., Casuarina) suggests this dual evolution to be recurrent 3 . This will be an exciting evolutionary issue to determine how different Glbs adapted to Lbs, and if these Lbs have any specific functional specificity.
The discovery of alternative mechanisms underpinning the nitrogen-fixing symbiosis strengthens A. evenia as a valuable model for the study of nodulation. The successful development of a forward genetic approach supported by a reference genome and companion resources also shows this legume is amenable for genetic research, this research being complementary to the one performed on M. truncatula and L. japonicus. The acquired knowledge will contribute to characterize the diversity of the symbiotic features occurring in legumes. It is also expected to benefit legume nodulation for agronomic improvement and, ultimately, it could provide leads to engineer nitrogen-fixation in non-legume crops.

Methods
Plant material for genome sequencing. We sequenced an inbred line of Aeschynomene evenia C. Wright (evenia jointvetch) obtained by successive selfings from the accession CIAT22838. This accession was originally collected in Zambia and provided by the International Center for Tropical Agriculture (CIAT, Colombia) (http://genebank.ciat.cgiar.org). A. evenia was previously shown to be diploid (2n = 2x = 20) and to have a flow cytometry-estimated genome size of 400 Mb (1 C = 0.85 pg) 11,12,16 .
Genome sequencing and assembly into pseudomolecules. High-quality genomic DNA was prepared from the root tissue of 15-day-old plants cultured in vitro using an improved CTAB method 12 , followed by a high-salt phenol-chloroform purification according to the PacBio protocol. DNA was further purified using Ampure beads, quantified using the ThermoFisher Scientific Qubit Fluorometry, and fragment length was evaluated with the Agilent Tapestation System. A 20-kb insert SMRTbell library was generated using a BluePippin 15 kb lower-end size selection protocol (Sage Science). In all, 55 SMRT cells were run on the PacBio RS II system with collections at 4-hourly intervals and the P6-C4 chemistry 49 by the Norwegian Sequencing Center (CEES, Oslo, Norway). A total of 8,432,354 PacBio post-filtered reads was generated, producing 49 Gb of single-molecule sequencing data, which represented a 78× coverage of the A. evenia genome. PacBio reads were assembled using HGAP (version included in smrtpipe 2.3.0), the assembly was polished using the Quiver algorithm (SMRT Analysis v2.3.0) and then the SSPACE-LongRead (v1.1) program scaffolded the contigs when links were found (Supplementary Table 1). MiSeq reads were also generated to correct the sequence and estimate the genome size based on k-mer analysis (Supplementary Note 1). The de novo genome assembly contains 1848 scaffolds, with a scaffold N50 of 0.985 Mb and with 90% of the assembled genome being contained in 538 scaffolds. Then, we performed the A. evenia chromosomal-level assembly using serial analyses (fully described in Supplementary Note 1). The anchored scaffolds were joined with stretches of 100 Ns to generate 10 pseudomolecules named Ae01 to Ae10 according to the linkage group nomenclature for A. evenia 12 (Supplementary  Tables 3 and 4).
Gene prediction and annotation. First, repeats were called from the assembled genome sequence using RepeatModeler v1.0.11 (https://github.com/rmhubley/ RepeatModeler) (Supplementary Table 12). The genome was then masked using RepeatMasker v4-0-7 (http://www.repeatmasker.org/). Nine tissue-specific RNA-Seq libraries (sequenced by the GeT-PlaGe Platform, Toulouse, France) and fulllength transcripts generated from Iso-Seq (sequenced by the Cold Spring Harbor Laboratory, NY, USA) (details in Supplementary Note 1) were aligned on the unmasked reference with STAR 50 v2.7. The resulting BAM files were processed with StringTie 51 v1.3.3b to generate gene models in GTF format, which were merged with Cuffmerge from Cufflinks 52 v2.2.1 to produce a single GTF file. This GTF was used to extract a corresponding transcript FASTA file using the gtf_to_fasta program included in the TopHat 52 v2.0.14 package. The masked genome, the transcript fasta file, and the GFF files were used to train a novel AUGUSTUS 53 v3.2.3 model. This model was used to call the genes for all chromosomes. The AUGUSTUS prediction and the GTF files were then given to EVM 54 v1.1.1 to refine the model and remove wrongly called genes. This produced a new GFF file that was used to extract the corresponding transcripts using gtf_to_fasta. These transcripts were processed with TransDecoder 55 v2.1.0 in order to validate the presence of an open reading frame.
To check the completeness of the prediction, a master list of 100 nodulation genes was created and used for some additional manual annotation leading to the current annotation containing 32,667 gene models (Supplementary Table 7). Alignments of the Illumina RNA-seq clean reads from the nine samples with the STAR v2.7 software supported 25,301 of the 32,667 predicted genes (Supplementary Table 8). Finally, genome assembly and annotation quality was assessed using the Benchmarking Universal Single Copy Orthologs (BUSCO 56 v3) with the BLAST E-value cutoff set to 10 −5 (Supplementary Table 9). The BUSCO analysis includes a set of 1440 genes that are supposed to be highly conserved and single-copy genes present in all plants. Gene functions were assigned according to the best match of alignments using BLASP (1e-5) to SwissProt database. The InterPro domains, GO terms, and KEGG pathways database associated with each protein were computed using InterProScan with outputs processed using AHRD (Automated Human Readable Descriptions) (https://github.com/groupschoof/ AHRD) for selection of the best functional descriptor of each gene product (Supplementary Table 10).
Gene expression analysis. The normalized gene expression counts were computed using Cufflinks package based on the TopHat 51 output results of the RNA-Seq data analysis from the nine samples' analysis (Root N−, Root N+, Nodule 4d, Nodule 7, Nodule 14d, Stems, Leaves, Flowers, and Pods) performed for the A. evenia genome annotation. Gene expression was calculated by converting the number of aligned reads into FPKM (fragments per kilobase per million mapped reads) values based on the A. evenia gene models. RNA-seq data previously obtained from RNA samples of A. evenia IRFL6945 18 were also processed and converted into FPKM.
Orthogroup inference. We inferred orthogroups with OrthoFinder 57 v.0.4.0 to determine the relationships between A. evenia, the other diploid Aeschynomene taxa and several legume species. In the latter, proteomes were last obtained from the Legume Information System (https://legumeinfo.org/), the National Center for Biotechnology Center (https://www.ncbi.nlm.nih.gov), or from specific legume species websites in March 2020. They  Table 18). Phylogenies were created by aligning the protein sequences using MAFFT 58 v7.205 and genetic relationships were investigated in the trees generated with FastTree 59 v2.1.5 which is included in OrthoFinder. FigTree v1.4.3 (http://tree.bio.ed.ac.uk/) was subsequently used to further process the phylogenetic trees. A consensus species tree was also generated by OrthoFinder, based on alignment of single-copy orthogroups (i.e., an orthogroups with exactly one gene for each species).
Nodulation mutants. Nodulation mutants were obtained for A. evenia and characterized as fully described in the Supplementary Note 3. Briefly, a large-scale mutagenesis was performed by treating 9000 seeds from the CIAT22838 line with 0.30% EMS incubated overnight under gentle agitation. Germinated M 1 seedlings were transferred in pots filled with attapulgite. M 1 plants were allowed to self and 4-6 M 2 pods corresponding to approximately 40 seeds were collected from individual M 1 plants. Seeds collected from the same tray containing 72 M 1 plants were pooled and defined as one bulk. In all, 116 bulks of M 2 seeds were thus produced to constitute the EMS-mutagenized population. Phenotypic screening for nodulation alterations was conducted on 600 M 2 plants per bulk, 4 weeks after inoculation with the photosynthetic Bradyrhizobium strain ORS278. Plants with visible changes in their root nodulation phenotype were retained and allowed to self. The stability and homogeneity of the symbiotic phenotype was analyzed in the M 3 progeny. Whole inoculated roots of confirmed nodulation mutants were examined using a stereomicroscope (Niko AZ100; Campigny-sur-Marne, France) to identify alterations in nodulation and to establish phenotypic groups. The genetic determinism of the nodulation mutants was analyzed by backcrossing them to the CIAT22838 WT parental line according to the established hybridization procedure 11 and by determining the segregation of the nodulation phenotype in the F 2 population, 4 weeks post-inoculation with the Bradyrhizobium strain ORS278. These F2 plants were also used for additional analyses. Allelism tests were performed between selected nodulation mutants using the same crossing procedure 11 to define complementation groups.
Targeted sequence capture. For targeted sequence capture of symbiotic genes in nodulation mutants of A. evenia, 404 symbiotic genes known to be involved in the rhizobium-legume symbiosis or identified in expression experiments in A. evenia, were selected and their sequence extracted from the A. evenia genome to design custom baits with the following parameters: bait length 120 nucleotides, tiling frequency 2x. These probes were commercially synthetized by Mycroarray ® in a custom MYbaits kit (ArborBiosciences, https://arborbiosci.com/). DNA was extracted from roots of M 3 nodulation plants to construct genomic libraries using a preparation protocol developed at the GPTRG Facility of CIRAD (Montpellier, France) (Supplementary Note 3). The captured libraries were sequenced on an Illumina HiSeq 3000 sequencer at the GeT-Plage Facility of INRA (Toulouse, France) in 150 bp single-read mode. Read alignment and genome indexing were performed in the same way as for PoolSeq v0.3.3. Variations were called with Freebayes v1.1.0 with standard parameters and annotated according to their effect on A. evenia genes using SnpEff 62 (v4.3t and 'eff -c snpEff.config transcript' parameters). This file was then manually searched to identify the candidate gene variations able to explain the phenotypes.