An RNA-Seq based gene expression atlas of the common bean

Phaseolus vulgaris gene expression atlas (PvGEA)

To facilitate additional use of the RNA-Seq data generated for these analyses, we have developed a web-accessible database, The Pv Gene Expression Atlas (Pv GEA), available at http://plantgrn.noble.org/PvGEA/. This database was built using a similar database structure, web application, architecture and tools as the LegumeIP platform [14] to retrieve and visualize the gene expression patterns using RNA-seq data. To facilitate the mining of the data included in Pv GEA, we have provided the capability to: (i) visualize expression profiles of genes of interest, (ii) identify genes exhibiting certain expression patterns in specific tissues, (iii) identify genes and gene expression patterns based on http://www.phytozome.net annotation terms; and (iv) download the entire data set, either raw or normalized, in tabular form to facilitate the analysis of more complicated biological questions. Using the predicted gene calls of the G 19833 Pv genome to build the Pv GEA database means it can be easily expanded to integrate RNA-Seq data from future experiments. Currently, Pv GEA includes gene expression profiles from 24 samples isolated from roots, root nodules, stems, leaves, flowers, pods, and seeds at various developmental stages under ideal growth conditions. Included in this dataset are transcripts from eight samples including nodule, root, and leaf tissues for plants having either fix + or fix- root nodules; providing preliminary data on the impact of nodulation and N fixation on gene expression, an important biological process for legumes.

To ensure our transcriptome analysis was reliable we visualized the expression profile of the purine and ureide biosynthesis pathways (Additional file 5) and conducted qPCR on 85 additional genes (for details see experimental procedures and Additional file 6). In warm season legumes, such as soybean and Pv, purine biosynthesis is known to be highly up regulated in nodules [15]. Our data is consistent with this as the genes in this pathway are highly up regulated in nitrogen fixing nodules compared to all other tissues (Additional file 5). The enzyme uricase (Additional file 3c) degrades ureate to allantoin, which is the main supply of N for plant nutrition. In our Pv data allantoinase, the enzyme responsible for allantoin degradation, is highly expressed early in seed and pod development, likely providing N to developing seeds (Additional file 3d). Expression of uricase and allantoinase in aerial tissues suggests ureides are degraded after being transported from the nodules. Leaves, seeds, and pods can then utilize the released NH₃ and CO₂ in a variety of cellular processes. These results are consistent with reports of high ureide levels observed in developing Pv seeds [16] and high allantoinase enzyme activity throughout pod development and seed filling as measured by Thomas et al. [17]. Additionally, Pellissier et al. [18] reported allantoin transporters highly expressed in developing pods and seed coats of Pv. The results of these three studies combined with the gene expression patterns observed in our data highlight the importance of ureide metabolism in aerial tissues to provide N for developing tissues.

Gene expression analysis

Genes exhibiting tissue specific expression (Additional file 2 and Additional file 3a) are involved in a variety of gene ontology biological processes. Genes uniquely expressed in leaves are involved in amino acid phosphorylation, DNA and protein binding. Genes uniquely expressed in seed tissues include processes such as carbohydrate metabolism, as exemplified by a starch branching enzyme (SBE), which is important in amylopectin synthesis, a carbohydrate precursor [19]. In our data, Pv SBE expression (Phvul.005G040300.1) is highest in developing seeds (Additional file 3e). These results are consistent with the high carbohydrate composition of Pv seeds reported by Broughton et al. [1]. Nodule specific transcripts annotate as involved in oxidoreductase activity, amino acid phosphorylation and membrane transport/signaling are highly and uniquely expressed; reflecting the importance of nutrient transport to the nodule and the high energy cost of N fixation. In root tissues, genes involved in pectinecterase, carbohydrate metabolism, iron ion binding, oligopeptide membrane transport, and lipid metabolisms are uniquely expressed. Expression of these genes illustrates the role of the root in nutrient acquisition and the importance of root growth for plant health.

The 6,667 known transcription factors (TF) in soybean (downloaded from SoyDB) [20] were compared to Pv genes using TBLASTN (e-value 1e-30). This analysis identified 3,726 putative TFs in the Pv genome, representing 52 of the 64 transcription factor families in soybean (Additional file 7). The 3,726 TFs identified in Pv is almost exactly half the number identified in soybean, as expected since soybean has undergone a whole genome duplication event not experienced by Pv[12, 13]. The average expression of TFs in all seed tissues was much lower than that of other tissues, including developing pods. Fisher tests confirmed 26 TF families exhibited statistically enhanced or repressed expression by tissue type (Figure 1). Twenty-one TF families exhibited altered expression patterns in a single tissue. Five other TF families (AS2, AUX, NAC, SBP, and ZD-HD) showed expression patterns that were statistically significant in multiple tissues.

Seed development and metabolism

Pv seeds are an integral component of diets around the world. Unlike soybean seeds, which are valued for high oil and protein content, Pv beans provide high levels of protein and carbohydrates making it a highly nutritious food for human consumption. The seed and pod samples represent an extended time-course collection of the same tissue spanning great developmental changes. Approximately half (16,292) of the 31,638 predicted genes are expressed in Pv seeds with 12,182 genes expressed in all three stages of seed development examined. In Pv pods 17,248 genes are expressed in at least one of the four developmental stages examined. We identified genes differentially expressed between seeds and pods at the three stages of development (8,189 genes; Additional file 8, Additional file 9, and Additional file 10), and genes differentially expressed within seeds (2,694) and pods (13,125) throughout development (Additional file 11 and Additional file 12). Additionally, we identified 9,702 genes with consistently decreasing expression levels as the seed develops, including 1,196 TFs. By comparison, 753 genes were identified with increasing expression levels as the seed develops (Figure 2) including 70 TFs from 25 families, including HB (Phvul.007G064100.1), MYB (Phvul.001G025200.1), and ARF (Phvul.011G080100.1) (Additional file 3f). Members of three of these TF families (HB, MYB, and ARF) are seed specific in Arabidopsis, M. truncatula, and our Pv data [21, 22]. Additionally, members of these TF families are among the TFs identified by Le et al.[23] as differentially expressed between structures of developing Phaseolus coccineus (scarlet runner bean) seeds that are also differentially expressed in Pv seeds at different developmental stages. All these results are consistent with a study by Hajduch et al. [24], which determined the expression of proteins involved in primary and secondary metabolism, cell growth and division, signal transduction, and protein synthesis all decrease as the seeds develop. Conversely, in the pods we identified 39 genes with expression levels steadily decreasing and 1,236 genes with consistently increasing expression patterns (Figure 2).

Soybean and Pv, although closely related species, have distinct seed compositions. While soybean accumulates oil and protein, Pv accumulates carbohydrates and protein [1, 25]. Comparing the 500 most highly expressed genes in both Pv and soybean seeds [26] allowed us to identify genes important for general seed development (Figure 3a). Genes involved in carbohydrate biosynthesis are highly up regulated in Pv seeds while genes involved in fatty acid biosynthesis are highly up regulated in developing soybean seeds (Figure 3b and c). The starch synthase (STS) genes, particularly Phvul.001G082500.1, are highly expressed in Pv seeds (RPKMs: S1 = 240, S2 = 286) but not in soybean seeds (RPKM < 15) (Figure 3b and c). STS is required for amylopectin biosynthesis, a component of carbohydrates. Similarly, sucrose acts as a key regulator of seed carbon flux; high sucrose synthase (SS) activity in developing seeds may channel available carbon towards carbohydrate biosynthesis and away from fatty acid biosynthesis [27]. In our Pv RNA-Seq data, genes encoding SS are highly expressed in developing seeds (Figure 3b). In Arabidopsis, SS loss of function mutants favor fatty acid and protein biosynthesis over starch biosynthesis in seed development, resulting in a 55% increase of fatty acids and a near 70% reduction in starch content of mature seeds [28]. The role of SS in regulating carbohydrate synthesis is further supported by the low expression in developing soybean seeds (Figure 3c) (Severin et al. [26]), which are valued for oil and protein. The synthesis of polyunsaturated fatty acids is regulated by FATTY ACID DESATURASE 2 (FAD2) [29]. In soybeans, FAD2 is highly expressed in developing seeds while in Pv FAD2 is expressed early in seed and pod development, but at a much lower level (Figure 3b and c).

Seed development in multiple species is regulated by four master TFs: LEAFY COTYLEDON1 (LEC1), LEAFY COTYLEDON2 (LEC2), ABSCISIC ACID INSENSITIVE3 (ABI3), and WRINKLED 1 (WRI1)[21, 22, 30]. Using BLASTP, we queried the Arabidopsis protein sequences to identify homologous sequences in the Pv predicted genes (Figure 4). The homolog for LEC2 was only weakly expressed (RPKM = 4) mid-seed development in Pv. Seeds of Arabidopsis loss of function lec2 mutants accumulated 15% less protein and 30% less oil while the seed starch content increased five fold as compared to wild type plants [31]. The altered seed composition of lec2 mutant plants closely resembles that of Pv, suggesting down regulation of LEC2 may affect seed composition. LEC2 controls the gene expression of WRI1, which also exhibits low expression patterns in Pv developing seeds (RPKM: SH = 5, S1 = 9, S2 = 5). WRI1 expression modulates the expression of a set of genes controlling late glycolysis and fatty acid biosynthesis. The low expression of both LEC2 and WRI1 may relate to the lower oil composition of Pv.

Abscisic acid (ABA) is a key hormone in seed development, important in developing desiccation tolerance and entrance into dormancy [30]. ABA accumulation in seeds is both temporally and spatially regulated [27]. We found high expression of ABA biosynthesis genes (Phvul.002G018700.1 and Phvul.005G031500.1) in developing seeds, with expression decreasing as the seeds matured (Additional file 3g). ABA biosynthesis is regulated by 9-cis-EPOXYCAROTENOID DIOXYGENASE (NCED9) [32]. The expression pattern of NCED9 (Additional file 3 h) and ABA biosynthesis genes (Additional file 3 g) in Pv developing seeds is consistent with those from developing seeds in Arabidopsis [33–35].

Trehalose biosynthesis is important in regulating both seed composition and nodule metabolism [36, 37]. In seeds, TREHALOSE 6 PHOSPHATE SYNTHASE 1 (TPS1), the enzyme responsible for converting glucose-6-phosphate to trehalose-6-phosphate is thought to regulate sugar utilization [32]. In Arabidopsis thaliana tps1 null mutants, both sucrose and starch content of seeds dramatically increased [38]. In Pv, TPS1 expression drops dramatically as the seed develops (Additional file 3i), corresponding with increased SS expression. We hypothesize that, as in Arabidopsis, the reduced TPS1 expression promotes increased carbohydrate biosynthesis in Pv seeds.

Nodule analysis

Legumes have established a unique symbiotic relationship with Rhizobium, which allows legumes to fix atmospheric N₂ into biologically useful NH₃. For this experiment plants were provided with nutrients containing NO₃⁻ nitrogen for optimal growth conditions or inoculated with either effective fix + Rhizobium tropici CIAT899 or ineffective fix- Rhizobium giardini 6917 to induce nodulation. Plants inoculated with normal (fix+) R. tropici appeared green and healthy, though smaller than plants provided with nitrate fertilizer (Additional file 13a). This phenotype is consistent with previous studies reporting the overall growth of N₂ fixing Pv plants is restricted compared to fertilized plants, likely due to altered carbon partitioning [39, 40]. Fix + plants inoculated with R. giardini were nitrogen (N) deficient, exhibiting severe chlorosis and a stunted phenotype (Additional file 13a). Small, pre-fixing white nodules (N5) were isolated from root tissues of plants inoculated with effective fix + R. tropici five days after inoculation (DAI). At 21 DAI nodules were collected from plants inoculated with either fix + R. tropici (NE) or fix- R. giardini (NI) (Additional file 13b and c respectively). Microscopic imaging of fix-nodules 21 DAI (NI) revealed early senescing cells with few, if any, infected cells compared to fix + nodules formed 21 DAI (Additional file 13d and e). In situ hybridization analysis was used to visualize the localization pattern of leghemoglobin transcripts in these two nodule types. Fix- nodules collected 21 DAI (NI) exhibited little to no expression of leghemoglobin transcripts while fix + nodules collected 21 DAI (NE) exhibited high expression levels, likely mirroring the bacteroid colonization patterning (Additional file 13f and g) and directly reflecting the gene expression patterns observed in the RNA-Seq data (Figure 5). Nodule acetylene reduction assays failed to detect nitrogenase activity at 5 DAI in pre-fixing nodules and at 21 DAI with fix-nodules. Fix + nodules (and associated roots) from plants 21 DAI reduced 320 nm/hr/gfw (roots) of acetylene, indicating high N₂ fixation activity. Leaf tissue was also collected at 5 DAI for fix + plants and 21 DAI for both fix + and fix- inoculated plants (see Table 1).

Comparing 5 DAI pre-fixing nodules (N5) and 21 DAI fix + nodules (NE) revealed 2,932 differentially expressed genes (Additional file 14). Comparing 21 DAI fix + and fix- nodules (NE vs NI) identified 2,953 differentially expressed genes (Additional file 15). Additionally, we found 245 nodule specific genes; genes expressed in any and/or all the nodule tissues sampled, but not expressed in any other tissues (Additional file 4). Comparing these nodule specific genes to those identified in two soybean gene atlases [26, 41] identified 21 nodule specific homologs common to both species (Additional file 16) including seven TFs and four transporters. The conserved expression of these genes highlights the importance of regulating gene expression, but also the exchange of nutrients between nodules and the plant roots. Five of these sequences have no known annotation, though nodule specific expression in both species suggests these are important candidates for characterization in future nodulation and nitrogen fixation research.

Cognate genes involved in nodule development and the establishment of N fixation have been identified and cloned from multiple species [42, 43]. Using BLASTN, the homologous sequences in common bean were identified and gene expression patterns were visualized as a heat-map (Figure 5). Upon further analysis of these nodule cognate genes, we detected three expression profiles: those highly up regulated early in nodule development (N5), those highly up regulated in 21DAI fix + nodules (NE), and those highly up regulated in 21 DAI fix- nodules (NI).

The autoregulation of nodulation (AON) pathway mediates nodule formation [44]. ASTRAY and UFD1a proteins, both expressed in leaves, function in the AON pathway [44]. ASTRAY, encodes a bZIP TF that interacts with a nodulation autoregulation receptor kinase (NARK) [44, 45]. UFD1a expression indicates the presence of Q, a root derived signal induced upon compatible rhizobial infection [44]. In soybean, three candidates for Q have been identified, all of which are CLE peptides [46]. Pv encodes a single CLE homolog (Phvul.005G097000.1), which is highly up regulated in N5 and NE, but noticeably absent in NI (Figure 5). Surprisingly, we observed aerial AON genes (ASTRAY and UFD1a) expressed higher in leaves of plants inoculated with fix- rhizobia than in leaves of plants inoculated with fix + rhizobia (Figure 5). We hypothesize the fix- inoculated plant may up regulate the AON pathway to minimize resources allocated to nodules as part of a survival strategy.

Early in nodule development (N5) nod factor receptors (NFR) and nodulation signaling pathway (NSP) TFs are highly expressed (Figure 5). The early calcium spiking response induces both a calmodulin dependent protein kinase (CCaMK/DMI3), which is required and sufficient for nodule organogenesis [47–49], and nuclear porin proteins (NUPs) [44] (Figure 5). Additionally, genes involved in infection thread formation and elongation including ERN1, FLOT, VPY, PUB1 and RPG[50–58] are highly up regulated in N5 (Figure 5). VPY and PUB1, which are up regulated in N5 and NE, are involved in rhizobial recognition, attachment, entry, and initiation of the infection thread [53, 58]. Nodule organogenesis involves the altered differentiation and division of root cortical cells prior to the formation of the nodule primordia. In Medicago truncatula, these processes are dependent on ENOD40[59], which is highly expressed in N5 (Figure 5). HAP2, which promotes nodule development and the release of bacteria from the infection thread [60] is expressed highest in N5, but remains elevated in NE (Figure 5).

Among the nodule cognate genes, those most highly expressed in NE are: NIN transcription factors, CYCLOPS, and CERBERUS with expression profiles increasing from 5 DAI to 21 DAI in fix + nodules (2X, 2.4X, and 4X respectively, Figure 5). CYCLOPS expression is required for rhizobia infection. Nodule organogenesis is dependent on CERBERUS gene expression [49, 67–69].

The primary function of nodules is to fix N₂ to NH₃. Nitrogenase, the enzyme responsible for nitrogen fixation, requires sulfur and a near anaerobic environment to function. The gene encoding SYMBIOSOME SULFATE TRANSPORTER 1 (SST1), which transports sulfur into bacteroids [70], is expressed 15 fold higher in NE than in NI (Figure 5). The genes encoding leghemoglobin, which sequesters oxygen [71], are only expressed in NE (Figure 5). Once N₂ is reduced to usable ammonia it must be assimilated for use and distribution throughout the plant. NADH-dependent glutamate synthase (NADH-GOGAT) is a key enzyme in ammonia assimilation [72]. Two NADH-GOGAT genes (Phvul.009G053900.1 and Phvul.001G076400.1) are expressed 5 and 10-fold higher in NE than in NI (Figure 5), reflecting the difference in effectiveness. Transcripts encoding glutamine synthase, uricase, and allantoinase (Additional file 3 j, c, and d), each involved in primary ammonia assimilation, exhibit similar expression patterns. Nod41 was identified in Pv by Olivares et al. [73] as important in controlling defense responses during symbiotic interactions and maintaining the integrity of the uninfected root nodule cells. Our data is consistent with this hypothesis as Nod41 is expressed 7-fold higher in 21 DAI fix + nodules than in 21 DAI fix- nodules (Figure 5).

Genes up regulated in NI include those involved in the GO processes of autophagy and early senescence including ubiquitination, proteolysis, peptidyolysis and apoptosis. Respiratory burst oxidase homolog (RBOH) genes, which generate reactive oxygen species (ROS) [74], are up regulated four-fold in NI. Increased ROS production is a common defense response to pathogen attack (ie: ineffective rhizobia) and in response to abiotic stress, including nitrogen deficiency. Additionally, we observed high expression of leucine rich repeat (LRR) genes in NI (Additional file 19), likely reflecting a defense response as the plant reacts to invading bacteria [75]. Also highly expressed in NI are genes involved in oxidation-reduction processes, membrane transport, protein binding, and amino acid phosphorylation. Among the nodule cognate genes are those encoding a second group of NIN and NSP TFs (Figure 5). This result suggests a group of alternative TFs may be induced in NI versus NE.

Type A-response regulators (RRs) negatively regulate cytokinin signaling [76]. In lotus, RRs are rapidly induced following rhizobial inoculation in root hairs and dividing cortical cells [77], repressing the cytokinin signaling pathway [48]. Inhibition of the cytokinin-signaling pathway may contribute to plant and bacterial cell differentiation. In Pv, the gene encoding RR5 is more highly expressed in NI than in NE (Figure 5). Genes encoding both CASTOR and POLLUX, both of which are required for the activation of voltage gated calcium (Ca²⁺) channels [74], are highly expressed in NI (Figure 5). The high expression of CASTOR and POLLUX genes in NI may suggest that at 21 DAI the plant is still attempting to establish SNF or may reflect the induction of Ca²⁺ channels by ROS as described in Pisum sativum[78, 79].

Roots and nitrogen

Gene expression profiles of Pv roots were examined from plants grown under three conditions 1) those from NO₃⁻ fertilized plants (RF), 2) those from plants with fix + effective nodules (RE), and 3) those derived from plants having fix- ineffective nodules (RI). RF and RE roots had adequate N for growth, while RI roots were N deficient. We identified 1,714 genes differentially expressed between the three 21 DAI root samples (RI, RE, and RF) (Additional file 20). The majority of these genes (1,668 genes) are differentially expressed between fertilized roots (RF) and nodulated roots (either RE or RI).

Comparing gene expression patterns between RI (roots from fix- plants) and both RF and RE identified 426 and 46 genes differentially expressed between root samples respectively (Figure 6c). Additionally, 210 genes were differentially expressed between RE and RF. These 210 sequences represent genes differentially expressed due to the N source (Figure 6c). A similar comparison of leaf tissues collected from each of the plants revealed 116 genes differentially expressed between + N leaves due to different N sources (Figure 6a). These 116 genes indicate that the source of N (either via N₂ fixation or NO₃⁻ fertilization) has a long-term impact on plant gene expression. Among all 2,641 genes differentially expressed between samples due to the N source are 340 TFs, the majority of which are up regulated in –N tissues (Additional file 21 a and b).

In the presence of abundant NO₃⁻, plants will preferentially take up and utilize NO₃⁻ rather than develop SNF. NO₃⁻ transporters exhibit either low (NRT1) or high (NRT2 and NRT3) NO₃⁻ affinity [80–84]. Examination of the expression patterns of NO₃⁻ transporters in our root samples revealed plants provided with NO₃⁻ as a fertilizer induce NRT1 and NRT3 gene expression, reflecting the abundance of available NO₃⁻ (Figure 6d). In –N roots, only members of the high affinity NRT2 gene family are up regulated (Figure 6d). N-deficient plants may up regulate N transporters in an attempt to increase the N content of the plant to mitigate –N stress. Members of the low affinity NRT1 gene family are also up regulated in the roots of N₂ fixing plants (Figure 6d). The constitutive expression of NRT1.1 is consistent with the recent evidence suggesting that it serves as both an N sensor and transporter [85–87]. Expression of NRT1.1 in fix + plants may be involved in N sensing.

Once N is within the plant it must be assimilated. Glutamine synthetase (GS) functions as a primary enzyme for NH₄ assimilation produced from N₂ fixation or NO₃⁻ nutrition [88, 89]; synthesizing glutamine from NH₃ and glutamate (Additional file 3j). In SNF plants, the majority of glutamine is committed to the de-novo purine biosynthesis pathway. Alternatively, glutamine may be reduced by GOGAT. Consistent with previous studies [72], NADH-GOGAT expression is highest in roots (particularly YR) and NE while Fd-GOGAT is expressed highest in leaf tissues (Additional file 3 k). Plants provided with NO₃⁻ fertilization utilize glutamine in the synthesis of asparagine via asparagine synthetase (AS) [88], which is most highly expressed in fertilized root (RF) tissues (Additional file 3 l).

Our data shows increased expression of auxin response factors unique to –N leaves (Figure 6b). This gene expression pattern indicates increased auxin levels in –N leaves, supporting auxin as the N signal. The availability of N for proper growth and development is likely monitored throughout the plant. Auxin has been proposed as an N status mediator, signaling from root to shoot [80]. Under low N and other nutrient stress conditions in the shoot, increased auxin is transported to the roots to enhance lateral root development, a hallmark response of –N plants [46, 80, 87].

Plant materials and growth conditions

Phaseolus vulgaris cv. Negro jamapa seeds were grown as described by O’Rourke et al. [90]. At the emergence of the unifoliate, pots were assigned to one of three nitrogen (N) treatments; inoculated with Rhizobium tropici CIAT899 (fix+), Rhizobium giardini 6917 (fix-), or fertilized with a full nutrient solution. Pots assigned to the fertilization treatment were watered daily with a nutrient solution of 9 mM KNO₃, 2.5 mM Ca(NO₃)₂, 1.0 mM Ca(H₂PO₄)₂, 1.0 mM MgSO₄, 12 μM Fe (as FeEDTA), 4.0 μM MnCl₂, 22.0 μM H₃BO₃, 0.4 μM NaMoO₄, and 1.6 μM CuSO₄. Twenty-four tissue samples were collected throughout development and across all N treatments (for details see Table 1). For each nitrogen treatment, a representative plant was chosen and leaf, root, and nodule tissue samples were collected. Two plants were maintained in full nutrient solution fertilized pots. From these plants, root, leaf, flower, stem, seed and pod tissues were collected (see Table 1 for details). All tissue collected for RNA -Seq analysis was immediately flash frozen in liquid nitrogen.

RNA extraction and expression analysis

Total RNA was purified from 24 tissue samples using RNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA). For RNA-Seq analysis, RNA samples were shipped on dry ice overnight to National Center for Genome Resources (NCGR, Santa Fe, NM) for sequencing as described by Severin et al. [26]. Illumina reads generated from all 24 samples are available at the NCBI SRA browser, accession number SRP046307. Illumina reads passing quality control standards, approximately 25 million sequences per sample or 596 million 36 bp sequences total, were mapped to the Phaseolus vulgaris v1.0 genome available at http://www.phytozome.net using the program Bowtie [91]. Reads were also mapped to the predicted transcripts to account for splicing events. Reads that mapped to more than one location were counted at each mapping location. Of the 596 million reads generated, 406 million (89%) mapped to the genome with 14% of those mapping to non-coding regions. Raw gene expression counts were normalized using the RPKM (reads/Kb/Million) method [92, 93] using custom R scripts. To ensure expression profiles were not statistical artifacts as described by earlier studies [94, 95], we determined an RPKM of 3 represents a 2X coverage across the coding region, assuming equal distribution, and would be the minimum level at which a gene would be considered expressed; genes with an RPKM < 3 were considered silent. Transcripts differentially expressed between libraries were identified using NOIseq [96]. Differentially expressed transcripts were required to have > 2-fold change in expression between samples and a probability of differential expression > 0.9. Additionally, one of the two sequences was required to have an RPKM > 3. Heatmaps illustrating expression patterns of various subgroups of transcripts were generated in R as described by Severin et al.[26].

To identify genes exhibiting enhanced expression in NE, we determined the Euclidian distance between Z-scores for each gene. A threshold of two Euclidian distances was established as significant; genes within the threshold were identified as co-expressed.

Real time quantitative RT-PCR (qPCR) and housekeeping genes

RNA was extracted using the RNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA) from three biological replicates of tissues grown in growth chambers under the same conditions described above. Transcript specific primers were designed using Primer3 (Frodo.wi.mit.edu). The qPCR analysis was run as described by O’Rourke et al. [90] for 85 genes identified as differentially expressed by NOIseq (Additional file 6). 92% of the qPCR experiments confirmed the differential expression measured by NOIseq analysis in at least two of the three biological replicates.

Genes exhibiting stable expression profiles between tissues and across growth conditions were identified as described by Severin et al. [26] (Additional file 22). The 10% of transcripts with the lowest CV were selected as potential housekeeping genes. This suite of stably expressed transcripts may be useful in future experiments for normalizing gene expression patterns across a variety of experimental conditions, or tissues [97]. One of these genes (Phvul.006G165300.1) was successfully utilized as a housekeeping sequence in the qPCR analysis. Comparing the housekeeping genes proposed in this study to the eleven potential housekeeping genes identified for Pv under biotic and abiotic stress by Borges et al. [98] found seven genes common to both lists, illustrating the utility of this list for multiple experimental conditions.

Acetylene reductase assay

Acetylene reduction assays were performed as described by Vance et al. [99] with the following modifications. Plant roots from six biological replicates of each sample (roots inoculated with fix + R. tropici CIAT899 at 5 and 21 DAI and roots inoculated with fix – R. giardini 6917 21 DAI) were placed in 500 ml airtight glass containers equipped with serum stoppers. 50 ml of air was removed from each container and replaced with 50 ml of ethylene; samples were incubated at room temperature for one hour, at which time 10 ml of gas was withdrawn from the container for analyses as previously described by Vance et al. [99].

Nodule In Situs

A partial coding sequence for Pv leghemoglobin (645 bp) was PCR amplified and cloned into pBSSK+. Nodules inoculated with fix + R. tropici CIAT899 and fix – R. giardini 6917 were collected 21 DAI and analyzed as described by Sbabou et al. [100].

Availability of supporting data

The expression data used in this study is publicly available at the NCBI short read archive; accession SRP046307. Additionally, the raw and normalized datasets can be downloaded and explored at the Phaseolus vulgaris Gene Expression Atlas (Pv GEA) website, http://plantgrn.noble.org/PvGEA/.