RESEARCHFree Access icon

Transcriptomic Analysis With the Progress of Symbiosis in ‘Crack-Entry’ Legume Arachis hypogaea Highlights Its Contrast With ‘Infection Thread’ Adapted Legumes

    Authors and Affiliations
    • Kanchan Karmakar1
    • Anindya Kundu1
    • Ahsan Z Rizvi2
    • Emeric Dubois3
    • Dany Severac3
    • Pierre Czernic2
    • Fabienne Cartieaux2
    • Maitrayee DasGupta1
    1. 1Department of Biochemistry, University of Calcutta, Kolkata 700019, India;
    2. 2LSTM, Univ. Montpellier, CIRAD, INRA, IRD, SupAgro, Montpellier, France; and
    3. 3Montpellier GenomiX (MGX), c/o Institut de Génomique Fonctionnelle, 141 rue de la cardonille, 34094 Montpellier Cedex 05, France


    In root-nodule symbiosis, rhizobial invasion and nodule organogenesis is host controlled. In most legumes, rhizobia enter through infection threads and nodule primordium in the cortex is induced from a distance. But in dalbergoid legumes like Arachis hypogaea, rhizobia directly invade cortical cells through epidermal cracks to generate the primordia. Herein, we report the transcriptional dynamics with the progress of symbiosis in A. hypogaea at 1 day postinfection (dpi) (invasion), 4 dpi (nodule primordia), 8 dpi (spread of infection in nodule-like structure), 12 dpi (immature nodules containing rod-shaped rhizobia), and 21 dpi (mature nodules with spherical symbiosomes). Expression of putative ortholog of symbiotic genes in ‘crack entry’ legume A. hypogaea was compared with infection thread–adapted model legumes. The contrasting features were i) higher expression of receptors like LYR3 and EPR3 as compared with canonical Nod factor receptors, ii) late induction of transcription factors like NIN and NSP2 and constitutive high expression of ERF1, EIN2, bHLH476, and iii) induction of divergent pathogenesis-responsive PR-1 genes. Additionally, symbiotic orthologs of SymCRK, ROP6, RR9, SEN1, and DNF2 were not detectable and microsynteny analysis indicated the absence of a RPG homolog in diploid parental genomes of A. hypogaea. The implications are discussed and a molecular framework that guides crack-entry symbiosis in A. hypogaea is proposed.

    Nitrogen-fixing root-nodule symbiosis (RNS) allows plants to house bacterial diazotrophs in an intracellular manner (Kistner and Parniske 2002). The establishment of RNS involves rhizobial invasion in the root epidermis and nodule organogenesis in the root cortical cells. The most common invasion strategy is through root-hair curling and infection thread (IT) formation, in which the nodule primordia are induced from a distance (Gage 2004). Invasion through the IT is widespread among temperate legumes like Vicia, Trifolium, and Pisum species. Model legumes like Lotus japonicus and Medicago truncatula also undertake IT-mediated rhizobial invasion (Geurts and Bisseling 2002; Oldroyd and Downie 2004, 2006). The alternate mode of rhizobial invasion is known as ‘crack entry’, in which the rhizobia enter through natural cracks at the lateral root base in an intercellular manner. Approximately 25% of all legume genera are adapted to crack-entry and it is a characteristic feature of some subtropical legumes belonging to dalbergoid/genistoid clades, like Arachis, Aeschynomene, and Stylosanthes species (Sprent and James 2007). In these legumes, rhizobia directly access the cortical cells for development of their nodule primordia and the infected cells repeatedly divide to develop the mature nodule with an aeschynomenoid infection zone (IZ) devoid of uninfected cells (Boogerd and Rossum 1997; Fabre et al. 2015). Very limited information is available for crack-entry legumes from the dalbergioid/genistoid clade that are basal in their divergence within the papilionoideae.

    Investigations on model legumes that are adapted to IT-mediated rhizobial invasion have unraveled the molecular basis of RNS. The host responses are initiated by Nod factor (NF) receptors (NFRs) LjNFR1/MtLYK3 and LjNFR5/MtNFP (Arrighi et al. 2006; Madsen et al. 2003; Radutoiu et al. 2003; Smit et al. 2007). Another NF-induced receptor LjEPR3 was shown to monitor rhizobial exopolysaccharide (EPS) in L. japonicus, indicating a two-stage mechanism involving sequential receptor-mediated recognition of NF and EPS signals to ensure host symbiont compatibility (Kawaharada et al. 2015). Downstream to NFRs is the ‘SYM pathway’ consisting of the receptor kinase LjSYMRK/MtDMI2 (Endre et al. 2002; Stracke et al. 2002), the predicted ion-channel proteins LjCASTOR and LjPOLLUX/MtDMI1(Ané et al. 2004; Imaizumi-Anraku et al. 2005), the nucleoporins LjNUP85 and LjNUP133(Kanamori et al. 2006; Saito et al. 2007), the Ca2+/calmodulin-dependent protein kinase LjCCaMK/MtDMI3 (Lévy et al. 2004; Tirichine et al. 2006), and the transcription factor (TF) LjCYCLOPS/MtIPD3(Messinese et al. 2007; Singh et. al. 2014; Yano et al. 2008). Nodulation-specific TFs, such as MtNSP1/LjNSP1, MtNSP2/LjNSP2, MtERF1, and MtNIN/LjNIN function downstream of the SYM pathway and are involved in transcriptional reprogramming for initiation of RNS (Kaló et al. 2005; Middleton et al. 2007; Schauser et al. 1999; Smit et al. 2005).

    Transcriptome analyses during symbiosis in several legumes like M. truncatula, L. japonicus, Glycine max, and Cicer arietinum provided molecular insights into symbioses that use the ‘root hair’ mode of entry (Limpens et al. 2014; Kant et al. 2016; Takanashi et al. 2012; Yuan et al. 2017). These resources have been useful in determining gene expression associated with different stages of RNS, for example, symbiont recognition, onset of infection, meristem activity for nodule organogenesis, and functional differentiation of infected cells (Breakspear et al. 2014; Roux et al. 2014). The outcome, in most cases, highlighted numerous novel genes coexpressed with known markers for symbiosis. Nodule transcriptome analyses in M. truncatula highlighted that TF genes essential for the control of the root apical meristem were also expressed in the nodule meristem (Roux et al. 2014). Transcriptomic analysis of infected root hairs further revealed involvement of genes regulating cell cycle, flavonoid biosynthesis, and auxin signaling during IT initiation (Breakspear et al. 2014). It also enabled identification of genes and processes previously undetected in whole-root studies or in forward genetic analyses. In M. truncatula transcriptomic analysis of symbiotic mutants further improved our understanding of gene expression associated with specific processes of symbiosis (Jardinaud et al. 2016; Moreau et al. 2011). But unlike in model legumes, transcriptomic profiling for species that utilize crack entry has rarely been reported to date and the symbiotic pathway of crack entry nodulation is largely unknown.

    An earlier report listed several differentially expressed genes (DEGs) at an early stage of symbiosis in a nodulation-competent A. hypogaea in comparison with an incompetent variety (Peng et al. 2017). Our manuscript undertakes a systematic effort to study the transcriptomic dynamics with the onset and advancement of symbiosis at five different nodule developmental stages in A. hypogaea, using uninfected roots (UI) as a reference. Primarily, we focused on symbiotic genes that were characterized in model legumes and could be broadly classified according to their functional involvement in bacterial recognition, early signaling, infection, nodule organogenesis, functional differentiation, and nodule number regulation. Expression of putative orthologs of these symbiotic genes in crack-entry legume A. hypogaea is compared with the corresponding expression profiles in M. truncatula and L. japonicus that undertake IT-mediated symbiosis. And the outcome of the analysis is summarized in a model that can be further utilized for understanding nodulation dynamics in basal legumes.


    Progress of symbiosis in Arachis hypogaea.

    Within 3 weeks of infection with Bradyrhizobium sp. strain SEMIA 6144, A. hypogaea roots developed spherical functional nodules. We followed the progress of symbiosis in A. hypogaea for 21 days, to identify the distinct stages of development by ultrastructure analysis. There are rosettes of root hairs in the junction of taproot and lateral root that are reported to be important for bacterial invasion in A. hypogaea (Boogerd and Rossum 1997). Within 1 day postinfection (dpi), rhizobia was found to be adhered to these root hairs (Fig. 1A to C). Within 4 dpi, bump-like primordial structures were noted at the lateral root bases (Fig. 1D). The longitudinal sections of these primordia revealed one or more centrally located defined pockets of rhizobia-infected cells that were surrounded by uninfected cells (Fig. 1E). These pockets of rhizobia-infected cells were distinct by having reduced Calcofluor-binding ability, indicating that they are thin-walled. The intracellularized rhizobia within the infection pockets were undifferentiated and rod-shaped (Fig. 1F). The infection pockets observed at 4 dpi act as IZ founder cells, and it is their uniform division and differentiation that give rise to the distinct aeschynomenoid-type IZ in mature nodules. There has not been a single case in which uninfected primordium was noted, which is in accordance with the proposition of infection preceding development of aeschynomenoid nodules (Fabre et al. 2015). By 8 dpi, there was a visible nodule-like structure at the lateral root base (Fig. 1G). Ultrastructure analysis revealed that, by 8 dpi, the compactness of the primordial structure with defined pockets of infected cells was lost and the IZ started growing by division of the infected cells (Fig. 1H and I). By 12 dpi, there were white spherical nodules (Fig. 1J). At this stage, the tissue organization turned aeschynomenoid where there were no uninfected cells in the IZ and the endocytozed rhizobia remained undifferentiated and rod-shaped (Fig. 1K and L). At 21 dpi, the nodules were mature and functional where the rhizobia differentiated within the plant-derived peribacteroid membranes to develop spherical symbiosomes (Fig. 1M to O).

    Fig. 1.

    Fig. 1. Progress of symbiosis in Arachis hypogaea infected with Bradyrhizobium sp. strain SEMIA6144. A, Root was harvested at 1 day postinfection (dpi) at the root junction. B, Scanning electron microscopy of inoculated root hair at lateral root junction and C, confocal laser scanning microscopy of lateral root junction. D to F, Section of nodule primordia at 4 dpi (dashed lines indicate the infection zone in E and asterisks indicate the identical position in E and F) and G to I, 8 dpi (asterisks indicate the identical position in H and I). J to L, Section of nodules at 12 and M to O, 21 dpi. A, D, G, J, and M, stereoimages; C, E, H, K, and N, bright field and propidium iodide (PI) merged; F, I, and L, PI + Calcofluor merged; and O, PI + Calcofluor + syto9 merged. Arrows indicate rod shaped rhizobia in F, I, and L. Scale bar = 500 μm (A, D, G and J), 1 mm (M), 2 μm (B), 100 μm (E, H, K, and N), and 10 μm (C, F, I, L, and O).

    Download as PowerPoint

    Transcriptome analysis with the progress of symbiosis in Arachis hypogaea.

    Ultrastructural analysis revealed five distinct stages during the progress of symbiosis in A. hypogaea 1 dpi (recognition and invasion), 4 dpi (primordia formation), 8 dpi (nodule-like structure), 12 dpi (immature nodules with rod-shaped rhizobia), and 21 dpi (mature nodules with spherical symbiosomes). To probe into the expression of genes associated with the progress of symbiosis, RNA was extracted from these stages. As a reference (0 dpi), the UI roots were used. RNA-seq was done in triplicate for these six stages, using Illumina single-end sequencing technology (Illumina Hiseq 2000 SR50). The genomic data from Arachis duranensis (AA) and Arachis ipaensis (BB), which are two wild diploid parents of A. hypogaea, were used to assess the quality and coverage of the assembled transcriptomes. A total of 1,429,876,614 raw reads of 50 bp (approximately 71.5 Gb) were generated with an average of 88,029,386 reads per library. This was 600 times the total size of transcript sequences (109.0 Mb) of A. hypogaea for both AA and BB genomes and gave an average coverage of 36 times per library. The proportion of clean reads among the total acquired reads was more than 91.34% (Table 1). The filtered reads were simultaneously mapped to the AA and BB genomes, in which the overall accepted mapping rate per library ranged from 80.15 to 89.98%, with an average mapping rate of 86.42% with A. duranensis (AA) and 86.65% with A. ipaensis (BB). For both AA and BB genomes, about 66% of reads aligned to a gene exon in an unambiguous way, whereas the remaining 33% of reads aligned outside exons.

    Table 1. Summary of raw Illumina sequencing and filtered reads after trimming and alignment of reads to AA (Arachis duranensis) and BB (A. ipaensis) genomes in each library

    The expression level of each assembled transcript sequence was measured through FPKM (fragments per kilobase per million reads) values. The DEGs in the five different stages of symbiosis were evaluated by the significance of differences in their expression with respect to UI roots, using a false discovery rate (FDR) < 0.05, P value < 0.05, and fold change log2 ratio ≥1. Since a single UI was used for all the stages, there could be representations of some DEGs that are associated with normal root development within the experimental timeframe. Comparison between upregulated and downregulated DEGs at different stages is shown in a Venn diagram in Figure 2A. A total of 2,745 genes were up-regulated (1,297 in AA, 1,448 in BB) and a total 20,259 genes are down-regulated (9,552 in AA, 10,707 in BB), with the progress of symbiosis in A. hypogaea of which 59 genes (29 in AA, 30 in BB) were commonly up-regulated and 2,054 genes (1,015 in AA, 1,039 in BB) were commonly down-regulated in all five of the stages of symbiosis. From the Venn diagram, we identified those genes that were first up- or down-regulated at a particular stage, though their subsequent regulations could be different. Differentially expressed genes show clear transcriptional shifts at these stages, and the diverse expression patterns of these genes are indicated in a heatmap (Fig. 2B). Hierarchical clustering (Fig. 2B) and principal component analysis (PCA) both suggested that there is substantial overlap between the transcriptome of infected roots at 1 dpi and at 4 dpi, when rhizobial invasion and primordia formation occurs. Similarly, there was overlap between transcriptome of immature nodules at 8 and 12 dpi, when the primordia structurally develop into a nodule. Finally, the transcriptome from mature nodules at 21 dpi was distinct from all other stages in which a nodule matures to its functional form. These expression waves indicated that there could be three major transcriptional programs during the inception and maturation of symbiosis in A. hypogaea. The actual number of genes upregulated or downregulated from either the AA or BB genome at 1 dpi was 337 and 8,046 (156 upregulated and 3,826 downregulated in AA;181 upregulated and 4,219 downregulated in BB), at 4 dpi was 388 and 4,643 (180 upregulated and 2,167 downregulated in AA; 208 upregulated and 2,476 downregulated in BB), at 8 dpi was 641 and 1,587 (297 upregulated and 735 downregulated in AA; 344 upregulated and 852 downregulated in BB), at 12 dpi was 502 and 1,548 (240 upregulated and 751 downregulated in AA; 262 upregulated and 797 downregulated in BB), and for mature nodules at 21 dpi was 877 and 4,435 (423 upregulated and 2,141 downregulated in AA; 454 upregulated and 2,294 downregulated in BB). This analysis shows that, within 1 to 4 dpi, the number of downregulated DEGs were approximately 18-fold higher than the number of upregulated DEGs, whereas in other stages the difference was approximately three- to fivefold. Previous reports on symbiosis-associated transcriptome of other legumes also indicated the prevalence of a higher number of downregulated transcripts (Boscari et al. 2013; Yuan et al. 2017). But the significantly high number of downregulated transcripts in early stages of nodule inception could be a unique feature for A. hypogaea RNS.

    Fig. 2.

    Fig. 2. Comparison of differentially expressed genes (DEGs) at five different stages during the progress of symbiosis in Arachis hypogaea. The total number of DEGs at 1, 4, 8, 12, and 21 days postinfection (dpi) was identified by EdgeR analysis from both AA and BB genomes. A, Venn diagram showing comparison between upregulated and downregulated DEGs. B, Heat map of the hierarchical cluster analysis of the same DEGs. Columns represent different timepoints and the arborescence above them indicates similarity. Scale bar represents log2-fold change in gene expression. Numbers beside the heat map indicate the number of DEGs at different timepoints.

    Download as PowerPoint

    Functional analysis of DEGs.

    Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) categories were searched for terms that are significantly enriched in our DEGs. Among the 1,248 enriched GO terms, there was a major representation of defense response genes; 470 and 31 such defense-related GO terms were enriched in down- and upregulated DEGs, respectively. Accordingly, KEGG analysis of plant-pathogen interaction pathways show that most genes involved in pattern-triggered immunity was notably down-regulated (Fig. 3A and B). The flagellin-sensing receptor FLS2 mediated the mitogen-activated protein kinase pathway (MAPK), however remained active along with a subset of cyclic nucleotide-gated ion channel genes (CNGCs) and genes encoding respiratory burst oxidase (Rboh) proteins. A subset of genes involved in the effector-triggered immunity also remained active during symbiosis, for example, genes encoding resistance proteins, like RPM1, RPS2, RPS5, Pti1 kinase, and pathway regulators, like SGT1, HSP90, and EDS1. Intriguingly there was a significant upregulation of genes encoding pathogenesis-responsive (PR)-1 proteins, which are members of the CAP (cysteine-rich secretory proteins, antigen 5, and pathogenesis-related 1 proteins) superfamily (Breen et al. 2017). The PR-1 proteins upregulated during symbiosis clustered away from the PR-1 proteins that were reported to be upregulated in defense responses, indicating the symbiosis-associated PR-1 proteins were divergent in nature (Fig. 3C). There are two PR-1 proteins that clustered with defense-responsive PR proteins and these PR-1 genes were not upregulated during symbiosis, further confirming the symbiotic PR-1 proteins to be distinct. PR-1 proteins harbor an embedded defense signaling peptide (CAP-derived peptides or CAPE) in which CNYxPxGNxxxxxPY is considered as a functional motif that marks cleavage of these bioactive peptides (Breen et al. 2017). The cleavage site is conserved in both classes of PR-1 proteins, suggesting the CAPE peptides could also be generated from the divergent PR-1 proteins synthesized during symbiosis (Fig. 3C and D). Since genes encoding CAP proteins are marker genes for the salicylic acid (SA) signaling pathway and systemic acquired resistance, we also checked the SA/JA (jasmonic acid) pathways to further understand the symbionts responsive signaling in A. hypogaea. As shown in Figure 3A, the JA pathway was completely downregulated but the SA-responsive genes, like TGA1 and NPR1, were up-regulated. Thus, symbiotic PR-1 gene expression could be justified by the activation of the SA-mediated signaling.

    Fig. 3.

    Fig. 3. Differentially expressed genes (DEGs) associated with plant-pathogen interaction and jasmonic acid and salicylic acid signaling pathways. A, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of upregulated and downregulated DEGs in the respective pathways. B, Heat map of the DEGs that are annotated in KEGG pathway analysis. C, Neighbor-joining phylogenetic tree of all the annotated CAP (cysteine-rich secretory proteins, antigen 5, and pathogenesis-related 1 proteins), using nontruncated amino acid sequences when a branch denotes divergent upregulated CAP-PR1 proteins. D, CLUSTALW sequence alignment of CAPE (CAP-derived) peptides using Multalin; an arrow indicates cleavage site for CAPE peptide; CBM and CAP domain are annotated by shaded boxes, and bullets indicate Arachis hypogaea CAP peptides.

    Download as PowerPoint

    Several genes are reported to be expressed in nodulating roots by comparing the transcriptome profiles of nonnodulating and nodulating lines of A. hypogaea (Peng et al. 2017). The list includes known symbiotic genes, like NIN and CLE13, and other genes, like NF-YA10, Myb127, a receptor kinase, a soluble kinase, a F-BOX protein, TFs of the SHI family, and a lectin. All these genes were represented in our upregulated transcriptome, which thereby revalidates the importance of expression of these genes during the onset of symbiosis in A. hypogaea.

    Expression profiles of putative orthologs of symbiotic genes.

    Our final objective was to understand the expression of the putative orthologs of symbiotic genes in A. hypogaea that are characterized in model legumes like M. truncatula and L. japonicus. A total of 71 such genes were chosen and were functionally classified on the basis of their demonstrated role during nodulation in model legumes (Fig. 4). A BLAT search on A. ipaensis and A. duranensis genomes identified 63 genes that were annotated in PeanutBase and six genes that are yet to be annotated, i.e., NSP2, ERN1, ERF1, VPY, ENOD40, and RSD. A homolog for MtRPG could not be detected by BLAT search in either of the two parental genomes of A. ipaensis and A. duranensis, which could be due to an incomplete genome sequence or the quality of the gene model predictions. Microsynteny analysis in Phytozome indicated the absence of MtRPG orthologs in A. ipaensis, whereas, in A. duranensis, the entire gene block corresponding to MtRPG was not detectable. In a recent report, Griesmann et al. (2018) also indicated the loss of RPG gene in the A. ipaensis genome. Similar microsynteny analysis also revealed the absence of MtDNF2 in A. ipaensis genome but indicated the presence of a singleton in A. duranensis. This singleton has 16% amino acid sequence homology with MtDNF2 and phylogenetic analysis (maximum likelihood) reveals that the A. duranensis sequence is highly divergent from both legume and nonlegume sequences. MtDNF2 codes for a phosphatidyl-inositol phospholipase C-like protein, but a BLAST search indicated this singleton to code for an unknown protein. While the singleton showed very low transcript density, we could not identify any homologous contigs corresponding to MtRPG and MtDNF2 in our RNA-seq data, using Burrow-Wheeler Aligner (Fig. 4G).

    Fig. 4.

    Fig. 4. Comparative expression pattern of 71 symbiotic genes in Arachis hypogaea, Medicago truncatula, and Lotus japonicus. For A. hypogaea, histograms represent normalized RNA-seq reads (fragments per kilobase per million reads) of symbiotic orthologous genes aligned to A. duranensis and A. ipaensis genomes. For M. truncatula and L. japonicus, histograms represent microarray data retrieved from the respective Gene Expression Atlas Affymetrix database (i.e., MtGEA and LjGEA) (Benedito et al. 2008; Verdier et al. 2013). Expression values of symbiotic genes (in arbitraty units [a.u.]) during nodulation are grouped as A, bacterial recognition, B, SYM pathway and early signaling, C, interactors, D, early transcription factors, E, infection, F, nodule organogenesis, G, nodule differentiation, H, nodule number regulation (dpi: days postinfection by Bradyrhizobium sp. strain SEMIA6144 for A. hypogaea, Sinorhizobium meliloti for M. truncatula, and Mesorhizobium loti R7A for L. japonicus). a, b, c, d, e, and f in the x axis represent 0, 1, 4, 8, 12, and 21 dpi, respectively for A. hypogaea; 0, 3, 6, 10, 14, and 20 dpi, respectively, for M. truncatula; 0, 1, 3, 7, 14, and 21 dpi, respectively, for L. japonicus. An asterisk indicates genes that are differentially expressed in A. hypogaea; a diamond indicates expression of a nonsymbiotic ortholog in A. duranensis and A. ipaensis, and an X indicates a gene absent in the diploid genome.

    Download as PowerPoint

    Orthology for the identified symbiotic genes were checked by reciprocal BLAST and sequence alignments. Symbiotic orthologs were identified for 65 of 70 genes, for which the A. hypogaea homologs clustered with other legumes in the corresponding gene trees. Exceptions were the homologs for MtSymCRK, LjROP6, MtRR9, and LjSEN1, for which the A. hypogaea genes clustered with nonlegumes, and MtDNF2, which was highly divergent and did not cluster with either legumes or nonlegumes. For comparison of expression between the symbiotic or nonsymbiotic orthologs in Arachis sp. with their counterpart in model legumes, we used the microarray data derived from the M. truncatula gene expression atlas (MtGEA) (Benedito et al. 2008) and L. japonicus gene expression atlas (LjGEA) (Verdier et al. 2013). For the symbiotic genes that are characterized from one of these model legumes, reciprocal BLAST was done to identify the ortholog in the other. The pattern of expression of these genes with the progress of symbiosis was compared between A. hypogaea and model legumes. Both absolute expression (Fig. 4) and relative (log2-fold) expression values were analyzed so that high and constitutively expressed genes were not ignored. The notable exceptions are mentioned below, and the analyzed data are summarized in a model in Figure 5. The contrasting genes are highlighted along with the genes that are absent or divergent in sequence.

    Fig. 5.

    Fig. 5. Simplified model for the symbiotic signaling pathway in Arachis hypogaea. Symbiotic genes identified in Medicago truncatula and Lotus japonicus are listed and classified according to their main symbiotic functions. The homologs of these genes in A. hypogaea that are distinct from model legumes are marked as follows: an X (absent in diploid genome), shaded diamonds (symbiotic ortholog absent), an upward arrow (upregulation), and black diamonds (divergent expression).

    Download as PowerPoint

    In the recognition module, expression of lipo-chitooligosaccharide (LCO)-binding AhLYR3 and EPS-binding AhEPR3 were significantly high as compared with the orthologs of canonical NF LysM domain-containing receptors AhNFR1 and AhNFR5 (Figs. 4A and 5). Most of the SYM pathway members and their interactors and the early TFs showed similar trends in expression in both A. hypogaea and in model legumes for both inducible and constitutive genes (Fig. 4B to D). Exception was the gene encoding an ortholog of cyclic nucleotide-gated channel MtCNGC (Charpentier et al. 2016), which was significantly upregulated in A. hypogaea (Fig. 4B). Among the interactors, the ortholog of an E3 ubiquitin ligase, MtPUB1, was constitutively expressed in A. hypogaea, whereas in the model legumes its expression is up-regulated only upon bacterial invasion (Fig. 4C). In model legumes, TFs like NIN and NSP2 are upregulated with bacterial invasion in early stages (Kaló et al. 2005; Schauser et al. 1999). In contrast, upregulation of AhNIN was noted at 8 dpi onwards, and AhNSP2 expression was induced only at 21 dpi in mature nodules (Figs. 4D and 5). The expression of ethylene responsive TF, AhERF1 appeared to be very high and constitutive with the progression of symbiosis. On the other hand, in model legumes, the expression of ERF1 was significantly reduced (L. japonicus) or upregulated (M. truncatula) upon bacterial invasion.

    Among the genes associated with rhizobial invasion, the homolog for MtRPG (Arrighi et al. 2008), a factor responsible for rhizobium-directed polar growth of Its, was absent in A. hypogaea (Fig. 5). Additionally, symbiotic orthologs of LjROP6, a Rho-like GTPase that is responsible for IT elongation in the cortex was not detected in Arachis sp., though the expression of its nonsymbiotic homolog was significantly high during symbiosis (Fig. 4C). Symbiotic orthologs of flotillins-like MtFLOT2 and MtFLOT4 (Haney and Long 2010; Ke et al. 2012), which are required for lipid raft formation during infection, were present in A. hypogaea, but only the expression of AhFLOT2 was detectable during symbiosis, suggesting its adaptation during crack-entry nodulation (Fig. 4E). The expression of the ortholog of LjCERBERUS, a WD-40 repeats U-box protein essential for the progress of infection in L. japonicus (Yano et al. 2009), was significantly high in A. hypogaea as compared with its expression pattern in model legumes (Fig. 4E).

    Among the factors associated with nodule development, the transcript for MtDNF2 homolog encoding a phosphatidyl-inositol phospholipase C-like protein (Bourcy et al. 2013), was not detectable in A. hypogaea, but syntenic analysis revealed a highly divergent copy of DNF2, represented by a singleton in our RNA-seq data from a progression of symbiotic samples (Fig. 4G). Using the Burrow-Wheeler Aligner, we found that the singleton showed very low transcript density in our RNA-seq data, and intriguingly, PeanutBase also does not show any expression from this locus. In addition, symbiotic orthologs of non-RD receptor kinase MtSymCRK (Berrabah et al. 2014), type A response regulator MtRR9, and metal transporter LjSEN1 was also not detectable in A. duranensis and A. ipaensis genomes (Fig. 4F and G). But significant expression of the nonsymbiotic homologs of these genes was noted during the progress of symbiosis in A. hypogaea. On the other hand, the expression of orthologs for cytokinin inducible TF MtbHLH476 (Ariel et al. 2012) and the ethylene-responsive TF MtEIN2 (Penmetsa et al. 2008) (sickle) was found to be significantly higher in A. hypogaea as compared with their expression pattern in model legumes (Figs. 4F and H and 5). Apart from the indicated differences, relative expression patterns of all other symbiotic genes were found to be comparable between A. hypogaea and model legumes during progress of symbiosis, indicating their conserved roles.

    Quantitative reverse transcription-polymerase chain reaction (qRT-PCR) validation of the RNA-seq results.

    To validate the RNA-seq results, we did qRT-PCR analysis using AhActin as the reference gene. The stability of expression of AhActin during symbiosis in A. hypogaea was reported earlier (Kundu and DasGupta 2018; Morgante et al. 2011). The expression was evaluated for 10 symbiotic genes from different modules that are described in Figure 5. Additionally, expression was evaluated for five more genes that were described to be important for A. hypogaea symbiosis (Peng et al. 2017). The results were in agreement with the transcriptional profile data for 83 of 90 (92%) datapoints. The fold change of expression did not always match between RNA-seq and qRT-PCR, but the trend of expression was similar in most of the timepoints for all the genes. No change of expression was noted for any of these genes in uninoculated roots corresponding to different infection timepoints, further confirming the expression of these genes to be specifically associated with symbiosis.

    PCA of symbiotic gene expression in Arachis hypogaea and model legumes.

    PCA was done to check if there was a signature, in the expression pattern of symbiotic genes in crack-entry legume A. hypogaea, that contrasts with the model legumes in which rhizobial entry is IT-mediated. Figure 6 is a projection of differential expression of symbiotic genes from A. hypogaea, M. truncatula, and L. japonicus for all the infection timepoints into the first two principal components. Altogether, expression of around 87% of the genes were accommodated along dimension 1 (PC1) and dimension 2 (PC2) in the analysis. Expressions of symbiotic genes that are constitutive or show minimal change are clustered near the origin for both Arachis sp. and the model legumes. These include the SYM pathway and most early signaling genes that are likely to be regulated at the posttranscriptional level. The genes that were noted in Figure 4 to have a contrasting trend in expression between A. hypogaea and the model legumes were all placed in opposing or adjacent quadrants in the PCA. For example, AhCNGC, AhNIN, AhNSP2, AhPUB1, and AhCERBERUS were clustered separately from their counterparts in model legumes. On the other hand, there are some factors whose contrasting trends in expression between A. hypogaea and model legumes was not obvious in Figure 4, but PCA highlighted their differences. For example, AhVPY, AhENOD40, AhWOX5, AhRR1, and AhRR4 were placed in opposing or adjacent quadrants from their counterparts in model legumes. These factors highlighted by PCA can also be considered to be differentially adapted for crack-entry mediated rhizobial invasion and aeschynomenoid nodule development in A. hypogaea.

    Fig. 6.

    Fig. 6. Principal component analysis (PCA) of differentially expressed symbiotic genes. Principal components of the variability in expression values for the symbiotic genes at all the infection timepoints in Arachis ipaensis, Medicago truncatula, and Lotus japonicus was calculated using R-Package. The two principal components and their fraction of the overall variability of the data (%) are shown on the x axis (PC1 68.40%) and the y axis (PC2 19.70%).

    Download as PowerPoint


    This is the first systematic effort toward transcriptome profiling with the progress of symbiosis in a crack-entry legume, A. hypogaea. We could correlate the hierarchical shift in upregulated DEGs to the morphological changes occurring during different timepoints of nodule development, and the process appeared to be governed by three major transcriptional programs. The first program is associated with rhizobial recognition and generation of nodule primordia by 1 to 4 dpi, the second program is associated with structural development of nodules by 8 to 12 dpi, and the third program is associated with functional maturation of nodules at 21 dpi (Figs. 1 and 2). The comparison between expression patterns of putative orthologs of symbiotic genes with progress of symbiosis in A. hypogaea with model legumes highlighted the genes that are important or dispensable for the crack entry–mediated root nodule symbiosis (Figs. 4, 5, and 6). Additionally, absence of symbiotic orthologs of certain symbiotic genes in A. hypogaea was highlighted. The probable implications of these distinguishing features in A. hypogaea are discussed in the light of their importance during crack-entry nodulation.

    The most significant observation in the A. hypogaea symbiotic transcriptome was the significantly high expression of a group of genes encoding a divergent form of cysteine-rich PR-1 proteins (Fig. 3). Upregulation of AhCNGC and AhERF1 in the symbiotic transcriptome (Fig. 4) could explain the expression of PR1 genes, as their role in SA-induced PR gene expression is well-documented (Asamizu et al. 2008; Moeder et al. 2011; Solano et al. 1998). PR-1 proteins are ubiquitous across plant species and are among the most abundantly produced proteins in plants in response to pathogen attack (Breen et al. 2017). As a result, it is used as a marker for SA-mediated disease resistance in plants. Although differential expression of defense response genes belonging to GO:0006952 (defense related) and PR-1 and PR-10 protein families has previously been reported for M. truncatula RNS (Jardinaud et al. 2016), this is the first report in which a divergent group of PR-1 proteins is shown to be associated with nodule development (Fig. 3). PR-1 proteins harbor a caveolin binding motif (CBM) that binds sterol and the embedded Prorich C-terminal peptide (CAPE) is involved in plant immune signaling (Breen et al. 2017). All the symbiotic PR-1proteins in A. hypogaea have both these conserved features, but whether these CAPE peptides are actually derived from PR-1 proteins during symbiosis remains to be understood. It may be noted that, in the IRLC clade of legumes as well as in the dalbergioid legume Aeschynomene evenia, another group of cysteine-rich peptides belonging to the NCR family is responsible for bacterial differentiation associated with morphological change of symbionts (Van de Velde et al. 2010). Intriguingly, morphological change of symbionts from rod to spherical shape is also noted in A. hypogaea (Kundu and DasGupta 2018). Whether the antimicrobial CAPE peptides are enrolled as symbiosis effectors in A. hypogaea in a manner similar to NCRs remains to be understood.

    Bacterial recognition is NF-dependent in A. hypogaea and NF-independent for Aeschynomene evenia, though rhizobial invasion in both these legumes is mediated through epidermal cracks (Chandler et al. 1982; Giraud et al. 2007; Ibáñez and Fabra 2011). In accordance, the canonical NF receptor NFR1 was absent in Aeschynomene evenia (Chaintreuil et al. 2016) but was expressed in A. hypogaea (Figs. 4A and 5). Intriguingly, the expression of AhLYR3, ortholog of LCO-binding LysM receptor MtLYR3 (Fliegmann et al. 2016), was higher in A. hypogaea as compared with canonical NFR1/5, which is unlike model legumes, in which the expression of canonical NFR genes were much higher than orthologs of MtLYR3 (Fig. 4A). Since LjLYS12, an ortholog of MtLYR3, is important for recognition of pathogens (Kelly et al. 2017), it is tempting to suggest a possible involvement of AhLYR3 in symbiont recognition in A. hypogaea. PUB1, a substrate of LjNFR1, has a role in determining partner specificity (Mbengue et al. 2010). The constitutive expression of AhPUB1 in A. hypogaea suggests that it might have such a role during the choice of symbionts (Fig. 4C). Another step in bacterial recognition is through perception and recognition of surface molecules present on the bacterial cell wall, like EPS, LPSs, cyclic glucans, and capsular polysaccharides (Via et al. 2016). EPS-binding LysM receptor MtEPR3 has a role in bacterial recognition followed by IT elongation (Kawaharada et al. 2015), and its ortholog in A. hypogaea showed significant expression with progress of symbiosis (Fig. 4A). This receptor in Aeschynomene evenia transcriptome (Chaintreuil et al. 2016), indicating that it may not be important during crack-entry but, rather, it might have a role in symbiont recognition during A. hypogaea nodulation. Thus, similar to IT-adapted model legumes, expression of LysM receptors may also be important for symbiont recognition during intercellular invasion of rhizobia in a crack-entry legume like A. hypogaea.

    Most of the genes that are important for IT progression in model legumes are significantly expressed during the progress of symbiosis in the crack-entry legume A. hypogaea (Fig. 4E). An apparent exception is MtRPG, which has a role in polar growth of ITs in M. truncatula (Arrighi et al. 2008; Murray et al. 2011). Our results indicate the absence of a RPG homolog in Arachis sp., which is in accordance with a recent report that loss of RPG gene is indicated in the A. ipaensis genome (Griesmann et al. 2018). Since RPG is present in the nonlegumes, it is expected to be present in the legume ancestor that has been selectively lost in A. hypogaea. Whether this loss of RPG precedes or follows the evolution of nodulation in A. hypogaea remains to be understood. Along with the absence of RPG, the symbiotic ortholog of MtFLOT4, which is important for IT initiation and elongation (Haney and Long 2010), has very low expression in A. hypogaea (Fig. 4E) and was absent in Aeschynomene evenia (Chaintreuil et al. 2016). Other factors that are important for IT progression in model legumes are MtFLOT2 (Haney and Long 2010), MtVPY (Murray et al. 2011), LjCERBERUS (Yano et al. 2009), and LjARPC1 (Hossain et al. 2012). Symbiotic orthologs for all these genes are present in A. hypogaea and they demonstrated similar expression pattern during symbiosis in comparison with the IT-adapted model legumes, indicating the importance of their additional roles during crack-entry nodulation (Fig. 4). For example, LjCERBERUS is functionally associated with nodule organogenesis and LjARPC1 with cytoskeletal rearrangement during bacterial endocytosis (Hossain et al. 2012; Yano et al. 2009) and these additional functions of CEREBERUS and ARPC1 could be conserved in crack-entry legumes.

    Though cytokinin signaling is important for nodule development in both crack-entry and IT-adapted legumes (Fabre et al. 2015; Gonzalez-Rizzo et al. 2006; Kundu and DasGupta 2018), PCA revealed a distinct expression pattern for several downstream factors in A. hypogaea as compared with model legumes. For example, a homolog of MtRR1, a type B response regulator (RR) has a distinct expression pattern in A. hypogaea and, unlike model legumes, its downstream TF, MtbHLH476 (Ariel et al. 2012), was highly expressed with the progress of infection in A. hypogaea (Figs. 4F and 6). The symbiotic ortholog of MtRR9 (Op den Camp et al. 2011), a type A RR, is absent in A. hypogaea, but the expression pattern of its nonsymbiotic homolog indicates its possible recruitment during symbiosis (Figs. 4F and 5). Apart from the RRs, the cytokinin-responsive TFs like LjNIN (Tirichine et al. 2007) and MtNSP2 (Plet et al. 2011), which are central regulators for both nodule organogenesis and rhizobial infection, do not show early expression during nodulation in A. hypogaea as compared with model legumes. Thus, unlike in model legumes, the first transcriptional changes that are associated with rhizobial recognition and generation of nodule primordia in A. hypogaea may not involve NIN and NSP2. Thus, cytokinin signaling pathway and its effectors appeared to be differentially adapted for regulating nodule development in crack-entry legumes.

    Several contrasting features were also noted among the symbiotic genes involved in nodule differentiation. For example, symbiotic orthologs for MtDNF2 and MtSymCRK were not detectable in both A. hypogaea (Figs. 4G and 5) and Aeschynomene evenia transcriptome (Chaintreuil et al. 2016). As both these factors are required for suppressing defense responses during bacteroid differentiation (Berrabah et al. 2014; Bourcy et al. 2013), the divergence of DNF2 and SymCRK in crack-entry legumes further indicate that the developmental program for aeschynomenoid legumes is distinct from that of the IT-adapted legumes. LjSUNERGOS and LjVAG1 codes for topoisomerases and are required for endoreduplication of plant cells (Suzaki et al. 2014; Yoon et al. 2014). Although the orthologs of LjSUNERGOS and LjVAG1 were missing in Aeschynomene evenia transcriptome, they are expressed in A. hypogaea highlighting the differences in nodule differentiation program between the Nod-dependent (Fig. 4G) and -independent crack-entry legumes (Chaintreuil et al. 2016).

    A previous report on A. hypogaea transcriptomics (Peng et al. 2017) as well as our transcriptome highlighted the upregulation of several AP2 domain–containing ethylene-responsive TFs during nodulation. LjERF1 downregulates the expression of defense gene LjPR-10 during symbiosis (Asamizu et al. 2008), whereas, in A. hypogaea, high expression of AhERF1 is associated with high expression of a divergent group of PR-1 genes, indicating the differences in ethylene signaling between the IT-adapted legumes and A. hypogaea (Fig. 4D). Also, the expression of AhEIN2, the master regulator of ethylene signaling, was found to be constitutively high in A. hypogaea at all postinfection timepoints (Figs. 4H and 5). This predominance of ethylene signaling during A. hypogaea nodulation is consistent with the previous observations in Sesbania rostrate, in which, under submerged conditions, ethylene signaling inhibited intracellular infection via ITs while promoting intercellular infection via crack entry (D’Haeze et al. 2003; Guinel and Geil 2002).

    In summary, the comparative transcriptomic analyses have shed light on the contrasting expression profiles of DEGs in crack-entry versus IT-adapted symbiosis. Genes like VPY, FLOT2, CERBERUS, and ARPC1 that are primarily known for their role in IT elongation are significantly expressed in A. hypogaea, making it imperative to investigate their role in crack entry–adapted legumes. Similarly, the nonsymbiotic homolog of genes like SymCRK, ROP6, RR9, and SEN1 are significantly expressed during symbiosis in A. hypogaea, and it is important to know whether these genes were co-opted for symbiosis. Further investigations are also required to understand how cytokinin and ethylene signaling may have been differentially adapted in crack-entry legumes. Finally, in the absence of early induction of NIN and NSP2, it is important to identify the factors that coordinate the rhizobial invasion and nodule primordia formation in crack-entry legumes. The highlighted molecular intricacies can be further investigated for better understanding of crack-entry nodule development. But overall, our investigations have revealed a considerable overlap in expression profiles of most of the symbiotic genes between a crack-entry legume A. hypogaea and the IT-adapted model legumes, suggesting a functional conservation of most of the factors that govern the process of nitrogen-fixing symbiosis irrespective of the mode of rhizobial entry in their host plants.

    The supplementary files published online provide lists and expression patterns of i) all DEGs demonstrating major shifts in the transcriptome during symbiosis in A. hypogaea (additive to Figure 2), ii) nodulation-specific DEGs in A. hypogaea L. (Peng et al. 2017) and DEGs involved in plant-pathogen interaction and JA/SA pathway (additive to Figure 3), and iii) orthologs of symbiotic genes in A. duranensis, A. ipaensis, M. truncatula, and L. japonicus demonstrating their comparative expression profile (additive to Figures 4, 5, and 6). Orthology prediction of symbiotic genes in A. hypogaea and validation of their expression is also included.


    Plant materials and sample preparation.

    Five different developmental stages of A. hypogaea total infected roots, primordia/nodules, and uninfected roots were used in this study (UI at 0 dpi and at 1, 4, 8, 12, and 21 dpi). A. hypogaea JL24 strain seeds (ICRISAT, Hyderabad, India) were surface-sterilized with sodium hypochlorite solution (3% active chlorine) containing 0.01% Tween 20 for 15 min and were rinsed five times with sterile deionized water; the sterile seeds were soaked in sterile deionized water for germination. Germinated seeds were then transferred in pots containing sterile vermiculite and soilrite at 25°C in a growth room for 7 days before inoculation with Bradyrhizobium sp. strain SEMIA 6144 (Sinharoy et al. 2009). Bradyrhizobium sp. strain SEMIA 6144, obtained from MIRCEN, Porto Alegre, Brazil and gifted by A. Fabra, Universidad Nacional de Rio cuarto, Cordoba, Argentina. Samples were harvested, cleaned, and frozen in liquid nitrogen. Frozen samples were stored at −80°C for RNA isolation.

    Phenotypic analysis and microscopy.

    Images of whole-mount nodulated roots were captured using a Leica stereo fluorescence microscope M205FA equipped with a Leica DFC310FX digital camera (Leica Microsystems). Detached nodules were embedded in Shandon cryomatrix (Thermo scientific) and were sliced into 30-µm thick sections with a rotary cryomicrotome CM1850 (Leica Microsystems). For confocal microscopy, sample preparation was done according to Haynes and associates (2004). Sections were stained with Calcofluor (Life Technologies), propidium iodide (Life Technologies), and Syto9 (Life Technologies). Images were acquired with a Leica TCS SP5 II AOBS confocal laser scanning microscope (Leica Microsystems). For confocal and scanning electron microscopy, sample preparation was done according to (Kundu and DasGupta 2018). All digital micrographs were processed using Adobe Photoshop CS5.

    Isolation of total RNA.

    A total 100 mg of frozen plant root was ground in liquid nitrogen, and total RNA was isolated using Trizol reagent (Invitrogen). RNA degradation and contamination was detected on 1% agarose gels. RNA concentration was then measured using NanoDrop spectrophotometer (ThermoScientific). Additionally, RNA integrity was assessed using the Bioanalyzer 2100 system (Agilent Technologies). Finally, the samples with RNA integrity values above 8 were used for library construction.

    Library construction and sequencing.

    Eighteen RNA libraries were prepared, using an IlluminaTruSeq stranded mRNA sample preparation kit, by the MGX-Montpellier GenomiX core facility, Montpellier, France (MGX). The protocol first required the selection of polyadenylated RNAs on oligodT magnetic beads. Selected RNAs were chemically fragmented and the first-strand cDNA was synthesized in the presence of actinomycin D. The second-strand cDNA synthesis incorporated dUTP in place of dTTP in order to quench the second strand during amplification. A 3′ ends adenylation was used to prevent fragments from ligating to one another during the adapter ligation process. The quantitative and qualitative validation of the library was performed by qPCR, Roche Light Cycler 480 and cluster generation, and primary hybridization are performed in the cBot with an Illumina cluster generation kit. The sample libraries were sequenced on an IlluminaHiSeq 2000, sequencing by synthesis technique performed by MGX, and 50-bp single-end reads for each library were generated (Fuller 1995).

    Illumina reads mapping and assembly.

    Quality control and assessment of raw Illumina reads in FASTQ format were done by FastQC software (version 0.11.5) to obtain per-base quality, GC content, and sequence-length distribution. Clean reads were obtained by removing the low-quality reads, adapters, and poly-N–containing reads, using Trimmomatic v0.36 software (Bolger et al. 2014). Clean Reads were simultaneously aligned to the two wild peanut diploid ancestors A. duranensis (AA) and A. ipanensis (BB) reference genomes, using TopHat2 version 2.0.13, which is a fast splice junction mapper for RNA-Seq reads, for which up to three mismatches per readwere authorized (Bertioli et al. 2016; Trapnell et al. 2010). It aligns RNA-Seq reads, using the ultrahigh-throughput short read aligner Bowtie2 version 2.2.3, and then analyzes the mapping results to identify splice junctions between exons (Langmead et al. 2009). Alignment files were combined and analyzed into Trinity for genome-guided assembly (Grabherr et al. 2011).

    The reference-based assembly was compared with its respective transcript files from annotated reference genomes, using BLAT (Kent 2002). An e-value cutoff of 1e−05 was used to determine a hit. The annotated hits were further analyzed in this study. Genome annotation files in generic feature format (GFF) were downloaded from PeanutBase (Dash et al. 2016). Estimation of the gene expression level of each annotated transcript was performed by StringTie v1.3.3, which takes a sorted sequence alignment map or binary file for each sample along with genome annotation files (Pertea et al. 2015). Counting reads in unannotated genes is performed with HTseq-count after converting the coordinates of these genes (chromosomal positions) to a GFF file. As the data were from a strand-specific assay, the read was mapped to the opposite strand of the gene. The resulting gene transfer format, normalized gene locus expression level as FPKM, transcripts per million, and count files for each sample were further analyzed for fold change analysis in gene expression levels. RNA-seq read alignments have been performed using the Burrow-Wheeler Aligner (version 0.7.15-r1140) (Li and Durbin 2009) against MtRPG, MtDNF2, and singleton sequences..

    Identification of DEGs and functional GO and KEGG pathway analyses of the DEGs.

    Before statistical analysis, genes with less than two values lower than one count per million were filtered out. EdgeR 3.6.7 package was used to identify the DEGs (Robinson et al. 2010). Data were normalized using the trimmed mean of M values method. Genes with an adjusted P value less than 5% (according to the FDR method, using Benjamini-Hochberg correction) and log2-fold change > 1 was called differentially expressed. Venn diagrams were generated using InteractiVenn (Heberle et al. 2015) and the hierarchical heatmap was generated usingTM4MeV (Howe et al. 2011) and the values from the Venn diagram. PCA was performed by web-tool ClustVis (Metsalu and Vilo 2015). Principal components of the variability of expression data were calculated using the PCA method R-Package. Detailed functional annotation and explanations of DEGs were extracted from the GO database (Ashburner et al. 2000) and GO functional classification analysis was done using software WEGO (Ye et al. 2006). The GO terms for DEGs in genome annotation were also retrieved from the GFF file downloaded from the PeanutBase website. To identify important and enriched pathways involved by the DEGs, the DEGs were assigned to KEGG pathways, using the GenomeNet database KAAS job request webserver (Kanehisa and Goto 2000) against A. duranensis and A. ipaensis gene datasets. Enriched KO and GO terms were obtained by a developed Python script that uses hypergeometric test and Bonferroni corrected P value < 0.05.

    Identification of symbiotic orthologous gene in A. hypogaea.

    Candidate symbiotic genes were identified in A. duranensis, A. ipaensis, L. japonicus, and M. truncatula using BLASTN searches with reported nucleotide sequence of genes from L. japonicus and M. truncatula. The homologous genes were searched in A. duranensis and A. ipaensis at PeanutBase, the M. truncatula Mt4.0v1 genome was searched at MtGEA, and the L. japonicus v3.0 genome was searched at LjGEA. Initial searches were conducted with E value = e−5. The results were manually validated for the presence of an orthologous gene in an open reading frame and were searched for orthologs using BLASTP. Orthology of the genes were validated by generating maximum likelihood (SymCRK, RR9, ROP6, SEN1, FLOT4, NOOT, and DNF2) and a neighbor-joining phylogenetic tree, using amino acid sequences and nucleotide sequences (ENOD40) in MEGA 6.0, using 1,000 bootstrap replications (Tamura et al. 2013). Syntenic analysis of RPG and DNF2 between the legumes was performed using the Legume Information System genome context viewer.

    qRT-PCR validation.

    Total RNA (500 ng) was reverse-transcribed by using Super-ScriptIII RT (Life Technologies) and oligo(dT). RNA quantity from each sample in each biological replicate was standardized prior to first-strand cDNA synthesis. qRT-PCR was performed using Power SYBR green PCR master mix (Applied Biosystems), with primers as designed using software Oligoanalyser (IntergratedDNA Technology). Calculations were done using the ΔΔ cycle threshold method, using AhActin as the endogenous control. The reactions were run in the Applied Biosystems 7500 Fast HT platform using the following protocol: 1 cycle at 50°C for 2 min, 1 cycle at 95°C for 5 min, 40 cycles at 95°C for 30 s, 54°C for 30 s, 72°C for 30 s, followed by melt curve analysis for 1 cycle at 95°C for 1 min, 55°C for 30 s, and 95°C for 30 s. A negative control without cDNA template was checked for each primer combination that was designed using OligoAnalyzer 3.1. Results were expressed as means ± standard error of the number of experiments.

    Data availability.

    The raw FASTQ files for the 18 libraries were deposited in the Gene Expression Omnibus database of the National Center for Biotechnology Information under accession number GSE98997.



    The raw FASTQ files for the 18 libraries were deposited in the Gene Expression Omnibus database of the National Center for Biotechnology Information under accession number GSE98997.

    Kanchan Karmakar and Anindya Kundu contributed equally to this work.

    Funding: This work was funded by grants from the government of India: IFCPAR/CEFIPRA (IFC/5103-4/2014/543); DBT-CEIB (Centre of Excellence and Innovation in Biotechnology, BT/01/CEIB/09/VI/10); DBT-IPLS (BT/PR14552/INF/22/123/2010); DST (EMR/2015/001006); fellowship to Kanchan Karmakar and Ahsan Z Rizvi (IFCPAR/CEFIPRA: IFC/5103-4/2014/543); fellowship to Anindya Kundu (Council of Scientific and Industrial Research, CSIR-09/028[0756]/2009–EMR–I).

    Current address for Ahsan Z Rizvi: 112, INSERM U981, Bâtiment Médecine Moléculaire (B2M), Gustave Roussy, 114, rue Edouard Vaillant, 94805 Villejuif Cedex-France