Identification and Expression Analysis of Candidate Genes Involved in Carotenoid Biosynthesis in Chickpea Seeds

Plant carotenoids have a key role in preventing various diseases in human because of their antioxidant and provitamin A properties. Chickpea is a good source of carotenoid among legumes and its diverse germplasm and genome accessibility makes it a good model for carotenogenesis studies. The structure, location, and copy numbers of genes involved in carotenoid biosynthesis were retrieved from the chickpea genome. The majority of the single nucleotide polymorphism (SNPs) within these genes across five diverse chickpea cultivars was synonymous mutation. We examined the expression of the carotenogenesis genes and their association with carotenoid concentration at different seed development stages of five chickpea cultivars. Total carotenoid concentration ranged from 22 μg g−1 in yellow cotyledon kabuli to 44 μg g−1 in green cotyledon desi at 32 days post anthesis (DPA). The majority of carotenoids in chickpea seeds consists of lutein and zeaxanthin. The expression of the selected 19 genes involved in carotenoid biosynthesis pathway showed common pattern across five cultivars with higher expression at 8 and/or 16 DPA then dropped considerably at 24 and 32 DPA. Almost all genes were up-regulated in CDC Jade cultivar. Correlation analysis between gene expression and carotenoid concentration showed that the genes involved in the primary step of carotenoid biosynthesis pathway including carotenoid desaturase and isomerase positively correlated with various carotenoid components in chickpea seeds. A negative correlation was found between hydroxylation activity and provitamin A concentration in the seeds. The highest provitamin A concentration including β-carotene and β-cryptoxanthin were found in green cotyledon chickpea cultivars.


INTRODUCTION
Chickpea (Cicer arietinum L.) is one of the most important legume crops in the semi-arid tropics. Its worldwide production ranks second after common bean (FAOSTAT, 2012). It is considered as one of the most important food legumes in the developing countries because of its nutritional value and its capacity for symbiotic nitrogen fixation that can provide the entire crop demand for nitrogen (Jukanti et al., 2012). Two chickpea types, kabuli and desi, are commonly grown. Selection for traits like flower color and zero tannins in seed resulted in the evolving kabuli from desi chickpea in the Mediterranean basin (van der Maesen, 1972;Moreno and Cubero, 1978;Jana and Singh, 1993). Chickpea has a relatively small (740 Mb) diploid (2n = 2x = 16) genome (Arumuganathan and Earle, 1991) and the genome sequence of both kabuli (Varshney et al., 2013) and desi (Jain et al., 2013) types are available that makes chickpea a good case for legume genetic and genomic research.
Carotenoids are also precursors for apocarotenoids such as plant hormone abscisic acid, which is essential for plant growth and development (Kermode, 2005;Umehara et al., 2008). Apocarotenoid formation is mediated by carotenoid cleavage dioxygenases (CCD; Auldridge et al., 2006). Violaxanthin deepoxidase works as part of the xanthophyll (or violaxanthin) cycle and has a key role in the de-epoxidation of xanthophyll pigments such as violaxanthin (V) and antheraxanthin (A) into zeaxanthin (Z) (Misra et al., 2006;Chen et al., 2014). In the last step, violaxanthin is converted into allenic carotenoid neoxanthin by neoxanthin synthase (NSY; Welsch et al., 2008). The detail pathway of carotenoid biosynthesis can be seen in the report by Yan et al. (2010) andda Silva Messias et al. (2014). Each class of these enzymes seems to be responsible for the carotenoid metabolism at different and specific subcellular sites under both normal and stress conditions (Rubio et al., 2008).
Specific gene family members that are responsible for carotenoid content and composition during endosperm development have been well characterized in maize (Li et al., 2008;. Various approaches including metabolic and genetic engineering (Welsch et al., 2010;Kumar et al., 2012;Mintz-Oron et al., 2012;Nogueira et al., 2013;Ariizumi et al., 2014), advanced genomics and bioinformatics (Wurtzel et al., 2012), genomic-assisted selection (Campbell et al., 2014;Owens et al., 2014) and transcriptome analysis (Caroca et al., 2013;Frusciante et al., 2014) have been applied for studying carotenogenesis and improvement of carotenoid content in different crops. To develop chickpea cultivars with higher carotenoid concentration, information on the genetic basis of carotenogenesis in chickpea is substantial. The availability of genome assembly of chickpea (Jain et al., 2013;Varshney et al., 2013) provides a good source of information to identify the potential candidate genes involved in carotenogenesis.
The objectives of this research were first to identify candidate genes for carotenogenesis in chickpea through genome-wide analysis, and secondly to examine their expression pattern and their correlation with carotenoid concentration at different seed developmental stages across five diverse chickpea cultivars.

Candidate Gene Analysis
The sequences of genes involved in the carotenoid pathway were collected from Medicago truncatula and Arabidopsis thaliana genomes in the GeneBank database such as National Center for Biotechnology Information (NCBI BLAST R online database). The sequences were blasted against CDC Frontier genome assembly (Varshney et al., 2013) in order to retrieve the gene sequences from chickpea. The structure of the genes and proteins in the carotenoid and isoprenoid pathways were analyzed in chickpea, M. truncatula and A. thaliana. To evaluate the similarity of the genes that have more than two copy numbers, protein sequence alignment was done using Bio Edit sequence alignment editor software (Hall, 1999).

Re-sequencing and SNP Calling
Whole genome resequencing was done on four chickpea cultivars CDC Verano, CDC 441-34, CDC Jade, and CDC Cory. The Paired-end (PE) genomic DNA libraries were constructed from 1 µg of gDNA using Illumina TruSeq DNA PCR-Free HT Library Preparation Kits (Illumina, Inc) and sequenced on Illumina HiSeq 2500 using 2 × 125 chemistry.
We selected 32 genes from both carotenoid and isoprenoid pathways for SNP analysis. The genome assembly of CDC Frontier was considered as the reference and the sequences from the other four cultivars were compared to the CDC Frontier genome annotation. The 2 kb upstream and downstream of each gene sequence were chosen to cover the intergenic region for SNP discovery. Using SnpEff software, the SNPs among the five chickpea cultivars were identified and extracted (Cingolani et al., 2012).

Seed Sample Collection
Five chickpea cultivars with different cotyledon colors including CDC Frontier (yellow cotyledon kabuli), CDC Verano (green cotyledon kabuli), CDC 441-34 (yellow cotyledon kabuli), CDC Jade (green cotyledon desi), and CDC Cory (yellow cotyledon desi) were grown in 15 L pots filled with Sunshine mix #4 media (SunGrow, Seba Beach, Alberta, Canada) in the green house. The NPK fertilizer (20-20-20) was added to each pot (three g L −1 ) for three times after plants reached the height of 20 cm.
The flowers were tagged at the day of anthesis (before the flowers open completely) and three to four developing seeds (pods) were harvested at each growth stage at 8, 16, 24, and 32 days post-anthesis (DPA). Three biological replicates were planted for each cultivar and seed growth stage. Harvested pods were immediately floated in liquid nitrogen and were kept in −80 • C until RNA extraction.

Gene Expression Analysis
Different primer pairs were designed using Primer3 online program (http://bioinfo.ut.ee/primer3-0.4.0/primer3/) for q-PCR analysis and their accuracy was checked by the Primer-BLAST program (Table S4).
Total RNA was isolated from chickpea seeds of four developmental stages using hexadecyltrimethylammonium bromide (CTAB) according to the procedure described by Kannan et al. (2014). Extracted RNA samples were treated with DNase I (Life Technology, Invitrogen, USA) to remove any DNA contamination. The synthesis of first cDNA strand was performed using Sensi FAST TM cDNA Synthesis Kit (BIOLINE, USA). The real-time PCR assay was conducted using C1000 Touch TM Thermal Cycler (BIO-RAD, USA) using the Sensi FAST TM SYBR NO-ROX Kit (BIOLINE, USA). The relative expression was calculated using 2 (− Ct) method (Livak and Schmittgen, 2001). In order to find the most appropriate internal control, six housekeeping genes including actin 1 (Act1), elongation factor 1-alpha (Ef1α), glyceraldehyde-3-phosphate dehydrogenase (GAPDH), initiation factor 4a (IF4a), heat shock protein 90 (HSP90), and 18S ribosomal RNA (18SrRNA) ( Table  S4) were selected and examined for their expression (Garg et al., 2010). The 8 DPA stage of CDC Frontier was used as the reference sample and two technical replicates were used for each biological replication to minimize sampling errors. We applied UPMG method to develop a dendrogram based on K-means clustering with Cluster v3.0 program (Eisen et al., 1998). The gene expression patterns are presented as a heat map using Treeview v1.60 (Page, 1996). The gene expression levels are also presented as bar graphs ( Figure S2).

Carotenoid Measurement
The seed carotenoid was measured across three developmental stages including 16, 24, and 32 DPA of five chickpea cultivars using high-performance liquid chromatography (HPLC). We used 100 mg of fine powder from whole chickpea seeds for carotenoid extraction. The samples were premixed with 400 µl of (1:1 of methanol and DCM [dichloromethane]) and then were centrifuged for 15 min at 11,000 rpm. The supernatant was transferred to a new tube and 400 µl of 100% acetonitrile was added and centrifuged for 5 min at 11,000 rpm. All solution was mixed with 0.1% BHT (butylated hydroxytoluene) to minimize carotenoid oxidation.
Chromatography was conducted using the Agilent 1200 LC system with Chemstation software (Agilent Technologies, Santa Clara, CA, USA). Separation was done on Prodigy 5 µm (250 × 4.60 mm) column with the mobile phase 58:20:22 acetonitrile/dichloromethane/methanol flowing at 0.8 ml/min. One hundred microliter volume from each sample was injected in each run for 45 min. Detection of various components was done using a photodiode array detector monitoring at a 450 nm wavelength.
UV-visible spectra analysis and comparison of retention time with the authentic standard was used for carotenoid determination (Ashokkumar et al., 2014). The retention time was 3.8, 4.5, 5.6, 9.3, and 22.3 min for violaxanthin, lutein, zeaxanthin, β-cryptoxanthin and β-carotene, respectively. We used three biological replicates with two injections as technical repeats to improve the analysis by the HPLC. Results were converted into µg g −1 as carotenoid concentration.

Statistical Analysis
Pearson correlation analysis was done between transcript levels and carotenoid concentrations in chickpea seeds at different developmental stages. All the statistical analyses were done using SAS software (Version 9.1, SAS Institute Inc., Cary, NC, USA).

Carotenoid Biosynthesis Genes in Chickpea, Arabidopsis, and Medicago
The total number of genes for both isoprenoid and carotenoid pathways were 32 in chickpea and 26 in Arabidopsis. We found more copy numbers of some genes from both pathways in chickpea compared to Arabidopsis which can be explained by the differences in the genome size of these two species. Chickpea genome size of 740 Mbp (Varshney et al., 2013) is over five times larger than Arabidopsis (135 Mbp; Arabidopsis Genome Initiative, 2000).
The properties of carotenoid and isoprenoid genes in M. truncatula were similar to chickpea as both species share common legume family. Except for three genes, PSY, BCH, and ZEP, we found similar copy numbers for the rest of the carotenoid and isoprenoid genes in M. truncatula and in chickpea (Table S1). The estimated genome size of M. truncatula is 465 Mbp that is bigger than that of Arabidopsis (Bennett and Leitch, 2011).
In chickpea, the genes CCD1 and GGPPS2 represent the largest and smallest size with 15 Kbp and 1.8 Kbp in size, respectively. Four genes including GGPPS1, GGPPS2, LCYB, and NSY had only one exon, whereas ZEP1 and ZEP2 had the highest number of exons (16 exons). In Arabidopsis, the largest and smallest genes were PDS with 6.6 Kbp and GGPPS2 with 1.7 Kbp, respectively. Two genes GGPPS1 and GGPPS2 had only one exon, while ZEP1 was detected to have the highest number of exons (16 exons) in Arabidopsis (Table S1).
Domain analysis showed that the largest number of genes in the pathways have the Rossmann-fold NAD(P)H/NAD(P)(+) binding (NADB) domain in their structure. The main domain in PSY genes is Isoprenoid Biosynthesis Enzymes-Class 1 that is similar to GGPPS genes from isoprenoid pathway as observed in Arabidopsis and Medicago (data not shown). The properties of all the domains are listed in Supplemental information (Table S2).
Two genes including DXS and PSY each with four copy numbers were chosen for sequence similarity analysis in chickpea. Sequence alignment and similarity matrix indicated that DXS1 and DXS2 had the highest similarity (0.842) and DXS4 had the lowest level of similarity with other DXS genes. In case of PSY genes, PSY2, and PSY3 were highly similar (0.816) and PSY1 and PSY4 had the lowest similarity (0.484) within this group ( Figure S1, Table S3).

Sequence Analysis and SNP Identification
Sequence analysis was done on the 32 candidate genes (Figure 1) that play an important role in carotenoid biosynthesis. In total 476 SNPs were found in the upstream, exon, intron and downstream sequences of the genes across the five chickpea cultivars. A total of 17 SNPs (Table 1) was found in the coding region of 1-deoxy-D-xylulose-5-phosphate synthase (DXS1 and DXS2), lycopene β-cyclase (LBC), β-carotene hydroxylase (BCH1), zeta carotene isomerase (ZISO2), 1-deoxy-D-xylulose 5-phosphate reductoisomerase (DXR2), and geranylgeranyl diphosphate synthase (GGPPS2). Only two missense mutations, which code for different amino acid, were found in ZISO2, while the rest of the SNPs were synonymous mutation. Consequently, the amino acid Serine (S) changed to Proline (P) at position 15 and Phenylalanine (F) to Leucine (L) at position 147 in CDC Verano compared to the reference genome CDC Frontier.

Carotenoid Concentration
The highest concentration of lutein was observed in CDC Jade at 16 DPA. The lutein concentration dropped significantly at 24 and 32 DPA. The lowest concentration of lutein was observed in CDC Frontier. The concentration of zeaxanthin was highest in CDC Cory followed by CDC Jade. The zeaxanthin concentration was lowest and similar in three cultivars CDC Frontier, CDC 441-34 and CDC Verano ( Table 2). The highest concentration of β-carotene was observed in CDC Cory and CDC Jade at 16 DPA, but at the later stage, we observed the highest concentration in CDC Jade and CDC Verano. CDC Jade and CDC 441-34 showed the highest and lowest concentration of β-cryptoxanthin, respectively, in almost all stages. In addition, the highest and lowest concentration of violaxanthin was observed in CDC Jade and CDC Frontier, respectively. Total carotenoid was highest in CDC Jade and lowest in CDC Frontier. In general, the concentration of different types of carotenoid decreased from 16 to 32 DPA in all cultivars. Overall in the five cultivars, the concentration of lutein was the highest followed by zeaxanthin, β-carotene, β-cryptoxanthin and violaxanthin ( Table 2).

Gene Expression
One of the most important factors for the real-time PCR data analysis is finding a good internal control because all the results would be normalized based on the cycle threshold (CT) of the house-keeping gene. Among all the house-keeping genes tested in this study, the 18SrRNA showed constant expression within and among the cultivars at different developmental stages (data not shown); therefore, it was used as the internal control for the expression analysis.
The next two enzymes in the pathway, carotene isomerase (CRTISO) 1 and 2, were expressed predominantly in the kabuli type. Lycopene β-cyclase, lycopene ε-cyclase, and carotenoid dioxygenase were highly expressed at 8 DPA in all cultivars. On average the lowest expression of β-carotene hydroxylase 1 and 2 was observed in CDC Jade and CDC Verano cultivars at 8 DPA. Expression of violaxanthin de-epoxidase was higher in CDC Jade at 8 DPA than that of other cultivars and the lowest expression was observed in CDC Frontier. Expression pattern of the two zeaxanthin epoxidase was similar with higher expression at 8 DPA, except for CDC Verano in which its expression increased at 16 DPA. The highest expression of neoxanthin synthase was observed in CDC 441-34 at 16 DPA (Figure 2 and Figure S2).

Correlation Analysis
Correlation analysis revealed positive correlations between the expression of the genes in the primary steps of the carotenogenesis including zeta-carotene desaturase, zeta 1 | The list of 17 SNPs found in the coding region of DXS1 and DXS2, 1-deoxy-D-xylulose-5-phosphate synthase; LBC, lycopene β-cyclase; BCH1, β-carotene hydroxylase; ZISO2, zeta carotene isomerase; DXR2, 1-deoxy-D-xylulose 5-phosphate reductoisomerase; and GGPPS2, geranyl-geranyl diphosphate synthase with their chromosomal location, physical position (bp) and the differences between the reference (CDC Frontier) and the alternative across four other cultivars.   carotene isomerase, poly lycopene isomerase, and carotenoid concentration. A negative correlation was obtained between hydroxylation and cleavage activities and provitamin A concentration (Table 3).

DISCUSSION
The current study revealed that the majority of the SNPs within the genes involved in carotenoid biosynthesis resulted from synonymous substitutions. Similar patterns of variation were also found in genes of carotenoid biosynthesis in tomato, citrus, pepper, and carrot (Livingstone and Anderson, 2009). The highly conserved carotenoid biosynthesis genes across different species reflect that these genes are required for a wide range of end products essential for the plants. Any changes on these enzymes could have major deleterious effect on plant fitness. Traditionally, researchers mainly focused on the nonsynonymous mutations that can result in the changing of the amino acids and consequently protein function (Sauna and Kimchi-Sarfaty, 2011). However, recent studies indicated that the synonymous mutations may have an effect on mRNA splicing (Pagani et al., 2005), mRNA stability (Kimchi-Sarfaty et al., 2007), the efficiency of protein translation, and protein folding (Sauna and Kimchi-Sarfaty, 2011). Protein expression can also be influenced by synonymous SNPs as they are involved in regulating microRNA-mediated genes (Wang et al., 2015). Therefore, the synonymous SNPs identified in the diverse chickpea cultivars in this study may have potential functional significance in carotenoid biosynthesis. Further study might be needed to examine the detail of this mechanism. The two non-synonymous SNPs in ZISO2 in CDC Verano are good examples of mutations that changed the amino acids as a result. The enzyme ZISO2 is involved in isomerization activities in the carotenoid pathway. Its expression significantly correlated with various carotenoid concentrations. The green cotyledon kabuli, CDC Verano, has the highest carotenoid concentration among all the kabuli type which could be due to the mutations in ZISO2 in this cultivar. Variation of exon numbers among different copies of the carotenoid genes is interesting in this study. For example, among the different copies of the CRTISO gene the exon numbers varied from 5 to 13, while in PSY the exons varied from five to nine (Table S1). Exon-intron architecture is one of the mysterious issues in gene evolution (Zhu et al., 2009). It seems that the length and CG content of first exon and intron have an association with functional element in higher organism (Kalari et al., 2006). Kreimer and Pe'er (2013) also discussed that exon variant can affect gene expression. It is possible that the changes in gene expression between different copy numbers of a carotenoid gene like PSY in chickpea might be affected by the variation in exon number and size.
The 32 candidate genes involved in isoprenoid and carotenoid pathways were distributed across all 8 chromosomes of chickpea. Different copy numbers for some candidate genes exist in the chickpea genome (Figure 1). The PSY gene has been considered as a key regulator of the carotenoid pathway (Li et al., 2007(Li et al., , 2008. In domesticated maize, the PSY1 locus has been the target of selective sweep and it was reported that 6.6-27.2% variations of the seed carotenoid concentration are associated with the activity of this enzyme (Palaisa et al., 2004;Chander et al., 2008;Bai et al., 2009). Pozniak et al. (2007) showed that PSY1 is cosegregated with the 7B QTL in durum wheat (Triticum turgidum L. var durum) and confirmed the correlation of the gene with phenotypic variation for endosperm color. Cereals commonly have three homologs of PSY in their genome. Study on PSY family in maize showed that the expression of each member can be different among various tissues. Usually, the expression of PSY1 is higher in leaves and yellow endosperm, but PSY2 expression is significant in almost all tissues (Gallagher et al., 2004). The third family member, PSY3 is normally expressed in embryo and root especially under stress condition (Li et al., 2007). In chickpea, we found four members of PSY family, which may have a positive effect on carotenoid concentration across different cotyledon colors.
Significant and positive correlation between isomerization activities and carotenoid concentration ( Table 3) in chickpea indicated the key role of these enzymes in providing common precursor lycopene for carotenoid types. We also found that desaturation reaction had a positive correlation (tend to be significant) with carotenoid concentration (data not shown). In addition, the transcript levels of isoprenoid genes including 1-deoxy-D-xylulose-5-phosphate synthase (DXS3), 1-deoxy-D-xylulose 5-phosphate reductoisomerase (DXR), 1-hydroxy-2-methyl-2-(E)-butenyl 4-diphosphate reductase (HDR), geranylgeranyl diphosphate synthase (GGPS1) (Julliard and Douce, 1991;Julliard, 1992;Lichtenthaler, 1999), can significantly affect carotenoid concentration in plant , as carotenoid biosynthesis is a derivative of isoprenoid pathway (Cuttriss et al., 2011). The deficiency of ζ -carotene desaturase (ZDS) gene in sunflower (Helianthus Annuus L.) resulted in the accumulation of ζcarotene and the absence of β-carotene, lutein and violaxanthin in cotyledon (Conti et al., 2004). In many cases, genetic transformation or overexpression of phytoene synthase and desaturase to crop plants has resulted in higher provitamin A and total carotenoid concentration (Ye et al., 2000;Paine et al., 2005;Diretto et al., 2007;Naqvi et al., 2009;Kim et al., 2012). This study showed that CDC Jade cultivar with highest carotenoid concentration also had the highest expression level of ZDS.
The β-carotene hydroxylation converts provitamin A compounds into xanthophylls that have no provitamin A properties (Quinlan et al., 2012), so the challenge for breeders to develop cultivars with higher provitamin A concentration is limiting the β-carotene hydroxylation (da Silva Messias et al., 2014). Diretto et al. (2006) showed that silencing of β-carotene hydroxylase results in a higher β-carotene concentration in potato. In this study, we observed lower hydroxylation activities in CDC Verano and CDC Jade cultivars with high carotenoid concentration.
Cleavage activity is involved in apocarotenoid production that consequently reduces the carotenoid concentration (Auldridge et al., 2006). Also, the carotenoid cleavage dioxygenase has been considered as a negative regulator for β-carotene concentration in Arabidopsis and its expression had a direct effect on the degradation and turnover of carotenoids in seeds during maturity period (Rodríguez-Ávila et al., 2011;Gonzalez-Jorge et al., 2013). For example, functional analysis of carotenoid cleavage dioxygenase mutant in Arabidopsis (AtCCD1), showed higher carotenoid concentration in mature seeds of mutant cultivar than the wild type (Auldridge et al., 2006). The important effect of ZEP1 region on carotenoid concentration in 281 maize lines with various kernel colors has been discussed by Owens et al. (2014). It seems that the ratio between carotenoid production and its conversion into apocarotenoid, plays a key role in carotenoid concentration in chickpea. In CDC Frontier and CDC 441-34 cultivars this ratio is lower than CDC Cory and CDC Jade cultivars with higher carotenoid concentration. Based on our results the carotenoid levels decreased from earlier stage to later stage in chickpea seeds. The ABA concentration increases during seed development and it is positively correlated with embryo maturation (Wang et al., 1998). The legume seeds, which are mostly embryo, have a significant concentration of ABA (Goldberg et al., 1994). ABA is not the only cleavage product of carotenoids and there are more products like strigolactone (Booker et al., 2004). This process can explain the reduction of carotenoid during seed development. In addition, non-enzymatic activities including oxidative stress and lipid peroxidation involved in carotenoid degradation during seed desiccation (Sattler et al., 2004(Sattler et al., , 2006Mène-Saffrané et al., 2010). However, our knowledge regarding regulatory mechanism involved in carotenoid concentration and composition may not be complete as potential other mechanisms may occur (Shumskaya and Wurtzel, 2013;De Moura et al., 2013).
The two stay-green cultivars, CDC Jade and CDC Verano had the highest concentration of β-carotene among the Desi and the Kabuli types, respectively. Also, total carotenoid concentration is higher in pea and chickpea cultivars with green cotyledon color compared to cultivars with yellow cotyledons (Ashokkumar et al., 2014), and the same results were obtained as the case of lutein in pea (Holasovál et al., 2009).
Mutation on the stay-green gene (SGR) that involves in degradation of chlorophyll resulted in delayed senescence and consequently green cotyledon chickpea seed. In thylakoid membranes, carotenoids are in a complex of chlorophyll and protein of Photosystem I and II. For example, the PSI contains high β-carotene concentration and PSII is rich in lutein. In addition, the high correlation was found between lutein and chlorophyll concentration in green pea (Holasovál et al., 2009). Importantly the increased levels of carotenoid products normally associated with both the chloroplastic and the cytosolic pathways, including mevalonate, sterols, and squalene, as well as triacyl-glycerides (Kumar et al., 2012). For instance, PSY has a coexpression with genes involve in the synthesis of plastoquinone, NAD(P)H dehydrogenase, tiorredoxin, plastocianin, and ferredoxin (Meier et al., 2011). It seems that carotenoids in the complex of green tissues stay more stable than tissues with chlorophyll degradation. It is discussed by Ashokkumar et al. (2014) that higher lycopene cyclase activity is the main reason for higher total carotenoid in green cotyledon pea; however the cyclase activity in two green cotyledon chickpea was not higher than the other cultivars in this study. We believed that lower cleavage and hydroxylation activities have a significant effect on carotenoid concentration in green chickpea.
The total carotenoid concentration in chickpea dry seeds (9.08 µg g −1 ; Abbo et al., 2005) is higher than that geneticallyengineered "golden rice" endosperm (1.6 µg g −1 ; Ye et al., 2000), or red colored wheat (1.8-5.8 µg g −1 ; Kruger and Reed, 1988; United States Department of Agriculture, 2010). However, in Golden Rice2 the carotenoid concentration increased up to 23-fold (37 µg g −1 dry weight) compared to the original one (Paine et al., 2005). Carotenoid concentration is a trait with high heritability that mostly influenced by genetic (cultivar), and the effect of environment is not substantial for all carotenoid types (Abbo et al., 2010;Owens et al., 2014). Lutein is the main carotenoid type in chickpea (Abbo et al., 2005(Abbo et al., , 2010 and wheat (Ramachandran et al., 2010). Desi chickpea usually has higher lutein concentration than the Kabuli type likely due to higher grain weight in Kabuli type (Ashokkumar et al., 2015). It seems that the chromosomal linkage and pleiotropy can address the association between higher lutein concentration and low seed weight (Abbo et al., 2005). The carotenoid concentration obtained in our work was consistent with results from earlier studies for Kabuli and Desi chickpeas (Ashokkumar et al., 2014(Ashokkumar et al., , 2015. It was reported that carotenogenesis genes are active in photosynthetic organs under various light qualities and the levels of both chlorophyll and carotenoids will increase dramatically in de-etiolation period (Romer and Fraser, 2005;Toledo-Ortiz et al., 2010). In addition, the total concentration of xanthophylls in different plant species increases dramatically in a strong light condition that results in decreasing the ratio between " lutein (L) and the xanthophylls-cycle components, zeaxanthin, antheraxanthin, and violaxanthin (Z+A+V)" (Hirschberg, 2001). In order to control the light effect on carotenoid types, we grew all plants in the green house under controlled condition.
In conclusion, the carotenoid concentration in chickpea is under the control of various genes that are associated with the structure and function of the genes in the carotenoid pathway. New technology in genome sequencing has helped us to understand the details regarding the structure, polymorphism, copy number, and location of the genes involved in carotenoid biosynthesis in chickpea seeds. Along with genome sequence data, the variability of the expression pattern and carotenoid concentration in five cultivars revealed a logical relationship between genotypes and phenotypes in this study. We demonstrated that synonymous mutations may have functional effects on the expression pattern of the different genes involved in the carotenoid biosynthesis pathway.

AUTHOR CONTRIBUTIONS
MR conducted the experiments, analyzed, and summarized the results. MR, AD, and BT wrote and finalized the manuscript; BT conceived and directed the project.