Transcriptome Characterization of Developing Bean (Phaseolus vulgaris L.) Pods from Two Genotypes with Contrasting Seed Zinc Concentrations.

Dry bean (Phaseolus vulgaris L.) seeds are a rich source of dietary zinc, especially for people consuming plant-based diets. Within P. vulgaris there is at least two-fold variation in seed Zn concentration. Genetic studies have revealed seed Zn differences to be controlled by a single gene in two closely related navy bean genotypes, Albion and Voyager. In this study, these two genotypes were grown under controlled fertilization conditions and the Zn concentration of various plant parts was determined. The two genotypes had similar levels of Zn in their leaves and pods but Voyager had 52% more Zn in its seeds than Albion. RNA was sequenced from developing pods of both genotypes. Transcriptome analysis of these genotypes identified 27,198 genes in the developing bean pods, representing 86% of the genes in the P. vulgaris genome (v 1.0 DOE-JGI and USDA-NIFA). Expression was detected in 18,438 genes. A relatively small number of genes (381) were differentially expressed between Albion and Voyager. Differentially expressed genes included three genes potentially involved in Zn transport, including zinc-regulated transporter, iron regulated transporter like (ZIP), zinc-induced facilitator (ZIF) and heavy metal associated (HMA) family genes. In addition 12,118 SNPs were identified between the two genotypes. Of the gene families related to Zn and/or Fe transport, eleven genes were found to contain SNPs between Albion and Voyager.

Background

Zinc is essential to human health and nutrition. Zinc is an important enzyme cofactor and component of proteins, and is needed for DNA synthesis, RNA transcription, and cell division [1]. Human Zn deficiency symptoms are quite varied, including reduced immune function, fetal brain cell development, and reproductive and cognitive development [2]. Mild to moderate Zn deficiency is common, especially in populations consuming vegetarian diets rich in unrefined cereals [3]. Biofortification of staple foods such as wheat and dry beans with Zn is one agricultural science based approach being developed and applied to combat micronutrient malnutrition [4].

Dry beans (Phaseolus vulgaris L.) are a nutrient dense food crop and a dietary staple in East Africa and Latin America. Genotypic variability for seed Zn levels is relatively high within the species and Zn seed levels from 20 to 59 μg g-1 have been observed [5,6]. Understanding the genetic control of seed Zn content has the potential to improve the breeding process for this important nutritional trait by identifying candidate genes for marker assisted selection and also increase the overall Zn content levels achievable through breeding.

Numerous genes involved in Zn transport have been characterized in model plant species including Arabidopsis and Medicago [7]. Major gene families shown to play a role in transport of Zn include, ZIP, YSL, NAS, ZIF, HMA and MTP. The ZIP family is made up of ZRT (zinc related transporter) and IRT (iron related transporter) like proteins. The common feature of members of this family is eight transmembrane domains and a metal binding domain [8]. In addition, transcription factors that regulate ZIP genes include members of the basic region leucine zipper (bZIP) gene family. bZIP19 and bZIP23 genes have been shown to interact with ZIP genes in Arabidopsis [9]. Yellow stripe like (YSL) is a gene family that transport metal-NA complexes long distance. In Arabidopsis AtYLS2 is responsible for mobilization of micronutrients such as Mn, Zn, Cu and Fe from leaves and for loading of Fe-NA complex into seed [10]. Nicotianamine (NA) a non proteinogenic amino acid chelates Fe and Zn phloem movement to sink tissue [11]. Four NA synthase genes have been characterized [12] and are related in reproduction and seed Fe loading [13]. Zinc sequestration and detoxification occurs in the vacuoles when zinc is abundant providing a site of zinc storage [14-16]. Several gene transporters have been reported to be involved in zinc translocation to the vacuoles. Zinc induced facilitator (ZIF1) protein contributes to Zn and NA sequestration into the vacuoles thus removing the opportunity for both to be transported symplastically [17]. Two members of the heavy metal ATPase gene family (HMA) HMA2 and HMA4, mediate the release of zinc in root-to-shoot transport of zinc into the xylem [8,18]. Two members of the cation diffusion facilitator (CDF) family, the metal tolerance protein (MTP) MTP1 and MTP3, have been found to be involved in vacuolar zinc sub-cellular sequestration and contribute to zinc tolerance in A. thaliana [19,20]. NRAMP (natural resistance associated macrophage protein) are involved in transport of metals out of vacuoles [21]. Six members have been identified in Arabidopsis. NRAMP3 and NRAMP4 are involved in iron remobilization from the vacuole [22]. In a yeast complementation analysis only NRAMP4 was able to complement zrt1zrt2 growth under low zinc levels [23].

These gene families involved in mineral transport and sequestration represent some of the most obvious candidate genes for increased Zn seed levels in crops such as P. vulgaris. Very little is known about how Zn is transported from leaf xylem to phloem of developing seeds and ultimately unloaded into seeds [24].

Here we characterize the transcriptome of developing pod of two bean genotypes. These genotypes were shown to have similar Zn concentration in roots sampled during vegetative growth, leaves, and pods, but different levels of Zn in the seed. A total of 381 genes were differentially expressed including four genes that may play a role in Zn or other mineral transport to seeds including zinc-regulated transporter, iron regulated transporter like (ZIP), bZip transcription factor, zinc-induced facilitator (ZIF) and natural resistance associated macrophage protein (NRAMP) family genes. A total of eleven genes in the ZIF, NRAMP, YSL, and ferritin gene families contained SNPs between the two genotypes.

Materials and Methods

Plant Material

The two common bean genotypes used for this study both are navy beans from the Mesoamerican gene pool. Albion is a navy bean variety released by Asgrow in 1987, Voyager, is a navy bean released by Rogers Seed Company in 1995. These genotypes were selected based on their contrasting seed Zn concentration. Voyager has higher levels of seed Zn than Albion in diverse growing conditions [25]. In addition, in contrast to Voyager, Albion exhibits foliar Zn deficiency symptoms in low Zn and/or calcareous soils [26].

Plant Zinc Uptake Experiment

Seeds of Albion and Voyager were individually planted in 500 ml pots with 3:1 Sunshine Brand premium grade vermiculate and horticultural perlite grade (P.V.P. Industries, Inc). Treatments consisted of 0.5X strength modified Hoagland solution [27] with zinc added and without zinc (3 mMKNO3, 2 mM Ca (NO3)2 x 4H2O, sequestrene DTPA 10% Fe, 1.0 mM MgSO4 x 7H2O, 23.1 mM H3BO3, 0.38 mM ZnSO4 x 7H2O, 0.16 mM CuSO4 x 5H2O, 4.6 mM MoO4 x 2H2O, 1M KH2PO4 (pH to 6.0) was used as a fertilization treatment at a rate of 400 ml three times per week. Plants were grown in a growth chamber with a photoperiod of 16 hours light and 8 hours dark. Two replicates per plant were harvested as follows: roots and leaf tissue samples of vegetative plants were collected when the third trifoliate leaf had unfolded. Tissue of roots and leaves during flowering was collected when 30% of flowers were opened. Flowering was monitored daily and pods were collected 20 d after flowering. Seed was collected at physiologic maturity. Tissue was collected in liquid nitrogen and stored at -80°C. All samples were lyophilized and ground to powder with a Geno/Grinder 2000 (SpexCertiPrep, Metuchen, NJ) and zircon grinding balls. Plant tissue samples sent to A & L Laboratories (Fort Wayne, IN) for mineral analysis using induced coupled plasma spectroscopy. Mineral concentration was measured on 48 beans samples as follows: six tissue types, two Zn fertilization treatments, two genotypes, and two replications of each. Statistical significance was determined using proc glm and Tukey tests for pairwise comparisons in SAS for Windows v.9.2 (SAS Institute Inc., Cary, NC, USA).

Greenhouse RNA Seq Experiment

Seed of Voyager and Albion were planted in a greenhouse at Michigan State University. Two seeds were planted in 10 cm diameter clay pots filled with SureMix potting soil (Michigan Peat Company). Three pots were planted of each genotype and each pot was treated as a replication. Plants were watered as needed and fertilized with 0.5X Hoaglands solution (3 mM KNO3, 2 mM Ca (NO3)2 × 4H2O, sequestrene DTPA 10% Fe, 1.0 mM MgSO4 × 7H2O, 23.1 mM H3BO3, 0.38 mM ZnSO4 × 7H2O, 0.16 mM CuSO4 × 5H2O, 4.6 mM MoO4 × 2H2O, 1M KH2PO4 (pH to 6.0) biweekly starting at 20 d after germination. At anthesis, flowers were marked with a tag. At 12 days after anthesis individual pods were removed from plants and flash frozen in liquid nitrogen. Two pods per replication were ground to a fine powder with a mortar and pestle while completely frozen. Liquid nitrogen was continuously added to ensure tissue remained frozen throughout the grinding process. Total RNA was extracted from the samples using an RNA easy Plant kit (Qiagen, Cat. No. 74904). Following extraction, RNA samples were treated with RNase free DNase I (Qiagen Cat. No. 74254). RNA integrity and concentration was assessed for each of the samples using an Aligent 2100 Bioanalyzer (Agilent Technologies, Inc.). A subsample of the pod tissue was retained for mineral analysis. Following pod sampling, plants were grown to maturity and mature seeds were also analyzed for mineral concentration as described above and nitrogen concentration according to the Dumas method. Statistical significance was determined using proc glm and Tukey tests for pairwise comparisons (P<0.05) in SAS for Windows v.9.2 (SAS Institute Inc., Cary, NC, USA).

RNA Sequencing and Pre-Processing

Six RNA samples in total (3 replicates each of Albion and Voyager) were sequenced at the Michigan State University Research Technology Support Facility (RTSF) using an Illumina Genome Analyzer II (GA II). The library and flow cell preparation using kits and protocols from Illumina was conducted by the MSU RTSF. The sequencing was conducted as 75-bp paired-end reads. The RNA sequence was received from RTSF in FASTQ formatted files containing 75-bp paired-end reads. The file contained sequences and quality information about each sequence. The data were filtered using scripts from FASTX-Toolkit (FASTQ Quality Trimmer and FASTQ Quality Filter http://hannonlab.cshl.edu/fastx_toolkit/). FASTQ Quality Trimmer clipped the low quality ends with a quality threshold of 20 and removed the reads shorter than 64 bp. Subsequently, FASTQ Quality Filter script was used to remove low quality sequences with quality scores of 20 or less.

The P. vulgaris reference genome sequence v. 1.0 (DOE-JGI and USDA-NIFA and corresponding gene model annotation files were obtained directly from http://www.phytozome.net) [28]. A P. vulgaris genome index was built using Bowtie v. 0.12.7 [29]. Paired-end clean reads for each individual were aligned to the reference genome by TopHat v 1.4.1 [30] with the provided guidance of annotated gene models (GTF file). Bam files for each alignment are deposited in Gene Expression Omnibus (GEO).

Gene Ontology Enrichment Analysis

Mapped reads were evaluated for homology and annotated using Blast2GO software (http://www.blast2go.com). This process included three steps: 1) BLAST to find homologous sequences, with the following options, e-value threshold of E-10, non-redundant protein database (nr), high-scoring segment pairs (HSP) length cutoff 33. 2) MAPPING to retrieve (Gene ontology) GO terms and 3) ANNOTATION to select reliable functions, with e-value hit filter of 1E-6, cutoff 55, GO weight 5, Hsp-Hit coverage cutoff 0BLASTx sequence translation tool. We evaluated the functional enrichment level using the David Bioinformatics Database [31]. Significance of the enrichment analysis was declared when p-value<0.05.

Differential Expression Analysis

Abundance estimation was carried out using Cufflinks v1.3.0 with default parameters [32]. Normalization, estimated abundance, and tests for differential expression between tissue samples were performed using the program Cuffdiff with default parameters. The value used to compute significance of the observed change of transcript abundance was measured in fragments per kilobase of transcript per million mapped reads (FPKM) [33]. The differentially expressed genes were reported as fold change (on a log2 scale), and significance in terms of p-value and q-values. FPKM greater than zero was used as the threshold to consider a gene as being expressed. The significance of the gene expression difference was determined when the p-value was lower than the false discovery rate (q-value) (FDR<0.05). To validate the accuracy of expression profiles obtained by RNAseq, RT qPCR was performed on 6 genes with high or low expression levels as described in the supporting information. The genes belong to the ZIP gene family and the bZIP transcription factors involved in zinc homeostasis [34].

SNP Discovery and Validation

Single nucleotide polymorphisms between Voyager and Albion were discovered. Samtools mpileup [35] and bcftools were used for discovery and filtering using the following parameters: Alignments with a unique location, mapping quality >30, minimum depth = 8 and presence of SNP in >90% of reads on both forward and reverse strand. The options selected were–D (Output per-sample read depth),–u (Compute genotype likelihoods), and–f (The faidx-indexed reference file in the FASTA format). To validate SNPs called from the transcriptome sequence analysis, three genes YSL, HMA and ZIF were selected based on their role in Zn transport and presence of more than three SNPs. Primers were designed to amplify a template approximately 700 bp long. DNA was extracted from primary leaf tissue of Voyager and Albion using centrifugal filter “DNeasy Plant Kit” (QIAGEN) and quantified with Quant-iTPico Green dsDNA Assay kit (Invitrogen) following the manufacturer’s instructions. The mixture for each gene was optimized to contain 30 ng of DNA extract, 30 pmol of the primers and 0.5 U of Taq polymerase (AccuPrimePfxSuperMix, Invitrogen). After initial denaturation (95°C 5 min) 35 cycles (95°C 30 sec, 63°C 30 sec, 72°C 30 sec) of amplification were performed, followed by a final extension of 72°C for 5 min. Amplification products of both parents Albion and Voyager were visualized by electrophoresis on 2% agarose gel, stained with ethidium bromide and detected by ultraviolet trans illumination. PCR products with a single band were purified by ethanol precipitation, and directly sequenced via Sanger sequencing with the same primers used for PCR. Sequencing was conducted at the MSU RTSF. Sequenced products were compared and aligned with the reference common bean genome using Geneious (v. 5.6.2) (Biomatters).

Results and Discussion

Growth Chamber Experiment

Voyager and Albion are two navy bean cultivars. Voyager has been shown to contain higher seed Zn than Albion and these differences have been noted in both field and controlled environment experiments [25,26]. In order to determine if the difference in Zn levels between Voyager and Albion is limited to the seed, Zn levels were measured in roots, leaves, pods, and seeds of each genotype under two Zn fertilization treatments. Plants grown under the low Zn fertilization treatment had reduced Zn levels in the roots and leaves during vegetative growth, leaves during flowering, and pods for both genotypes. Zinc levels between Albion and Voyager were not significantly different for each tissue evaluated, except for seeds, where Voyager had 1.8 to 2.3 fold higher seed Zn levels (Fig 1). Zn content on a per seed basis was also calculated to determine if the seed Zn differences between the genotypes were due to seed size differences. On a per seed basis, Voyager had 3 times and 1.6 times more Zn than Albion under normal and no Zn fertilization respectively (data not shown). The higher seed Zn content in Voyager indicate that the differences in seed Zn concentration observed are not due to a dilution effect because of seed size. The importance of seed size in influencing seed micronutrient concentrations has been observed in genetic studies with Medicago truncatula and P. vulgaris [36,37]. Previous physiological and fertilization studies with Albion and Voyager showed that Albion accumulated higher Zn levels in stems, leaves and pod walls than Voyager. Voyager also had higher seed yield than Albion under low Zn fertilization and similar seed yield under normal and high Zn fertilization treatments [26]. Based on these findings, it appears that Voyager is better able to remobilize Zn within the plant and transport Zn to seeds. Since very little is known about genes involved in Zn transport into developing seeds, and it appears that root and leaf Zn levels are similar between the two genotypes, we decided to study the transcriptome of developing pods for potential clues on genes involved in Zn remobilization to the seeds. The seed Zn concentration differences between Voyager and Albion have been shown to be controlled by a single gene [25] thereby making these two genotypes excellent candidates for transcriptome analysis to identify genes responsible for the seed Zn differences.

Fig 1

Zinc concentration of roots, leaves, pods and seed of two bean genotypes grown under normal Zn and no Zn fertilization.

Pod Transcriptome Characterization

Voyager and Albion were grown in a replicated greenhouse experiment. Under these growing conditions Voyager had 1.4 times more seed Zn than Albion (Table 1). These values observed in greenhouse grown plants are similar to what has been observed for these genotypes in some field studies in Michigan (Cichy, unpublished). Seed Fe levels were variable in the greenhouse grown plants and therefore no significant differences were detected. However other studies have shown Voyager to have higher seed Fe than Albion [38]. Voyager also had 12.5% more N in the seed than Albion (Table 1). Positive correlations between seed Zn and N have been found in a number of crops including wheat and beans [39,40].

Table 1

Mean concentration of zinc, iron, and nitrogen in pods and seeds of Albion and Voyager plants from which RNA samples for sequencing were taken.

	Pod		Seed
	Zn (ug g-1)	Fe (ug g-1)	Zn (ug g-1)	Fe (ug g-1)	N (%)
Voyager	36a	92a	42a	101a	3.19a
Albion	28b	75b	29b	88a	2.79b

Means followed by the same letter in a column are not significantly different at P = 0.05

Developing pods were collected from the greenhouse grown plants at 12 days after flowering. At this developmental stage, 84% of pod weight was the pod wall and 15% was the developing seed in Voyager and 66% was the pod wall and 34% was the developing seed in Albion. This developmental stage has been characterized as the time prior to seed filling and when nitrogen is accumulating in the pods [41].

Voyager and Albion pod RNA was sequenced as 75 bp paired end reads. RNA reads were mapped to the P. vulgaris genome sequence, v 1.0 DOE-JGI and USDA-NIFA http://www.phytozome.net) [28] with an average of 31,505,836 high quality mapped reads per sample. A total of 27,197 genes were detected in the transcriptome of developing pods. Analysis of the number of transcripts per chromosome showed that chromosomes 1, 2, 3, 7, 8, and 9 contained the highest number of transcripts. Chromosome 10 had the lowest number of transcripts. In dry bean mapping populations, few quantitative trait loci have been found on chromosome 10 and it is often difficult to identify polymorphic markers on this chromosome, indicating a low rate of recombination [42]. A positive correlation (r = 0.5; P<0.05) between size of the chromosome and number of transcript was discovered. Chromosomes 2, 3, 7, and 11 had the highest number of highly expressed genes in the pods (Fig 2).

Fig 2

Number of transcripts and their distribution along of eleven chromosomes of common bean.

Gene annotation was achieved using Blast2Go (using BLASTx and E-value -6 as parameters). Using gene ontology, graphs were developed which classify gene expression in the bean pods. Transcripts expressed during this pod developmental stage were most related to oxidation reduction, auxin biosynthesis and amino acid phosphorylation (Fig 3a). Transcripts functioning in the nucleus and plasma membrane were the most represented cell types in the developing bean pods (Fig 3b). Transcripts related to ATP binding and protein binding were the most highly abundant gene types (Fig 3c). Enrichment analysis of differential expression genes determined their association to a specify pathway [31]. Genes were clustered into nine GO categories according to the David bioinformatics database as shown in Fig 4.

Fig 3

Histogram representation of GO classification.

Transcripts of developing pods were annotated in three categories: biological processes (3a), molecular function (3b) and cellular components (3c).

Fig 4

Enrichment analysis of 381 differential expressed genes between Albion and Voyager.

X axis indicates GO terms where genes were grouped, left Y axes indicate the number of genes and the right Y axis indicates p-values show significance of association of genes with GO term.

Genes highly expressed in developing pods included several seed maturation proteins, acid phosphatase and lipid synthesis related genes. The most highly expressed genes in the transcriptome of developing pods of dry beans were haloacid dehalogenase IIIB acid phosphatase (HAD), lipid transfer protein 3 (LTP3), PAP85 cupin protein, PRXR1 peroxidase protein, and phaseolin. The HAD IIIB acid phosphatase belongs to a family of plant phosphatases. Some members of this family have been annotated as vegetative storage proteins (VSP) highly expressed in soybean leaves and also abundant in A. thaliana flowers [43]. LTP3 (lipid transfer protein 3) transfer several different phospholipids and can bind fatty acids and could play a major role in membrane biogenesis [44]. PAP85 cupin family are metalloenzymes with two motif conserved sequences which act as ligands for the binding of an active-site metal ion, such as Fe, Mn, or Zn [45]. Additionally, they encode seed storage proteins and are involved in the regulation of nitrogen utilization [46]. Phaseolin is the seed storage protein most abundant in P. vulgaris seeds [47]. It is a glycoprotein formed by two genes, the α and β phaseolin genes with a relative electrophoretic diversity useful for discriminating geographical origin and wild and domesticated beans [47,48]. These results demonstrated that in P. vulgaris, storage protein transcripts were the most abundant process beginning at R5 and R6 growth stage which involve storage product accumulation, phases of cell expansion and synthesis of reserve metabolites [49].

Differential Expression

The expression levels of genes in developing pods were measured in fragments per kilobase of exon model per million mapped reads (FPKM). With this criterion,19,510 expressed genes were detected in Albion and 19,527 expressed genes in Voyager. We determined the correlation between normalized sequencing counts and RT qPCR gene expression data and there was a high correlation between the two methodologies (r = 0.85 p = 0.05) (Table A in S1 File).

The distribution of gene expression values in log10 was left-skewed; the median and mean FPKM values are 11.11 and 45.43 respectively. There were 381 differentially expressed genes in Albion and Voyager of which 6 were exclusively expressed in Albion and 3 genes were only expressed in Voyager. Genes with the highest differential expression patterns between the two genotypes included cysteine proteinases and MLP-like protein-43 more highly expressed in Albion pods (Table 2). These genes are related to growth and mobilization and accumulation of storage proteins in seeds in development [50]. In Voyager, the genes most highly differentially expressed as compared to Albion were cinnamoyl-CoA reductase (CCR-like) and 2Fe-2S ferredoxin-like. Both of these genes are related to metal ion transport in addition to abscisic acid biosynthesis and the electron transport chain. Differential expression analysis showed that the cinnamoyl-CoA reductase (CCR-like) and 2Fe-2S ferredoxin-like genes were 1.65 and 1.42 respectively more expressed in Voyager pods than Albion pods.

Table 2

Genes and function of the most highly differential expressed in Albion and Voyager.

Genotype	Function	Albion a	Voyager a	log2 (Fold Change) b , c
Albion
Cysteine proteinases	Extracellular proteinase probably having a crucial role during rapid cell growth and leaf expansion	791	203	-2
MLP-like protein 43	Associated with fruit and flower development and pathogen defense responses	732	98	-3
SCR-like 11	S locus cysteine-rich protein	392	66	-3
Serine carboxypeptidase-like	Serine-type carboxypeptidase activity involved in proteolysis	237	81	-2
Aspartic proteinase A1	Encodes an aspartic proteinase that forms a heterodimer and is stable over a broad pH range	231	76	-2
Low-molecular-weight cysteine-rich	Predicted to encode a PR (pathogenesis-related) protein.	215	31	-3
Voyager
CCR-like	Cellular cation homeostasis, divalent metal ion transport. Expressed in embryo axis, cotyledons.	146	460	2
2Fe-2S ferredoxin	Abscisic acid biosynthetic process, electron transport chain, pentose-phosphate shunt.	110	295	1
Zinc-binding ribosomal protein	DNA recognition, RNA packaging, transcriptional activation, regulation of apoptosis, protein folding and assembly, and lipid binding.	74	169	1
Adenine nucleotide alpha hydrolases	Function unknown. Involved in response to stress. Expressed during petal differentiation and expansion stage.	65	162	1
Bifunctional inhibitor/lipid-transfer protein/seed storage 2S albumin superfamily protein	Function in lipid binding. Located in endomembrane system. Expressed in shoot apex, embryo, flower, leaf, seed. Expressed during cotyledon expansion stage.	35	161	2
Basic chitinase	Defense response after wounding or pathogenic attack	10	150	4

a The value used to compute significance of the observed change of transcript abundance was measured in fragments per kilobase of transcript per million mapped reads (FPKM.)

b Log2 (Fold change): negative values indicate that the gene is more highly expressed in Albion and positive values indicate the gene is more highly expressed in Voyager.

c All genes were significantly differential expressed FDR<0.05

Gene families known to be involved in Zn transport were identified and their expression in the developing pods of Albion and Voyager was quantified (Fig 5). These families included ZRT and IRT–like protein (ZIP), basic region/leucine zipper motif (bZIP) transcription factors, vacuolar iron transport (VIT), natural resistance-associated macrophage protein (NRAMP), zinc induced facilitator (ZIF), yellow stripe (YSL), heavy metal ATPase (HMA), nicotianamine synthase (NAS), dehydrin, and metallothionein. The ZRT and IRT–like protein (ZIP) family is involved in uptake, transport to leaves and translocation to seeds, embryo, endosperm, and seed coat of zinc [51]. It was the largest family and 20 out of 23 members [34] were found in the developing pod transcriptome of which fifteen were expressed (Table 3). Expression analysis of this family in dry bean showed that some members were highly expressed in leaves and pods under two Zn treatments [34]. Additionally, in Arabidopsis and maize, some members of this family are preferentially expressed in the embryo and endosperm [52]. These results emphasize the importance of these genes in Zn transport into sink organs. bZIP transcription factors bZIP19 and bZIP23 in Arabidopsis, were associated to promoter regions of the zinc deficiency-induced ZIP4 gene of ZIP family [9]. These bZIPs belong to group F, have been described containing two DNA binding domains needed to respond to low zinc supply in Arabidopsis [9]. Two bZIP basic leucine-zipper transcription factor genes homologous to bZIP23 in Arabidopsis were identified and both were expressed in pods.

Fig 5

Expression levels of the ten gene families involved in Zn transport and accumulation identified in developing pod transcriptome.

Bars represent each gene family and colors correspond to different member within each family. Height of each bar represents gene fold change in Albion and Voyager. ZIP (zinc-regulated transporter, iron regulated transporter like), NRAMP (natural resistance associated macrophage protein), bZIP (bZip transcription factor), ZIF (zinc-induced facilitator), HMA (heavy metal associated), VIT (vacuolar iron transport), YSL (yellow stripe like), NAS (nicotiamine synthase), FRO (Iron reductase) and ferritin.

Table 3

Gene families involved in Zn and/or Fe transport and expression analysis in the developing pods of Albion and Voyager.

Gene family	Homologous A. thaliana	Chromosome	Position	Albion a	Voyager a	fold change b	FDR c
PvbZIP1	bZIP23 basic-leucine zipper	Chr05	3,212,438	22	23	0.0
PvbZIP2	bZIP23 basic-leucine zipper	Chr11	3,134,439	41	46	0.2
PvNRAM1	ATNRAMP, metal ion transporter	Chr01	44,116,444	26	31	0.3
PvNRAM2	ATNRAMP3, metal ion transporter	Chr02	1,609,575	61	62	0.0
PvNRAM3	ATNRAMP3, metal ion transporter	Chr03	46,129,963	12	11	-0.2
PvNRAM4	ATNRAMP6, metal ion transporter 6	Chr05	40,351,734	17	17	0.0
PvNRAM5	ATNRAMP, metal ion transporter	Chr07	37,134,084	30	34	0.2
PvNRAM6	ATNRAMP2, metal ion transporter 2	Chr09	11,751,007	19	23	0.3
PvNRAM7	ATNRAMP6, metal ion transporter 6	Chr09	18,914,511	14	12	-0.3
PvNRAM8	ATNRAMP2, metal ion transporter 2	Chr10	37,315,780	0	0	0.0
PvNRAM9	ATNRAMP6, metal ion transporter 6	Chr10	42,893,083	2	0	-2.3	*
PvZIF1	ZIFL2, zinc induced facilitator-like 2	Chr02	21,890,013	1	19	4.6	*
PvZIF2	ZIFL1, zinc induced facilitator-like 1	Chr05	1,050,386	3	6	0.9
PvZIF3	ZIFL1, zinc induced facilitator-like 1	Chr11	44,602,386	0	0	0.0
PvZIF4	ZIFL1, zinc induced facilitator-like 1	Chr11	44,656,239	0	0	0.0
PvZIF5	ZIFL1, zinc induced facilitator-like 1	Chr11	44,662,432	0	0	-0.4
PvZIF6	ZIFL1, zinc induced facilitator-like 1	Chr11	46,613,123	9	11	0.2
PvZIF7	ZIFL1, zinc induced facilitator-like 1	Chr11	46,625,565	0	0	0.0
PvZIF8	ZIFL1, zinc induced facilitator-like 1	Chr11	46,638,452	0	0	0.0
PvZIF9	ZIFL1, zinc induced facilitator-like 1	Chr11	46,652,668	4	3	-0.4
PvZIF10	ZIFL1, zinc induced facilitator-like 1	Chr11	46,667,766	18	8	-1.1
PvZIP1	ATZIP4, zinc transporter 4 precursor	Chr01	3,438,922	3	3	0.1
PvZIP2	ZIP10, zinc transporter 10 precursor	Chr02	19,642,778	0	0	0.1
PvZIP3	ZIP metal ion transporter family	Chr02	33,721,809	19	20	0.0
PvZIP4	ZIP10, zinc transporter 10 precursor	Chr03	49,001,484	0	0	0.0
PvZIP5	ZIP10, zinc transporter 10 precursor	Chr03	49,013,792	0	0	0.0
PvZIP6	ZIP11, zinc transporter 11 precursor	Chr05	37,425,474	0	0	0.0
PvZIP7	ZIP11, zinc transporter 11 precursor	Chr05	37,429,894	1	2	1.2
PvZIP8	ZIP1, zinc transporter 1 precursor	Chr05	5,642,976	6	8	0.4
PvZIP9	ZTP29, ZIP metal ion transporter family	Chr05	37,714,954	18	17	-0.1
PvZIP10	ZIP11, zinc transporter 11 precursor	Chr06	17,173,381	12	21	0.8
PvZIP11	ZIP1, zinc transporter 1 precursor	Chr06	199,508	1	1	-0.2
PvZIP12	ZIP5, zinc transporter 5 precursor	Chr06	1,040,877	2	8	2.2	*
PvZIP13	ATZIP6, metal ion transporter family	Chr06	18,953,200	7	8	0.2
PvZIP14	ZIP metal ion transporter family	Chr08	7,633,778	11	12	0.1
PvZIP15	ZIP5, zinc transporter 5 precursor	Chr08	59,348,008	8	15	0.9
PvZIP16	ATZIP6, metal ion transporter family	Chr08	57,181,379	21	24	0.2
PvZIP17	ATIRT3, iron regulated transporter 3	Chr09	12,668,955	31	36	0.0
PvZIP18	ZIP metal ion transporter family	Chr10	9,814,851	2	2	0.1
PvZIP19	ZTP29, ZIP metal ion transporter family	Chr11	5,068,287	1	1	0.0
PvZIP20	ZIP10, zinc transporter 10 precursor	scaff	1,071	0	1	0.2
PvVIT1	ATVIT1, vacuolar iron transporter 1	Chr02	48,170,585	2	1	-0.5
PvVIT2	ATVIT1, vacuolar iron transporter 1	Chr02	48,175,491	14	18	0.4
PvVIT3	ATVIT1, vacuolar iron transporter 1	Chr02	48,252,436	12	14	0.2
PvVIT4	Vacuolar iron transporter (VIT)	Chr02	23,134,245	0	0	0.0
PvVIT5	Vacuolar iron transporter (VIT)	Chr02	36,507,752	0	0	-1.5
PvVIT6	Vacuolar iron transporter (VIT)	Chr02	36,521,460	0	0	-0.8
PvVIT7	Vacuolar iron transporter (VIT)	Chr02	36,533,751	0	0	0.4
PvVIT8	Vacuolar iron transporter (VIT)	Chr02	36,541,077	0	0	0.0
PvVIT9	Vacuolar iron transporter (VIT)	Chr04	27,416,262	10	19	0.9
PvVIT10	Vacuolar iron transporter (VIT)	Chr07	7,508,398	28	38	0.4
PvVIT11	ATVIT1, vacuolar iron transporter 1	Chr08	6,284,802	6	7	0.4
PvVIT12	Vacuolar iron transporter (VIT)	Chr08	49,131,247	6	6	0.0
PvVIT13	vacuolar iron transporter (VIT)	Chr09	8,164,280	8	7	0.0
PvVIT14	Vacuolar iron transporter (VIT)	Chr10	3,221,741	1	1	0.1
PvVIT15	vacuolar iron transporter (VIT)	Chr10	3,229,195	0	0	1.0
PvYSL1	ATYSL1, Yellow stripe like 1	Chr01	13,421,083	18	28	0.0
PvYSL2	YSL6, Yellow stripe like 6	Chr01	16,152,062	60	65	0.1
PvYSL3	YSL7, Yellow stripe like 7	Chr03	626,298	16	18	0.1
PvYSL4	YSL7, Yellow stripe like 7	Chr03	631,299	2	2	0.0
PvYSL5	ATYSL1, Yellow stripe like 1	Chr04	21,588,269	3	4	0.4
PvYSL6	YSL7, Yellow stripe like 7	Chr04	41,773,229	17	15	-0.1
PvYSL7	YSL7, Yellow stripe like 7	Chr06	20,249,225	0	0	0.0
PvYSL8	ATYSL3, Yellow stripe like 3	Chr08	40,137,873	4	9	0.0
PvYSL9	ATYSL3, Yellow stripe like 3	Chr09	9,292,230	62	72	0.0
PvHMA1	HMA5, heavy metal atpase 5	Chr02	29,860,709	11	5	-1.1
PvHMA2	HMA5, heavy metal atpase 5	Chr02	29,878,666	1	0	-0.6
PvHMA3	HMA5, heavy metal atpase 5	Chr02	45,175,820	0	0	0.0
PvHMA4	HMA5, heavy metal atpase 5	Chr02	45,187,687	0	0	0.3
PvHMA5	HMA6, PAA1P-type ATP-ase 1	Chr02	36,870,747	5	6	0.3
PvHMA6	HMA7, copper-transporting ATPase	Chr02	34,600,085	13	13	0.0
PvHMA7	ATHMA1, heavy metal atpase 1	Chr03	5,628,284	30	26	-0.2
PvHMA8	ATHMA1, heavy metal atpase 1	Chr03	46,285,474	26	25	-0.1
PvHMA9	ATHMA2, heavy metal atpase 2	Chr03	33,726,455	20	19	-0.1
PvHMA10	ATHMA4, heavy metal atpase 4	Chr09	35,288,969	1	1	-0.5
PvHMA11	ATHMA8, type ATPase	Chr09	13,120,413	8	7	-0.2
PvHMA12	HMA7, copper-transporting ATPase	Chr09	35,544,425	33	35	0.1
PvHMA13	HMA5, heavy metal atpase 5	Chr10	3,512,059	3	3	0.1
PvNAS1	ATNAS2, nicotianamine synthase 2	Chr01	48,680,147	0	0	-0.2
PvNAS2	ATNAS4, nicotianamine synthase 4	Chr05	6,792,803	0	0	0.0
PvNAS3	ATNAS4, nicotianamine synthase 4	Chr06	23,217,021	0	0	-0.1
PvDehydrin	Dehydrin	Chr04	44,048,043	2	1	-1.8
PvDehydrin	Dehydrin	Chr09	921,414	24	14	-0.8
PvDehydrin	Dehydrin	Chr11	49,334,383	0	0	0.3
PvMT	metallothionein 2A	Chr08	11,131,986	3262	3989	0.3
PvMT	metallothionein 2A	Chr10	1,509,086	669	847	0.3
PvMT	metallothionein 2A	Chr10	1,905,781	310	246	-0.3
PvMT	metallothionein 2A	Chr10	1,900,737	0	1	1.1

a The value used to compute significance of the observed change of transcript abundance was measured in fragments per kilobase of transcript per million mapped reads (FPKM)

b Log2 (Fold change): negative values indicate that the gene is more highly expressed in Albion and positive values indicate the gene is more highly expressed in Voyager

c *indicates genes were significantly differential expressed FDR<0.05

We identified fifteen vacuolar iron transport (VIT) genes which nine members were expressed in developing pods. Relatively low levels of expression (0 to 38 FPKM) was observed. In Arabidopsis, VIT1 has been found highly expressed in the developing seeds and mediates iron storage in the embryo [53]. Nicotianamine synthase (NAS) is a metal chelator that produces a nonproteinogenic amino acid which binding a variety of transition metals [54,55]. Metallothionein proteins bind transition metals and play a role in the homeostasis and detoxification of non-essential minerals [56]. Four metallothionein genes were identified in the developing pods and three of them were highly expressed (246 to 3,989 FPKM) but none were differentially expressed. In Arabidopsis, AtMT4a and AtMT4b have been suggested to be involved in Zn storage in seeds [57].

Very few of the Zn and/or Fe transport related genes were differentially expressed. Those that were differentially expressed included one gene each from NRAMP, ZIP and ZIF families (Table 3). A member of NRAMP was expressed in Albion and not expressed in Voyager (2.3 fold change). NRAMP family genes have been involved mainly in iron translocation [21,22]. However, complementation analysis showed that a member of NRAMP was able to complement growth in yeast [23]. PvZIF1and PvZIP12 were more expressed in Voyager than Albion in 4.6 and 2.2 fold change respectively. ZIP and ZIF family genes have been shown to be involved in transportation of minerals to the vacuole and transport to seeds in Arabidopsis [8,23]. This suggests that one possible reason why Albion has lower seed Zn that Voyager is because it is being moved to the vacuoles and not transported to the seed. However, for each Zn/Fe related family analyzed transcripts were from 1 to 72 of FPKM value which was relatively low as compared to those genes related to lipid synthesis and storage protein. Therefore future studies may need to target earlier and later developmental time points to better characterize genotypic differences in Zn remobilization.

The transcript sequences from Voyager and Albion were also analyzed for SNPs. A total of 12,118 SNPs were identified between the two genotypes (S2 File). In total 3,401 genes contained SNPs and on average there were 3.6 SNPs per gene. Sanger sequencing was used to validate SNPs in nine genes (Table 4). Of the gene families related to Zn and/or Fe transport, eleven genes were found to have SNPs with a total of 47 SNPs average of four SNPs per gene. In total, 15 of the SNPs resulted in an amino acid change (Table 5). Of the genes with SNPs, the same ZIF gene (Phvul.002G108300) which was more highly expressed in Albion than Voyager contained 10 SNPs in the coding region. This gene maps to chromosome 11 and while numerous studies have reported QTL for seed Fe and Zn concentration [37,58,59], this ZIF gene is not within those regions. It is interesting to note that the HMA gene (Phvul002G288300 and Phvul002G19000) maps to chromosome 2 in a region where a major QTL for seed Zn concentration has been identified in bean RIL populations from both Mesoamerican and Andean intra gene pool crosses [37,59]. The seed Zn differences in Albion and Voyager with single marker genetic analysis indicated that this trait associated with SSR markers BM154 and BM184 found on chromosome 9 [38]. Based on the physical position of these markers on chromosome 9 (1,856,660 1,718,891bp respectively) genes such as dehydrin and bZIP44 were found in the surrounding region (922,386 and 1,006,283). Dehydrins are responsible for osmotic stress from drought, cold, and high salinity but also binds metals reducing metal toxicity in plant cells under water-stressed conditions [60].

Table 4

Forward and reverse sequence for all primer pairs used to validate putative SNPs in genotypes Albion and Voyager.

P. vulgaris genome ID	Gene	Forward primer	Reverse primer	Product Size (bp)
Phvul.008G157800	PvYSL	CGCTATGTCGTAACACTTCTGCACC	TTTGTGCTTGCTGCCTTAGGTGGG	741
	PvYSL	GGGAAGGGCAGAAAAGCCTTCGAC	TTGCCTTGATCCTCGGTGATGGGT	874
Phvul002G19000	PvHMA	AACCTCTCACCGCGACCTCACTAC	CCCACACAACAACCCCATCGGAAG	1483
	PvHMA	AACCTCTCACCGCGACCTCACTAC	CATGTTCGCAGATCCACGGCGTAA	1131
	PvHMA	GACACGGCGGTTTTGCTGACTTTG	GCACTCTCTAATGCCTGCCCTCCT	746
	PvHMA	GTTGGTGCATCTCAGGGTGTGCTC	GGCCAATGGATGCTCACTATTCACCT	853
Phvul.011G173100	PvZIF	GCCCAGCATTGGGAGGCTATTTGG	AAGCCACATTCGGAACATGACCGC	1211
	PvZIF	CGTGACGTGTGCAATGATGCCACT	CGCACCAACAACACAAAACAGGGA	561
	PvZIF	AGGTGGTGCAGTGTGAGTGTTCCT	AGGCAAGTTACAGATTGAATTGGTTCCCT	886

Table 5

Identification of SNPs in genes that are members of Zn and/or Fe transport-related families, followed by the length of the CDS, genomic length, number of SNPs between Albion and Voyager, whether those SNPs validated via PCR and if the SNPs resulted in an amino acid change.

P. vulgaris genome ID	Family	CDS Length	Genomic Length	SNPs in CDS	Chr	SNPs confirmed a	AA change
Phvul002G288300	HMA	2,982	5,648	5	2		Syn, Syn, Syn, Syn, Phe/Ser
Phvul002G19000	HMA	2,958	4,150	6	2	*	Syn, Syn, Syn, Syn, Syn, Syn
Phvul.002G208800	HMA	2,835	24,712	1	2		Ala/Thr
Phvul.003G047300	HMA	3,564	9,174	4	3		Syn, Syn, Syn, Syn
Phvul.002G014300	NRAMP	1,524	2,952	4	2		Val/Leu, Syn, Syn, Syn
Phvul.008G157800	YSL	1,908	2,751	3	8	*	Met/Ile, Syn, Syn
Phvul.001G081600	YSL	2,031	6,538	2	1		Gly/Ser, Syn
Phvul.011G173100	ZIF	1,470	5,330	10	11	*	Gln/Leu, Asp/Glu, Arg/Pro, Ala/Val, Syn, Syn, Syn, Thr/Ile, Val/Ile, Gln/His
Phvul.011G189500	ZIF	1,470	4,785	5	11		Syn, Ile/Val, Syn, Syn, Syn
Phvul.008G221200	Ferritin	891	2,550	1	8		Glu/Lys

a: * indicates SNPs were confirmed by PCR amplification and sequencing.

Accumulation of minerals in the seed involved several complex and still unknown mechanisms. The ability to uptake and accumulate minerals, how much mineral is absorbed by roots, transfer into the shoots and leaves via the xylem, and translocation to seeds via the phloem are all potentially important genetic regulation points. It still unclear which mechanism is the most important step in terms of uptake, transport, remobilization and accumulation to determine where our effort to increase concentration of Zn in seed should focus. In Pisum sativum, Zn remobilization from vegetative tissues to the seeds has been measured and 75–95% of mineral content in pods was remobilized to the seed tissue [61]. In this study we reported the main genes that likely are related to Zn remobilization during the seed filling period.

RNA sequencing was used to identify members of mineral transporter gene families expressed during bean pod development. The comparative analysis of two closely related bean genotypes with different levels of seed Zn indicate which genes are differentially expressed and which contain SNPs. This information is useful to identify candidate genes for seed mineral biofortification and the most promising candidate from this study is the ZIF gene (Phvul.002G108300).

This file contains the description of validation of the RNA-seq results with quantitative RT-PC and supplementary tables.

Table A. Comparison of fold change of PvZIP and Pv bZIP obtained from RNA-seq and qRT PCR validation of PvZIP and Pv bZIP genes in developing pods. Correlation between normalized sequencing counts in RNA-seq and expression determined by RT qPCR showed high correlation between two methodologies (r = 0.85 p = 0.05). Table B. Primer list for gene expression analysis via RT-qPCR. Four member of the ZIP gene family and two member of the bZIP family were amplified via RT-qPCR to validate the expression profiles obtained by RNA-seq.

(DOCX)

Click here for additional data file.

List of SNPs found in the pod transcriptome of Albion and Voyager.

(XLSX)

Click here for additional data file.