Four parental genotypes included wild P. elatius JI64 from Turkey and cultivated Afghan landrace P. sativum JI92 both from John Innes Pisum Collection (Norwich, UK); wild P. elatius VIR320 (Bogdanova et al., 2012) from Vavilov Institute Research of Plant Industry (St. Petersburg, Russia) and cultivated P. sativum cv. Cameor from INRA France. Furthermore, 126 F_5:6 recombinant inbred lines (RILs) derived from JI64 and JI92 cross (North et al., 1989) were used to establish respective phenotypically contrasting (dormant vs. non-dormant, dehiscent vs. indehiscent) bulks. P. elatius VIR320 differs from other wild peas in relation to the absence of gritty and testa pigmentation, possibly as the result of being either semi-domesticate or hybrid between wild and cultivated pea with unknown origin (Bogdanova et al., 2012).

The seeds of four parental and 126 RILs of F₆ generation of JI64 (wild) × JI92 (cultivated) and reciprocal (RILs) were harvested from glasshouse grown plants (February–May 2015). Twenty five seeds per line were incubated in petri dishes (9 cm diameter) over two layers of medium speed qualitative filter papers (Whatman, grade 1) wetted with 3 ml of tap water and incubated in a 25°C incubator with darkness. Imbibition was scored at 24 h intervals based on changes in seed swelling and germination was determined based on the radicle breaking through seed coat. The percentage, Mean germination time (MGT), Timson index (TI), and Coefficient of Velocity (CV) were calculated over 7 days period. We have used these various mathematical measurements in order to more precisely describe germination process as shown by Ranal and Santana (2006).

Pod dehiscence was measured either by direct observation of the pods on the plant or by drying harvested pods at room temperature (Weeden et al., 2002; Weeden, 2007). In case of parental lines JI92 (domesticated, indehiscent pod) and JI64 (wild, dehiscent pod) are dehiscence/indehiscence obvious after pod maturing. On the other hand evaluation of RILs was difficult in some cases, as slight pressure on pods by fingers is necessary for opening. If slight pressure was enough to complete fruit opening the line was evaluated as dehiscent, if not this line was evaluated as indehiscent.

Samples of seed coat (JI64, VIR320, Cameor and JI92; at least five seeds per genotype) were dissected from dry seed and saturated with 2% sucrose solution under vacuum. Equal volume of cryo-gel (Cryomatrix Shandon) was added to samples and shake overnight. Saturated samples were mounted into cryo-gel on the alum chuck, frozen down to −25°C and cut in cryotome (Shandon SME, Astmoor, UK) into 12 μm transversal section (Soukup and Tylová, 2014). Sections were stained with toluidine blue (0.01%, w/v in water), alcian blue (0.1%, w/v in 3% acetic acid), aniline blue fluorochrome (Sirofluor; 0.01%, w/v in 100 mM K₂HPO₄ with pH 9), or Sudan Red 7B (0.01%, w/v) according to Soukup (2014). The presence of proanthocyanidins was evaluated by staining with vanillin (Gardner, 1975) and DMACA (Li et al., 1996). Callose immunodetection was performed according to Soukup (2014) using primary antibody toward (1,3)-β-glucan (1:100; Biosupplies Australia PTY Ltd) and anti-mouse IgG Alexa Fluor 488 secondary antibody (1:500; Invitrogen). Control samples were processed without the primary antibody. Sections were observed with an Olympus BX 51 microscope (Olympus Corp., Tokyo, Japan) in bright field, blue (Olympus WB filter—callose immunodetection) or UV (Olympus WU filter—aniline blue fluorochrome and DMACA) excitation. Unstained control sections were surveyed in bright field or UV-excited autofluorescence. Figures were documented with an Apogee U4000 digital camera (Apogee Imaging Systems, Inc., Roseville, CA, USA). Dry intact seeds were vaccuum dried and gold coated (Sputter Coater SCD 050, Bal-Tec) before imaging with scanning electron microscope (JSM-6380LV; JEOL, Tokyo, Japan). Twenty days (after flowering) old pods of parental lines (JI64 and JI92) as well as RILs of F6 were fixed in 2% formaldehyde and stored in 4°C for latter observation of pod suture. Samples were cut on hand microtones at thickness of 100 μm, and the resulting segments were stained 1% phloroglucinol (Sigma, USA) in 12% HCl (Soukup, 2014).

Testa was separated from the rest of the seed, crushed, and extracted using mixture of acetone:water (70:30, v/v) with addition of 0.1% ascorbic acid to achieve efficient extraction of polyphenolic compounds in wide range of polarity and structural diversity (adapted from Amarowicz et al., 2009). 0.5 ml of extract was dried under a stream of nitrogen and solid residue was dissolved in 0.5 ml of methanol. The samples were then analyzed by ultra-performance liquid chromatograph Acquity UPLC I-Class coupled to high resolution tandem mass spectrometer Synapt G2-S with ion mobility separation capability (Waters, Milford, USA). Chromatographic column Raptor ARC-18 (100 × 2.1 mm, dp = 2.7 μm, Restek) and mobile phases (MP) A: water + 0.1% formic acid, B: acetonitrile + 0.1% formic acid was used for separation of components present in seed coat extracts. Flow rate of mobile phase 0.2 ml/min was applied. Electrospray was used as ion source. Spray voltage 2.5 kV in positive and 1.5 kV in negative ion mode, were used, respectively. Process of the LC/ESI-MS method optimization and detailed setup of mass spectrometer will be provided in Válková et al. (in preparation).

Seeds of each genotype/line were mechanically disrupted and the seed coats were separated and pooled (four seeds per genotype). Description of the studied RILs is given in Table S1. Small pieces ~2 mm were fixed on MALDI plate using a common double sided adhesive tape. The samples were analyzed directly without application of a matrix. The prepared samples were analyzed using high resolution tandem mass spectrometer Synapt G2-S (Waters) equipped with vacuum MALDI ion source. For desorption/ionization a 350 nm 1 kHz Nd:YAG solid state laser was used. Details of LDI-MS setup and analytical parameters of hydroxylated fatty acids can be found in Cechová et al. (in preparation).

The obtained LC/ESI-MS and LDI-MS data were processed by MarkerLynx XS a software extension of MassLynx platform (Waters). The processed data matrix, i.e., after extraction, normalization and alignment of retention times (in case of LC/ESI-MS data), m/z-values and intensities of signals, were transferred to Extended Statistics (XS) module, EZinfo (Umetrics, Malmo, Sweden), and studied by principal component analysis (PCA) and orthogonal projections to latent structures discriminant analysis (OPLS-DA). Both PCA and OPLS-DA were used for reduction of data dimensionality. OPLS-DA is a multivariate statistical method employing latent variable regression developed as an extension of more frequently used partial least squares method (Trygg and Wold, 2002). Coordinates of particular samples and RT_m/z pairs (or m/z-values in the case of LDI-MS data) in appropriate biplots and S-plots were used for evaluation of dormant and non-dormant genotypes mutual segregation and significance of detected signals of metabolites. The procedure was adopted and modified from Kučera et al. (2017). The most significant markers were further studied by targeted MS/MS experiments to reveal their identity (Cechová et al., in preparation; Válková et al., in preparation).

For Massive analysis of cDNA Ends (MACE) each sample of parental genotype was composed by several pooled developmental stages of seed coat (2, 3, 4 weeks and older) as it is unknown at which stage putative candidate genes are expressed. Seven selected RILs forming each of dormant resp. non-dormant bulk were previously (at F₆ generation) and after mature seed harvest (of F₇ generation used for RNA isolation) tested for germination behavior (Table S1). Seed coat were dissected under stereomicroscope, immediately frozen in liquid nitrogen and stored at −70°C until use. Frozen seed coats or dorsal and ventral sutures of pods were ground to a fine powder with liquid nitrogen using sterile mortars and pestles. Total RNA was isolated from seed coat (~100 mg) using the BioTeke Plant Total RNA Extraction Kit (China) or NucleoSpin RNA Plant kit (Macherey Nagel) for pods, according to the manufacturer's instructions. Yield/quantity and purity was determined by using NanoDrop 2000 spectrophotometer (Thermo Scientific) and diluted in DEPC-H₂O to 100 ng/μl. Isolated RNAs were treated with DNase according to Baseline-ZERO™ DNase protocol (Epicenter). In case of parental genotypes (JI64, JI92, Cameor, and VIR320) four consecutive developmental stages of seeds (14–37 DPA) were taken each represented by 1.25 μg of total RNA. The RIL bulks were made of 1.425 μg of total RNA of each of seven lines. In case of RNA samples used for pod dehiscence study, two parental lines (JI64 and JI92) and two bulks of contrasting RILs (with dehiscent or indehiscent pods) were used. The bulk of dehiscent RILs was established from eight and bulk of indehiscent RILs from five lines using excised pod sutures of 10 and 20 days after flowering. Each of these four final samples contained ~1 μg of total RNA each. The integrity of the RNA samples was examined with an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, USA).

MACE libraries were generated using GenXPro's MACE kit (GenXPro GmbH, Frankfurt, Germany) as described in Zawada et al. (2014). Briefly, cDNA from 5 μg of total RNA was randomly fragmented and biotinylated 3′ ends were captured after binding to a streptavidin matrix. A library ready for high-throughput sequencing was prepared using TrueQuant adapters included in the kit. The library consisted of 50–700 bp-long fragments derived from the 3′-end of the cDNAs. The 5′-ends of the libraries were sequenced on a HiSeq 2000 machine (Illumina) with 100 cycles to generate the MACE tags, each tag representing one single transcript molecule. In total, 6 cDNA libraries were prepared and sequenced for seed dormancy, while 4 libraries for parents and two contrasting bulks for pod dehiscence, each providing over 10 million reads (Table S2).

After sequencing the reads are in 5′–3′ orientation. To remove PCR-bias, all duplicate reads detected by the TrueQuant technology were removed from the raw datasets. Low quality sequence-bases were removed by the software Cutadapt (https://github.com/marcelm/cutadapt/) and poly(A)-tails were clipped by an in-house python-script. The reads were aligned to reference sequences using Novoalign (http://www.novocraft.com). This tool maps reads to reference sequences depending on certain parameters (i.e., quality) and calculates thresholds for each assignment. The reference sequences consisted of all Pisum mRNA sequences from NCBI. We annotated these sequences to all Fabaceae proteins from Uniprot “http://www.uniprot.org/” by BLASTX to Swissprot (“sp|.,” good annotation) and afterwards to Trembl (“tr|…,” less good annotation) protein sequences. All reads that could not be mapped to Pisum mRNA sequences from NCBI were used for a de novo assembly to generate contigs denoted as “noHitAssembly_xxx” and annotated in the same way as the Pisum mRNA sequences from NCBI. Normalization and test for differential gene expression between the bulks were calculated using the DEGSeq R/Bioconductor package (Wang et al., 2010). Differential gene expression was quantified as the log₂ ratio of the normalized values between two libraries (log₂ FC). The p-value and correction for multiple testing with the Benjamini–Hochberg false discovery rate (FDR) were computed due to determining significance of gene expression differences in pairwise comparisons of libraries. Lists of Differentially Expressed Genes (DEGs) for three comparisons of contrasting phenotype (dormancy: wild × cultivated, RILs and their parents, dehiscence: RILs and their parents) were made based on combination of pairwise comparisons of log₂ FC ratio of the normalized values (log₂ FC > 2, log₂ FC < −2) and FDR (FDR <0.01) between all libraries of these groups.

Gene Ontology (GO) enrichment analysis and normalized gene expression data were used to identify function and relationships of differentially expressed genes (Young et al., 2010). The results of the GO analysis were then exported into the Blast2go for the final annotations. The annotations provided the fragments with blast hit with the appropriate gene ontology terms which were classified into three categories: biological process (BP), cellular components (CC), and molecular function (MF). The DEGs were subjected for their presence in the different Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The various enzyme activities and the DEGs involved in the KEGG pathways were revealed for each of the combinations.

Transcripts containing SNP with at least five reads in both samples that are homozygous distributed in dehiscent vs. indehiscent RILs bulks with only one false allele read in 100 in either bulks were considered dehiscence specific. SNPs were discovered using Joint-SNV-Mix (Roth et al., 2012). The output given by Joint-SNV-Mix was furthermore processed by GenXPro in-house software to filter the SNPs. A minimum coverage of 10 bp was needed to be identified as an SNP. To identify the genomic localization of the SNP the surrounding region of the SNP was assigned per blastn to the genome of M. truncatula “JCVI.Medtr.v4.20130313” from http://jcvi.org/medicago. The “snpviewer” a webtool from the “http://tools.genxpro.net” was used to visualize the data.

Gene-specific oligonucleotide primers were designed (Table S3) based on MACE consensus sequences using the FastPCR software (Kalendar et al., 2014). The expression of selected candidate genes was validated by quantitative real time PCR (qRT-PCR). RNA samples (treated with DNase) were reverse-transcribed with Oligo(dT)₁₅ primer (Promega) in a two steps reaction in final volume 40 μl. The qRT-PCR analysis was run on the CFX96™ Real-Time Detection System (Bio-Rad) using the SensiFast SYBR® No-ROX kit (Bioline) or LightCycler® 480 SYBR Green I Master kit (Roche) in case of pod dehiscence study. Primers were designed using FastPCR or Oligo Primer Analysis Software (Molecular Biology Insights, USA) and produced amplicons ranging from 77 to 220 bp (Table S3). Every PCR reaction included 2 μl cDNA (1:10 diluted cDNA), 5 μl 2 × SensiFAST SYBR mix or LightCycler® 480 SYBR and 400 nM of each primer in final volume 10 μl. The expression was studied at two developmental stages in four contrasting parental genotypes (JI64, JI92, Cameor, and VIR320) and in case of dehiscence study also contrasting RILs with dehiscent or indehiscent pods. The conditions for PCR were: 95°C for 2 min; 45 cycles of 95°C for 10 s, 55°C for 30 s, and 72°C for 20 s; followed by a melting curve of 65–94°C (recovered every 0.5°C held for 0.5 s). The specificity of primers was confirmed by melting curve and gel analysis of products. Quantification of transcript level was determined by CXF Manager Software (Bio-Rad). Actin or β-tubulin gene was used as a reference to normalize relative quantification (Ferraro et al., 2014) using the comparative Ct (2^−ΔΔCt) method. Changes in transcript were estimated as fold change relative to the expression in the genotype Cameor (younger stage) in case of seed dormancy study and genotype JI92 (younger stage) in case of pod dehiscence study.

The structure of pea pericarp follows the common arrangement in Fabaceae. The exocarp consists of thick walled epidermis, the relatively thick mesocarp is arranged in several layers of parenchyma and endocarp composed of lignified sclerenchyma on inside of which is thin-walled epidermis. Both exocarp and mesocarp are rich in pectins, as indicated with metachromatic staining of toluidine blue (Figure 11).

Histological analysis of the seed coat in M. truncatula revealed changes in cell wall thickness in the outer integuments throughout seed development (Verdier et al., 2013a). In Arabidopsis and Melilotus (legume), seed permeability was modulated by mutations affecting extracellular lipid biosynthesis (in Verdier et al., 2013a). Similarly, in M. truncatula, cells of the outer integument showed abundant accumulation of polyphenolic compounds (Figure 4); which upon oxidation may impact seed permeability (Moïse et al., 2005). Current knowledge about physical dormancy mainly comes from studies on morphological structure, phenolic content, and cuticle composition in legume species (reviewed in Smýkal et al., 2014). Morphological observation indicated that seed hardness was associated with the structure of palisade and cuticular layer (Vu et al., 2014) and presence or absence of cracks (Meyer et al., 2007; Koizumi et al., 2008). Other authors have proposed that the compositions of carbohydrates, hydroxylated fatty acids, or phenol compounds in seed coats control the level of permeability (Mullin and Xu, 2000, 2001; Shao et al., 2007; Zhou et al., 2010). Mullin and Xu (2001) found that the seed coat of an impermeable genotype had a high concentration of hemicellulose, essentially composed of xylans, which would reduce the hydrophilicity of the seed coat. We have found considerable structural and functional differences in testa properties between wild and domesticated peas. Contrary to Ma et al. (2004), who found small cuticular cracks in soft but not hard seeds of soybean, the surface was similar among used lines with only small discontinuities over the whole surface (Figure 3). However, we cannot exclude that these small fissures resulted from the SEM sample preparation, similarly to the above mentioned work. In the non-dormant genotypes subjected to imbibition, large fissures appear preferentially in the hilar region and strophiole (not shown). However, those are the most likely consequence of embryo imbibition and its volume increase and thus we do not expect those as primary sites of water entrance. There is a lack of detailed description of the primary pathway of water entry in pea as well as in the whole Fabaceae family, although the topic is thoroughly discussed (e.g., Baskin et al., 2000; Meyer et al., 2007; Ranathunge et al., 2010; Smýkal et al., 2014). It is thus not clear whether the hilar or strophiolar region is the primary entrance of water in the non-dormant genotypes as suggested also by McDonald et al. (1988) or Korban et al. (1981) in soybean and common bean (Agbo et al., 1987) or if the minor fissures present in the cuticle and properties of the outer part of palisade macrosclereids make the difference as suggested by Ma et al. (2004). Interesting and generally neglected feature of macrosclereids is a presence of autofluorescent, phenolics containing lipidic material in the terminal caps of macrosclereids above the light line, which is not directly connected with the cuticle (Figure S2). There is no detailed information on the nature of this material or its possible functional significance. Obviously, detailed structure of outer part of the palisade macrosclereids deserves future attention. Dormant genotypes have thicker macrosclereids palisade layer, which might contribute to the water impermeability of coats of dormant pea genotypes as suggested by Miao et al. (2001). However, thickness alone does not necessarily account for water impermeability (de Souza and Marcos-Filho, 2001). Our results also suggest different thickness of palisade layer among the dormant genotypes with the most dormant JI64 having the thickest palisade layer (Figure 2). Metachromatic staining with toluidine blue revealed that the non-dormant non-pigmented genotype Cameor exhibits high level of polyanionic pectins with exposed free carboxyl groups (Figure S3). On the other hand, the non-dormant, but well-pigmented JI92 showed lower abundance of pectins related anions, similarly to pigmented dormant genotypes. There is clear connection between toluidine blue stainability and seed coat pigmentation—anthocyanidin presence. Taken together results from toluidine and Alciane blue with vaniline and DMACA suggest that tannins in the seed coats of pigmented peas are probably bound to other compounds of the cell walls, changing the staining properties of cell walls. Nature of these linkages is unknown, but covalently bound tannins might be indicated as alkaline hydrolysis releases proanthocyanidins from the cell walls (Krygier et al., 1982). Such expectation might be further supported by presence quercetin-3-rhamnoside fragments from LC/MS and MALDI-MS study. The detected increase in metachromatic staining of testa cell walls after weak acid treatment might be indicative of a crosslinking of pectins with other compounds in dormant as well as JI92 genotype, which is released in acid environment. Interestingly, it was reported that condensed tannins might be released from linkages in acid environment (Porter, 1989). We can speculate whether the intensity of pectin—tannin crosslinking is associated with physical dormancy of pea. There are scarce references indicating possible role of tannins in cell wall polymer network and its properties (Pizzi and Cameron, 1986). There is some supposition for callose deposition in the light line area (e.g., de Souza et al., 2012). However, strong staining with aniline blue fluorochrome staining in the upper part of macrosclereids including the light line (Figure S2) cannot be attributed to callose. Such assumption is consistent with the specific callose antibody localization which led us to the conclusion that the signal for aniline blue is probably signal for one or more other compounds structurally similar to callose. The interaction of aniline blue fluorochrome with other 1,3- or 1,4-β-D-glucans was described by Evans et al. (1984) and this interaction depends on the degree of polymerization, and nature of substitution of the 1,3-β-D-glucan chain as well as on the concentration of phosphate in the staining solution (Evans et al., 1984). Deeper anatomical and histochemical analysis of seed coat and more reliable detection of primary entrance point of water during early rehydration phase is needed.

Seed development has been thoroughly studied in number of crops including legumes, especially with focus on embryo development (Bewley et al., 2013). In our study, we have made comparative transcriptomics analysis in order to dissect candidate genes/pathways associated with domestication imposed changes on seed coat properties. Temporal transcriptional changes during seed and pod development were studied in pigeonpea (Pazhamala et al., 2016), soybean (Aghamirzaie et al., 2015; Redekar et al., 2015), Medicago (Gallardo et al., 2007; Benedito et al., 2008; Verdier et al., 2013a; Righetti et al., 2015), peanut (Zhu et al., 2014; Wan et al., 2016), and pea (Liu et al., 2015) seeds but no study was made on comparison of wild progenitor and cultivated crop. Moreover, these studies analyzed either entire seed/pod or developing embryos, while in our work we have used excised seed coat or dissected pod suture. During the RNA isolation from the seed coat tissue, we have experienced great difficulties when working with wild pea seed samples. As reported earlier for pigmented soybean seeds (Wang and Vodkin, 1994) proanthocyanidins binds to RNA and prevent its extraction. We failed when using standard phenol/chloroform (McCarty, 1986) or guanidium thiocyanate (Chomczynski and Sacchi, 1987) methods, as well as common plant tissue RNA isolation kits.

There is problem of ambiguously mapped reads, which are the major source of error in RNA-Seq quantification (Robert and Watson, 2015). Short-read alignment is a complex problem due to the common occurrence of gene families. In contrast to RNA-seq, MACE methodology is derived from 3′ UTR end of transcript and each is represented by single molecule. The choice of quantification tool also has large effect, as these also differ in the way they handle aligned data and multi-mapped/ambiguous reads (Robert and Watson, 2015). In order to exclude or minimalize that identify DEGs are solely due to the genetic differences between contrasting parental genetic stocks we have used phenotypically classified bulks derived from RILs (Table S1, Figure S1). The concept of Bulk segregant analysis (BSA) was established as a method to detect markers in a specific genomic region by comparing two pooled DNA samples of individuals from a segregating population (Michelmore et al., 1991). Coupling BSA with the high throughput RNA sequencing has been shown to be an efficient tool for gene mapping and identification of differentially expressed genes (Chayut et al., 2015; Bojahr et al., 2016). One possible bottleneck of our analysis was bulking of several developmental stages into single MACE sample, where temporal and spatial expression might be hidden, resulting in differences between MACE and qRT-PCR data. Dynamic nature of gene expression both in spatial and temporal levels is clearly seen at qRT-PCR analysis of selected DEGs (Figures 9, 10). The bulking of various developmental stages, moreover from different genotypes (in case of RILs) is the source of imprecision which can result in masking of DEGs. The key to the successful use of BSA is precision of phenotypic assignment. Although some imprecision in phenotypic classification and comparable low number of RILs used for bulking (Table S1), transcriptomics analysis has provided valid results. As shown on heat-map (Figure 8), RIL bulks are indeed genetic mixture of parental genotypes. MACE method (Kahl et al., 2012) detects allele-specific SNPs and indels associated with the defined genotypes that can be instantly used in genetic mapping (Bojahr et al., 2016). As result, there was significant clustering of homozygous SNPs associated with seed dormancy (e.g., respective parental alleles) on Mt chromosomes 3 and 4, or chickpea chromosomes 5 and 7, respectively (not shown). These correspond to pea linkage groups (LG) III and IV. Using identical RIL mapping population and DARTseq markers, we mapped seed coat thickness to LG I, III, IV, and VI (unpublished) and percentage of seed germination to LGII. These indicate that there is likely more than single major gene involved in seed dormancy, acting at different stages (testa thickness, permeability).

Despite that detected DEGs between dormant and non-dormant pea seeds belongs to various GO and KEGG pathways, the largest number of annotated ones was found within phenylpropanoid and flavonoid pathways (Figure S4). These are involved in various activities such as UV filtration, fixing atmospheric nitrogen, and protection of cell walls (Zhao et al., 2013). Analysis of soybean mutant defective in seed coat led to identification of differentially expressed proline-rich and other cell wall protein transcripts (Kour et al., 2014). Moreover, this single gene mutation has resulted in differential expression of 1,300 genes, pointing out that complex series of events, many manifested at the transcript level, lead to changes in physiology, and ultimately structure of the cell wall. Similarly, we speculate that gene/-s causative of pea seed dormancy results in complex transcriptional and metabolomics changes. Recently, KNOTTED-like homeobox (KNOXII) gene, KNOX4, was found responsible for the loss of physical dormancy in the mutant Medicago seeds (Chai et al., 2016) resulting in differences in lipid monomer composition. These findings are in agreement with our data obtained by laser desorption-ionization mass spectrometric comparative analysis between dormant and non-dormant pea genotypes. Especially long chain hydroxylated fatty acids such as mono and dihydroxylated hexacosanoate, heptacosanoate, and octacosanoate were found in higher concentration in dormant peas compared to non-dormant ones (as shown in Figure 8 for dihydroxyheptacosanoate and dihydroxyoctacosanoate), implying that the presence of a greater proportion of hydroxylated fatty acids may provide a greater interconnectivity of cutin hydrophobic components improving its stability and impermeability for water as discussed also by Shao et al. (2007). As downstream targets of KNOX4 gene, several key genes related to cuticle biosynthesis were identified, such as the cytochrome P450-dependent fatty acid omega-hydroxylase and fatty acid elongase 3-ketoacyl-CoA synthase (Chai et al., 2016). We have not found homologs genes when searched within our DEGs set. It can be hypothesized that different genes have been altered in independently domesticated crops, although possibly acting on identical pathways. There are limited studies combining of transcriptomic and metabolomics analysis (Enfissi et al., 2010) or a combination of both techniques including proteomics (Barros et al., 2010; Collakova et al., 2013). Such integrative approach enables not just to identify transcript and metabolite changes associated with given process, but also to focus on biochemical pathways relevant to studied trait and possibly also delimit candidate genes. As shown in Medicago (Verdier et al., 2013a) and soybean (Ranathunge et al., 2010) the gene expression in seed coat is complex and dynamic. Since we could not currently annotate 10% of detected transcripts, it can be expected that with available pea genome this could be further improved. They might include putative target candidate/-s for seed coat permeability.

In legume seeds, there are three parts: the seed coat, the cotyledon, and the embryonic axis which, on average, represent 10, 89, and 1%, respectively, of the seed content. Seed coat pigmentation was shown to correlate with imbibition ability in several legumes, including common bean (Caldas and Blair, 2009), chickpea (Legesse and Powell, 1996), yardlong bean (Kongjaimun et al., 2012), faba bean (Ramsay, 1997), and pea (Marbach and Mayer, 1974; Werker et al., 1979). The presence of proanthocynidins (PAs) in seed coats can be assessed by the appearance of brownish coloration, which is the result of PA oxidation by polyphenol oxidase (Marles et al., 2008). In soybean (Glycine max), the recessive i allele results in high anthocyanin accumulation in the seed coat, resulting in dark brown or even black color (Tuteja et al., 2004; Yang et al., 2010). In contrast, the dominant I allele, which silences chalcone synthase (CHS) expression and hence blocks both anthocyanin and PA biosynthesis, results in a completely colorless seed coat. However, there is not simple relationship between the testa pigmentation imposing dormancy, as numerous cultivated pea varieties have colored testa yet do not display seed dormancy as illustrated by this study used JI92 landrace. Mendel's A gene beside flower color has pleiotropic effect including seed coat pigmentation (Hellens et al., 2010), yet these traits can be decoupled by recombination (Smýkal, unpublished). Second, B gene of pea encodes a defective flavonoid 3′, 5′-hydroxylase, and confers pink flower color, by control of hydroxylation of flavonoid precursors (Moreau et al., 2012). Neither this mutation results in alteration of seed dormancy. Comparably more is known on Arabidopsis and Medicago seed development owing to available mutants. Many of these mutations indicate the important role of proanthocyanidins and flavonoid pigments in testa development (Graeber et al., 2012) including effect on seed dormancy. The proanthocyanidins (PAs) received particular attention due to their abundance in seed coats (Dixon et al., 2005; Zhao et al., 2010) including pea (Ferraro et al., 2014). PAs are also known as the chemical basis for tannins, which are considered to be important part of physical dormancy in some species (Kantar et al., 1996; Ramsay, 1997). Flavan-3-ol-derived PA oligomers and anthocyanins are derived from the same precursors, proanthocyanidins (Lepiniec et al., 2006) and chemical diversity is introduced early in the pathway by cytochrome P450 enzymes (reviewed in Li et al., 2016a). Anthocyanidin synthase, anthocyanidin reductase, and leucoanthocyanidin reductase were studied at transcriptional level by Ferraro et al. (2014) in cultivated pea varieties and showed to be developmentally regulated. In our comparative transcriptome profiling we have not found any of these genes to be among DEGs, suggesting that differences in metabolites (quercetin, gallocatechin) found by chemical analysis are not at these steps of PA biosynthesis. Anthocyanidins are either immediately modified by glycosylation to give anthocyanins by anthocyanidin 3-O-glycosyltransferases (UGTs) or reduced to generate flavan-3-ols (such as epicatechin) by anthocyanidin reductase for PA biosynthesis (Xie et al., 2003, 2004). There are several described Medicago mutants defective in respective genes, resulting in reduced testa pigmentation (Li et al., 2016b), although the relationship to seed dormancy was not specifically investigated. Indeed, we have detected several differentially expressed UDP-glycosylases, two of them studied by qRT-PCR (Figure 9). The glycosyltransferase superfamily consists of 98 subfamilies and only few have been characterized so far. Only a few members of the UGT72 family been shown to have activity toward flavonoids; such as the seed coat-specific UGT72L1 from Medicago (Zhao et al., 2010) and several seed specific UGTs in L. japonicus (Yin et al., 2017). The UGT72L1 catalyzes (Zhao et al., 2010) formation of epicatechin 3′-O-glucoside (E3′OG), the preferred substrate for MATE transporters. MATE1 mutant display altered seed coat structure and PAs accumulation (Zhao and Dixon, 2011) and also has significantly lower seed dormancy levels (Smýkal, unpublished). The mechanism of PA polymerization is still unclear, but may involve the laccase-like polyphenol oxidase (Zhao et al., 2010). Notably, Medicago myb5 and myb14 mutants exhibit darker seed coat color than wild-type plants, with myb5 also showing deficiency in mucilage biosynthesis, and accumulating only of the PA content of wild-type plants. When myb5 seeds are exposed to water, they germinate readily without dormancy typical of wild type Medicago seeds (R. Dixon, personal communication and P. Smýkal, unpublished). All these observations suggest that PA oligomers play indeed a role in seed coat mediated dormancy.

Our LC/ESI-MS experiments confirm the presence of significantly higher contents of dimer and trimer of gallocatechin (i.e., soluble tannins of prodelphinidin type) in dormant compared to non-dormant pea genotypes. Catechin dimer and trimer was also found in the pea seed coat extracts but their differences between dormant and non-dormant peas are much less significant than their gallocatechin counterparts. This fact points out the significance of hydroxylation of B-ring of PAs in relation to dormancy. The insoluble PAs are the result of oxidative cross-linking with other cell components. Variation in PA content in the pea seeds has been reported (Troszyńska and Ciska, 2002) but not in comparison of wild vs. cultivated peas. PAs play also important roles in defense to pathogens, and because of the health benefits are of industry and medicine interest. PA biosynthesis and its regulation have been dissected in Arabidopsis using transparent testa (tt) mutants, which regulate production, transport or storage of PAs (Lepiniec et al., 2006), and 20 genes affecting flavonoid metabolism were characterized at the molecular level (reviewed in Bradford and Nonogaki, 2009). Many of these flavonoid biosynthesis pathway genes have been found to affect dormancy of Arabidopsis seeds, indicating the role of pigments in this process (Debeaujon et al., 2000). Similarly Medicago also synthesizes PAs in the seed coat, which consists essentially of epicatechin units (Lepiniec et al., 2006; Zhao et al., 2010). Polymerization of soluble phenolics to insoluble polymers is promoted by peroxidases (Gillikin and Graham, 1991) and catecholoxidases (Marbach and Mayer, 1974; Werker et al., 1979), which are abundant in legume seed coats. Positive correlation in content of phenolics, the requirement of oxidation and the activity of catechol oxidase in relation to seed dormancy (germination) in wild vs. domesticated pea seeds have been shown by Marbach and Mayer (1974) and Werker et al. (1979). Recently, epicatechin, cyanidin 3-O-glucoside, and delphinidin 3-O-glucoside were isolated in wild compared to cultivated soybean seed coats (Zhou et al., 2010) with epicatechin being in significant positive correlation with hardseededness.

Beside proanthocyanidin also flavonols, first of all quercetin derivatives, are frequently found in legumes including pea (Dueñas et al., 2004). We found significantly higher content of quercetin rhamnoside in dormant JI64 genotype compared to non-dormant JI92. Similarly, its hydroxylated analog, i.e., myricetin-3-rhamnoside appeared to be a marker of dormancy. This compound was found in many legumes including chickpea, horse gram (Sreerama et al., 2010) and pea (Dueñas et al., 2004). As described in review of Agati et al. (2012) the antioxidant properties of flavonoids represent a robust biochemical trait of organisms exposed to oxidative stress of different origin during plant-environment interactions (regulation of the action of reaction oxygen species, ROS). The effect of ROS on the plant developmental processes including seed germination was described by Singh et al. (2016) including increase of free radical scavenging during pea seed germination (Lopez-Amoros et al., 2006). Presence of phenolic compounds in seed coat might help to protect against fungal diseases during germination as shown in lentil (Matus and Slinkard, 1993).

Pod maturation might terminate with pod shattering, which is an important trait for seed dispersal of wild species but generally unwanted trait in crops (Fuller and Allaby, 2009). Central to the ballistic seed dispersal in Pisum is the dehiscent pod (single carpel fused along its edges) where the central pod suture undergoes an explosive rupturing along a dehiscence zone (Ambrose and Ellis, 2008). During pod shattering, the two halves of the pod detach due to a combination of the diminished cell walls adhesion in the dehiscence zone, and the tensions established by the specific mechanical properties of drying cells of endo and exocarp of the pod shell. These two principal aspects are shared among families producing dry dehiscent fruit, such as Fabaceae and Brassicaceae (Grant, 1996; Dong and Wang, 2015). The spring-like tension within the pod shell is generated during differential drying-induced shrinkage of endocarp and the outer part of the shell—exocarp (Armon et al., 2011). Properties of both the obliquely arranged rigid and lignified inner sclerenchyma as well as the pectin rich exocarp cell with longitudinal orientation are crucial factors of springing in Fabaceae pods. Regulators of their development affecting geometrical arrangement of the layers and their histological properties (cell wall thickness and composition, lignification, hydration) will be key factors generating required tension. Composition and characteristics of pod cell shell cell walls correlate with shattering of yardlong bean and wild cowpea (Suanum et al., 2016) and thickness of the shell and extend of sclerenchymatous dorsal bundle caps was connected with shattering in soy (Tiwari and Bhatia, 1995). The major QTL controlling pod dehiscence in soybean is qPDH1 (QTL for Pod Dehiscence 1). qPDH1 had been recently cloned and shown to encode a dirigent-like protein expressed in the sclerenchyma of differentiating endocarp and modulating the mechanical properties of the pod shell. Lignin biosynthesis is the most likely process affected by qPHD1 (Suzuki et al., 2010; Funatsuki et al., 2014), which might be connected with modulation of torsion within drying pod walls (Funatsuki et al., 2014). However, the precise biochemical activity of qPHD1 is still unclear. Similarly the sclerenchyma differentiation and lignification of the endocarp and valve margin cells of Arabidopsis are central to silique dehiscence with NAC SECONDARY WALL THICKENING PROMOTING FACTOR 1 (NST1) and SECONDARY WALL ASSOCIATED NAC DOMAIN PROTEIN 1 (SDN1) being identified as master regulators of their differentiation (Zhong et al., 2010). Lignification of endocarp and valve margin cells is lost together with dehiscence in the nst1snd1 double mutant (Mitsuda and Ohme-Takagi, 2008). The other decisive point of pod dehiscence is mechanical stability/instability of dehiscence (suture) zone which might trigger the explosive release of pod shell tension (Grant, 1996). Decreased cell to cell adhesion, might be carried out by action of endo-1,4-glucanases and endopolygalacturonases disintegrating the middle lamella in the separation layer (Christiansen et al., 2002). Degradation of pectin in the middle lamella of abscission zone is common theme in fruit shattering (Dong and Wang, 2015). Contrary, mechanism of pod shattering resistance due to reinforcement of the suture has been described in domesticated soybeans NST1/2 homologous transcription factor SHATTERING1-5 (SHAT1-5) from NAC family has been unveiled, inducing excessive secondary cell wall deposition and lignification in the outer part of the suture (fiber cap cells; Dong et al., 2014). Up-regulated expression of SHAT1-5 in domesticated soy “locks” the dehiscence zone interconnecting the vascular bundle caps of sclerenchymatous fibers in the suture vicinity preventing shattering.

The fundamental elements of fruit shattering regulatory network is being uncovered recently in Arabidopsis and homologous genes were identified also in other species and crops (Zea mays, Triticum aestivum, Oryza sativa, Glycine max, Sorghum bicolor, Sorghum propinquum, or Solanum lycopersicum; for review see Dong and Wang, 2015; Ballester and Ferrandiz, 2016). MADS box genes of Arabidopsis SHATTERPROOF1 (SHP1) and SHATTERPROOF2 (SHP2) participate in the dehiscence zone specification (Liljegren et al., 2000). INDEHISCENT (IND) b-HLH transcription factor act down-stream to SHP1/2 as regulator sclerenchyma differentiation in endocarp and valve margin. The shp1/2 double mutant as well as ind produces indehiscent siliques devoid of proper cell specification and differentiation in the dehiscence zone (Liljegren et al., 2004; Dong and Wang, 2015). SHP1 and SHP2 are required for the proper specification of the different cell types within the valve margin and the DZ and both genes probably represent the top of the hierarchy regulating DZ formation (Liljegren et al., 2000). SHATTERPROOF genes along with INDEHISCENT (IDEH) are the main regulators of establishment of lignified layer, which causes pod dehiscence. Another MADS box gene involved in dehiscence zone formation is FRUITFULL (FUL), which expression appears at the inception of the carpel primordia, and soon after becomes restricted to the cells that will give rise to the valves, a pattern that is complementary to that of the SHP genes (Liljegren et al., 2004; Dong and Wang, 2015). Significant up-regulation of SHATTERING1–5 (SHAT1–5) in the fiber cap cells (FCC) of cultivated soybean was shown to be responsible for the excessive cell wall deposition in the FCC, which in turn prevents the pod from committing dehiscence after maturation (Dong et al., 2014). Homologous genes defining dehiscence zone identity and its differentiation INDEHISCENT, SPATULA, SHATTERPROOF, bZIP, and SHATTERING (Ferrándiz et al., 2000; Girin et al., 2011; Dong et al., 2014) were identified in Pisum genome and their abundance was tested in RNA isolated from pod suture tissue of wild and domesticated pea as well as of contrast RILs incurred from crossing of wild and domesticate parent with dehiscent or indehiscent pods. In case of SHATTERPROOF and SHATTERING homologous genes, we found differences in expression between parental lines of pea but we didn't find the similar results in case of contrast RILs. The other homologous genes (INDEHISCENT, SPATULA, and bZIP) did now exhibit any significant difference in expression between dehiscent and non-dehiscent phenotypes.

In the legumes, the pod shattering trait is controlled by one or two dominant genes or QTL. In pea and lentil, genes controlling pod shattering map to a syntenic region, suggesting that the same genes may have been modified during the domestication of the two cool-season legumes (Weeden et al., 2002; Weeden, 2007). Single locus control of pod dehiscence was found in lentil (Ladizinsky, 1998), while two loci in mungbean (Isemura et al., 2012), yardlong bean (Kongjaimun et al., 2012), one controlling the number of twists along the length of the shattered pod, and second the percentage of shattered pods, similarly to two loci found in pea (Weeden et al., 2002; Weeden, 2007), and common bean (Koinange et al., 1996). Bordat et al. (2011) localized Dpo locus responsible for loss of pea pod dehiscence on LGIII. We obtained the similar result by the genome-wide DArTseq analysis. Based on comparison of our candidate gene position in M. truncatula genome and SNPs map (Tayeh et al., 2015) we localized our candidate genes in pea genome. In total of our 25 candidate genes 7 are localized on LGIII, 6 on LGVI, 4 on LGII, 3 on LGV, 2 on LGVII, 2 on LGIV, and 1 on LGI (not shown). MACE-P015, the main candidate gene possibly responsible for pod dehiscence localized on LGIII, is a homolog of peptidoglycan-binding domain protein (PGDB) of M. truncatula (Medtr2g079050). These proteins may have a general peptidoglycan binding function and this motif is found at the N or C terminus of a variety of enzymes involved in bacterial cell wall degradation. Many of the proteins having this domain are so far uncharacterized. Matrix metalloproteinases (MMP), which catalyze extracellular matrix degradation, have N-terminal domains that resemble PGBD (Seiki, 1999). On the other hand our candidate MACE-P015 has also 80% match with Cicer arietinum proline-rich extensin-like protein EPR1 (XM_004488673). Extensins are plant specific structural cell-wall proteins (Lamport et al., 2011); they can account for up to 20% of the dry weight of the cell wall and can significantly modulate mechanical cell wall properties through linkages to other cell wall component, which can play a role in pod dehiscence.

Study of biochemical and molecular mechanisms underlying plant domestication process is important area of research in plant biology. In the current study we used a comparative anatomy, metabolomics, and transcriptome profiling of pods and seed coats in wild and domesticated pea in order to identify genes associated with loss of seed dormancy as well as pod dehiscence. We have identified genes showing differential expression in respective parents as well as phenotypically contrasting RILs. Among others, there were number of genes belonging to phenylpropanoid pathway, which was also identified by metabolomics analysis of seed coat. Our results support the role of proanthocyanidins and their derivatives in physical seed coat mediated dormancy. One of the identified differentially expressed gene involved in pod dehiscence showed significant down-expression in dorsal and ventral pod suture of indehiscent genotypes. Moreover, this homolog of peptidoglycan-binding domain or proline-rich extensin-like protein mapped correctly to predicted Dpo1 locus on PsLGIII. This integrated analysis of the seed coat in wild and cultivated pea raised new questions associated with domestication and seed dormancy. Having underlying gene(s) in hands for various independently domesticated legume crops it would help our understanding of genetic and molecular processes involved in seeds dormancy. Moreover, extended knowledge on control seed dispersal and seed dormancy is necessary for diverse applications—biodiversity conservation as well as breeding.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.