Transcriptome Analysis of Stigmas of Vicia faba L. Flowers

1. Introduction

In angiosperms, pollination depends on a complex communication between the male (pollen grains) and the female (stigma/style) reproductive organs. In the compatible pollen-pistil interaction several events are involved: pollen capture, adhesion, germination, penetration of pollen tube into the stigma, growth of the pollen tube through the style and final entry of the pollen tube into the ovule. Stigmas can be generally classified into two main groups according to the presence (wet stigmas) or absence (dry stigmas) of a viscous secretion on the stigma surface [1,2]. Once pollen grains are transferred to the stigma by abiotic (e.g. water, wind) or biotic vectors (e.g insects, birds); or directly by contact between the anther and the stigma, pollen-pistil interactions differ between species. In species with wet stigma (e.g. Nicotiana tabacum, Lilium longiflorum), the cuticle is disrupted due to the presence of exudates, which can be composed by lipids, proteins, carbohydrates, phenols, glycoproteins, ions and enzymes such as esterases and peroxidases [3,4]. Unspecific pollen grains adhere to the stigma surface thanks to exudates and pollen hydration occurs passively transferring water from the stigmatic exudates [5]. By contrast, in species with dry stigmas (e.g. Arabidopsis thaliana, Zea mays, Oryza sativa) the events following pollination are species-specific and highly regulated [6]. Pollen adhesion and germination have been well studied in species with self-incompatibility systems such as Brassicaceae, Poaceae or Papaveraceae, and a high diversity of molecules and processes have been discovered (reviewed in [4,7]).

Previous proteomic and transcriptomic studies in species with wet and dry stigmas indicate that both strategies express unique as well as common genes and proteins during stigma maturation. It is expected that the evolution of genes involved in sexual reproduction occurs at a higher rate than those in charge of background processes; moreover, those genes responsible for the maintenance of species boundaries will be species-specific and therefore different between species [5,8]. Allen et al. [9] found that certain gene families were consistently found in pistil tissues of different species such as cytochrome P450, ATP-binding cassette (ABC) transporters, lipid transfer proteins (LTPs), zinc finger proteins, extensin-like proteins, receptor protein kinases, disease resistance proteins or nodulin/mtn3 genes. Similarly, Sang et al. [10] found at a broad level, that the proportion and abundance of stigma proteins in different functional categories (e.g. “defense and stress response”, “carbohydrate and energy metabolism”, “protein metabolism and folding”) was similar between maize (dry) and tobacco stigmas (wet), indicating that in general, similar processes occurs in both types of stigmas. However, the specific proteins found in “signal transduction” and “lipid metabolism” categories showed low protein homologies between wet and dry stigmas [10].

Fabaceae is the third largest plant family after Asteraceae and Orchidaceae [11] that is found to be globally distributed but for which few stigma composition studies have been carried out [12,13,14]. Most studies performed so far have been focused on a few species and families such as maize, rice, Lilium, Arabidopsis, Brassica, Crocus, Petunia and Nicotiana. Allen et al. [9] conscious of this reality added a new clade to the pool of studied species: Senecio squalidus. It belongs to the Asteraceae family and possesses a “semidry” stigma, which shows intermediate characteristics between dry and wet stigmas. This condition is characterized by having secretory cells with exudate retained by cuticle or protein pellicle that can be ruptured by pressure or physical friction [15,16]. The stigma of the Fabaceae has been classified as wet or semidry, although some cases of dry stigmas have been reported (e.g. Cassia grandis, Caesalpinia echinata) [17]. The semidry stigma is particularly characteristic of the Papilionoideae subfamily [17], which comprises ~14.000 species [18]. Some of its members are economically and culturally important legume crops such as pea, lentil, chickpea and faba bean. Legumes fix atmospheric nitrogen into available ammonia, promoting nitrogen fertilization of natural soils. Many of them are used for food or forage because of their high content in protein, starch, fibre and other essential nutrients, but they also can be exploited for industrial processes (dyes, gums) or have medicinal properties [19].

In a global climate change context it is expected that the reproductive success of plants, including those involved in agriculture will be affected [20,21]. In addition to physiological alterations caused by abnormal climatic conditions, the reproductive success of entomophilous plants can be also affected by changes in plant-pollinator interactions such as variations in the population distribution of pollinators or the uncoupling of flowering phenology and insect life cycles [22,23,24]. Hence, it is important to extend the knowledge about the mechanisms that promote self-fertilization, since pollinator dependence could restrict plant reproduction under climate change scenarios.

Faba bean (Vicia faba L.) is a partially allogamous species with both cross- and self-fertilization happening in the same plant [25]. Cross-fertilization depends on pollinator activity and unstable yield as well as low fruit and seed sets are related to low visitation rates. On the other hand, self-fertilization occurring by spontaneous selfing could ensure pod and seed set in absence of pollinators [26,27,28]. The ability of a flower to self-fertilize in absence of pollinators or mechanical disturbance is termed as autofertility [29]. The degree of autofertility differ among faba bean genotypes, and it has been related to some floral features like wider style-ovary angle, shorter style, shorter stigmatic papillae, few and shorter stylar hairs, thinner intervening cuticles with rupture previous to anthesis, or lower quantities of pollen grains [27,30,31]. Despite the importance of the rupture of the stigmatic cuticle for successful fertilization in faba bean flowers, little is known about the underlying processes taking place on the stigmas. Recently, a highly saturated genetic map was built and several quantitative trait loci (QTLs) associated with different autofertility traits were detected. Three of the QTLs were related with the rupture of the stigmatic cuticle in chromosome VI, although the function of the associated marker was not clearly related with autofertility [32].

Advances in faba bean breeding have been slow and costly due to its large genome (13Gbp) and its mixed breeding system. RNA-Seq analysis is a relatively inexpensive method and provides data for single-nucleotide variations, clarifying transcriptional and post-transcriptional gene regulation and transcript rearrangements. Differentially expressed genes (DEGs) can be identified with this method to facilitate our in-depth understanding of key biological and physiological mechanisms. Although some comparative transcriptomic analyses have been performed in faba bean using different tissues to understand stress responses such as drought [33], frost [34] or disease resistance [35] no transcriptional information on the genes involved in the fertilization process are still available and the molecular basis of this essential process are still unknown. Herein, we have performed a transcriptome analysis of styles and stigmas of faba bean flowers from lines contrasted for autofertility with two main objectives: (i) amplify the genetic information available for stigmas in a different plant species and family, (ii) identify differentially expressed genes between autofertile and autosterile lines to better understand functional biology underlying this important trait.

2. Results

2.1. Transcriptome sequencing and De novo assembly

A total of 1,189,079,630 raw reads were obtained from the 18 libraries. After quality control and filtering the total number of reads was 1,077,199,910. A summary of the transcriptome de novo assembly data is shown in Table 1. The assembly of sample Vf27.18 in Trinity produced 76,269 contigs with a N50 of 2,387 bp and an average contig length of 982.9 bp.

2.2. Annotation and differential expression analysis

From the 76,269 contigs of the whole transcriptome, 45,1% of the sequences (34,421 contigs) could be annotated with at least one GO term against the PLAZA 4.5 dicots database in TRAPID and 34,379 of them were assigned to 7,720 gene families. Besides 29,1% of the contigs showed full-length or quasi full-length sequences although more than 55% of the transcripts provided no information.

Gene Ontology analyses retrieved a total of 8,439 different GO terms assigned, which were summarized according to GO Slims categories for plants in Figure 1. A total of 115,920, 75,489 and 70,801 annotations were assigned to Biological Process (BP), Cellular Component (CC) and Molecular Function (MF) categories, respectively. The top two within BP were ‘cellular process’ and ‘metabolic process’. Some other numerous terms revealed by the analysis were ‘nucleobase-containing compound metabolic process’, ‘response to stress’, ‘anatomical structure development’, ‘reproduction’, response to different stimuli or ‘transport’. Among the CC category, “intracellular” followed by “cytoplasm” and “membrane” were the most abundant terms. Regarding the MF category, the majority of contigs were annotated within ‘binding’ and ‘catalytic activity’. In the binding category, ‘protein binding’, ‘nucleic acid binding’, ‘nucleotide binding’, ‘DNA binding’ and ‘RNA binding’ are the most numerous categories. On the other hand, ‘hydrolase activity’, ‘transferase activity’, ‘kinase activity’ and ‘transporter activity’ are also important categories (Figure 1).

Differential expression analyses performed in edgeR revealed 5,918 differentially expressed genes (DEGs) between the autofertile lines (AF) and the autosterile lines (AS). Of them, 3,443 genes were upregulated (higher expression values in AF than in AS) and 2,475 genes were downregulated (with significant lower expression values in AF than in AS). The KEGG pathway enrichment analysis using KOBAS-i indicated that the up- and downregulated genes were significantly enriched in 39 functional groups, being “Biosynthesis of secondary metabolites”, “Metabolic pathways” and “Starch and sucrose metabolism” the most significantly enriched terms in both groups. Upregulated genes were enriched in “Selenocompound metabolism”, “One carbon pool by folate”, “Monoterpenoid biosynthesis”, “Nitrogen metabolism” or biosynthesis of certain aminoacids like arginine, valine, leucine and isoleucine. On the other hand, downregulated genes were particularly enriched in “Limonene and pinene degradation”, “Phosphatidylinositol signaling system”, “Inositolphosphate metabolism”, “AGE-RAGE signalling pathway in diabetic complications”, “Phagosome”, “ABC transporters” or “Glycerolipid metabolism” (Figure 2). Some of these significant KEGG terms were exclusive of up- or downregulated genes. Thus, “Nitrogen metabolism” and “Monoterpenoid biosynthesis” were exclusive of upregulated genes, whereas “Limonene and pinene degradation”, “AGE-RAGE signaling pathway in diabetic complications” and “Glycerolipid metabolism” were exclusive of the downregulated genes (Figure 2).

The GO annotation analysis of the DEGs between AF vs. AS performed in TRAPID showed that 2,802 out of 5,918 contigs (47,3%) were annotated with at least one GO term and 2,793 transcripts were assigned to 1,285 gene families. Based on GO slims categories for plants, the general GO term classification of the DEGs showed a similar distribution to the one exhibited by the whole transcriptome (Supplementary file 1). Thus, ‘cellular process’ and ‘metabolic process’ were the most abundant subcategories in BP, ‘binding’ and ‘catalytic activity’ in the MF and ‘intracellular’, ‘cytoplasm’ and ‘membrane’ in the CC categories.

The GO enrichment analysis performed in TRAPID showed that more than 60, 10 and 50 GO terms were significantly enriched in the BP, CC and MF categories, respectively (Supplementary file 2). Among the BP component, the most enriched general terms for the upregulated genes in AF vs. AS were “glycoside metabolic process”, “aminoglycan metabolic process”, “chitin metabolic process”, “glucosamine-containing compound catabolic process”, “cell wall macromolecule catabolic process”, “nucleotide catabolic process”, “salicylic acid catabolic process” or “lignin catabolic process” with log2 enrichment values > 2. In the downregulated genes, GO terms like “negative Rho protein signal transduction”, “negative regulation of Ras protein signal transduction”, “pollen tube adhesion”, “protein homotetramerization” or “phosphatidilinositol-mediated signaling” were the ones showing log2 enrichment values > 2. In the CC category only one GO term (“plant-type cell wall”) was enriched in the upregulated genes (log2 enrichment value of 0.5). However, among the downregulated genes, terms like “exocytic vesicle”, “secretory vesicle”, “apical plasma membrane” or “pollen tube” were revealed (log2 enrichment values > 1.6). Finally, within the MF category GO terms like “(+)-neomenthol dehydrogenase activity”, “(-)-menthol dehydrogenase activity”, “xanthoxin dehydrogenase activity”, “chitinase activity”, “bis (5’-nucleosyl)-tetraphosphatase (assymetrical) activity”, “serine-type endopeptidase inhibitor activity”, “bis (5’-adenosyl)-penthaphosphatase activity”, “phenylalanine ammonia lyase activity”, or “diphosphoric monoester hydrolase activity” were the sub-functional categories found for the upregulated genes (log2 enrichment values > 2). On the other hand the downregulated genes showed enriched GO terms related with “pectate lyase activity”, “Rho GTPase binding”, “phosphatidylinositol kinase activity”, “Rab GTPase binding”, “B-fructofuranosidase activity”, “sucrose alpha-glucosidase activity”, “solute: proton antiporter activity” or “1-phosphatidylinositol binding” (log2 enrichment values > 2).The text continues here.

2.3. Search of DEGs within QTLs intervals for autofertility

In the recent high-density genetic map built for the RIL population Vf6 x Vf27 [32], used as well in this study, the authors reported 31 QTLs for traits related with autofertility. The number of markers within the QTLs intervals ranged from four to 105, and most of them (>70%) matched with a known genome sequence. We checked the differentially expressed genes (DEGs) between AF and AS lines and overlaid these DEGs on to the QTLs confidence intervals to find candidate genes associated with the trait. One DEG (including two transcripts) matched with the significant marker associated with the QTL NORMALQ located in chromosome VI and related with the quantity of pollen with regular size. Besides, up to 18 QTLs had at least one DEG falling within the QTLs confidence intervals (Supplementary file 3). Thus, 31 of the DEGs (transcriptome sequences) matched with at least one QTL marker. The corresponding genome sequences were selected and blasted (BLASTx) to determine the putative function of these genes (Table 2).

Regarding the rupture of the stigmatic cuticle QTL ascribed to chromosome VI, no DEGs matched with any of the markers in the genetic map (see Supplementary file 3). Nevertheless, when considering the 45 genes included in the QTL confidence interval, four DEGs matched with two genomic sequences: Vf6g026360 and Vf6g026920. The DEG that matched with Vf6g026360 is a full length sequence of a SMALL AUXIN UP-REGULATED RNA 54-like protein (AT2G28085 in Arabidopsis thaliana) and was upregulated in autofertile lines. On the other hand, the three transcripts matching with Vf6g026920 correspond to partial sequences of an ABC transporter G family member 28 (AT5G60740), downregulated in autofertile lines (Supplementary file 4).

3. Discussion

3.1. Sequencing, assembly and annotation.

In this study more than 1,000 million of clean reads have been obtained from the transcriptome sequencing of styles and stigmas of faba bean flowers using the Illumina Novaseq platform. A minimum of 42 millions of clean reads per sample and a maximum of 67 million were obtained (with the exception of the sample Vf27.18 sequenced at a higher depth that produced more than 121 millions of clean reads). The de novo assembly resulted in 76,269 contigs, with a N50 of 2,387 bp and an average contig length of 982.9 bp. Compared to other recent faba bean transcriptome studies [36,37,38,39,40], the number of unigenes acquired in our study was intermediate, but the N50 value was higher. All these previous studies used the Illumina platform to sequence the libraries and most of them performed the assembly with the Trinity software.

Concerning the annotation, only 34,421 out of 76,269 sequences (45.1%) could be annotated with at least one GO term, using the TRAPID software. This percentage was relatively low despite the several softwares assayed to analyse the data. The percentage of annotation for the above mentioned faba bean transcriptomes varied depending on the databases used, ranging from 40.8% [39] to 71.5% [38] for unigenes annotated against the NCBI non-redundant (Nr) protein sequence database, which is usually the one to obtain the highest annotation percentages. Poor annotation rates were as well reported in the stigma transcriptomes from different species. Thus, Wang et al. [41] annotated 43.7% of the sequences in jasmine, He et al. [42] annotated 53.8% of the sequences in Camellia oleifera and Quiapim et al. [43] reported no hits or known functions matches in 52.1% of the Nicotiana novel sequences after BLASTx searches. Besides, in a recent proteome and transcriptome analysis using stigmas and pollen from Brassica, Robinson et al. [7] pointed out that many of the proteins revealed in their study still have no known biological roles. This relatively low percentage of annotation or identification may be a result of the scarce genetic information available for flower stigmas and the underlying processes beyond the genetics of incompatibility systems. In addition, many of these sequences may correspond to rapidly evolving species specific genes involved in sexual reproduction which display high diversity in order to maintain species boundaries [5,8].

Previous works on the pistil composition highlight broad similarities between species with wet and dry stigmas for some functional categories such as ‘defense and stress response’, ‘carbohydrate and energy metabolism’, ‘protein metabolism and folding’, ‘cell wall remodeling’, ‘signal transduction’, ‘photosynthesis’ or ‘lipid metabolism’ [10]. Our analysis also showed numerous GO annotations related to these broad categories like ‘response to stress’ or other types of stimulus, ‘carbohydrate metabolic process’, ‘protein metabolic process’, ‘signal transduction’ or ‘lipid metabolic process’ (Figure 1). Similarly the KEGG enrichment analysis also showed pathways related to these categories like: ‘plant-pathogen interaction’, ‘pyruvate metabolism’, ‘biosynthesis of amino acids’, ‘MAPK signaling pathway’, ‘Phosphatidylinositol signaling system’, ‘fatty acid degradation’, etc. ( Figure 2).

One of the goals of this study was to identify genes differentially expressed between autofertile and autosterile lines. The KEGG enrichment analysis revealed several statistically enriched pathways in this set of genes. The upregulated DEGs in the autofertile lines, were mostly enriched in pathways related to the synthesis of amino acids such as: ‘Selenocompound metabolism’, ‘Valine, Leucine and Isoleucine biosynthesis or ‘Arginine biosynthesis’. Another highly enriched pathway was ‘One carbon pool by folate’. Folates act as donors and acceptors in one-carbon transfer reactions and are involved in the synthesis of important biomolecules such as aminoacids, nucleic acids and vitamin B5 (reviewed in [44]). But it has been also related to stress responses, and among them, the response to oxidative stress.

Metabolic pathways related with the synthesis or degradation of terpenoids were also highly enriched in the DEGs. ‘Monoterpenoid biosynthesis’ was enriched in the upregulated genes, whereas ‘Limonene and pinene degradation’ (two monoterpenes) were notably enriched in downregulated genes. Besides, ‘Diterpenoid biosynthesis’ was significantly enriched in both, up- and downregulated genes. Many monoterpenoids are volatile compounds and can be found in the essential oils of many plants. The biological functions of many of them are related to attraction or repellent of insects such as pollinators or herbivores [45]. For example, three volatile monoterpenes (linalool, limonene and β-pinene) can be identified by wasps from receptive female flowers of figs, which is the only stage receptive to pollinators [46]. Arabidopsis thaliana mutant plants that lacked the emission of a volatile sesquiterpene showed greater bacterial growth on their stigmas than did flowers of wild-type plants [47]. On the other hand, it has been seen that the beetle Bruchus rufimanus, an important faba bean pest, responds to floral volatiles in physiological and behavioral experiments, though the beetle did not necessarily pollinate the flowers [48]. Therefore, terpenes play important roles in defense against biotic interactions, as could be the case in faba bean flowers, with the already receptive stigmas emitting volatiles of monoterpenes for both, protecting against pathogens and attracting pollinators. In addition to terpenoid metabolism, the ‘Phenylpropanoid biosynthesis’ was also enriched in styles and stigmas of faba bean flowers. Phenylpropanoids are also part of the secondary metabolism of plants, contributing to all aspects of plant responses to abiotic and biotic stimuli [49].

Regarding signal transduction, several routes stand up in the enrichment analysis, highlighting the importance of recognition of different stimuli in stigmas. The ‘MAPK signaling pathway’ was enriched in both, up- and downregulated genes, whereas the ‘Phosphatidylinositol signaling system’ and ‘AGE/RAGE signaling pathway’ were notably enriched in the downregulated genes. A Mitogen-activated protein kinases (MAPKs) cascade is required for maintaining the stigma receptivity to accept compatible pollen in Arabidopsis. MAPKs converges in the receptivity factor Exo70A1, a member of the exocyst complex. Phosphorylation of Exo70A1 by MAPKs regulates pollen hydration and germination through exocytosis in Brassica and Arabidopsis species [50]. As reported by McInnis et al. [51], the accumulation of reactive oxygen species (ROS) in mature stigmas in a constitutive way suggests that ROS might be an upstream candidate signal as they are known to activate these kinases. Regarding the AGE-RAGE signaling pathway, it is better known in animals than plants. AGEs is the acronym of advanced glycation end products. AGEs are involved in pathogenesis of diabetes mellitus, Alzheimer’s disease, aging and are also involved in thermal processing of foods. Multiple membrane and soluble proteins have been annotated as receptors for glycation products in mammals (e.g. RAGEs). Upon interaction with receptors, AGEs trigger inflammatory response by activation of mitogen-activated protein kinase (MAPK-), janus kinases (JAC-) and mitogen-activated protein kinases/extracellular signal-regulated kinases MAPK/ERK -signaling pathways [52]. However the role of glycation in plants is poorly understood, and two main aspects are proposed: glycation as a marker of aging, senescence and tag for protein degradation and as a possible mechanism of signaling (reviewed in [53]).

On the other hand, the ‘Phosphatidylinositol signaling pathway’ was enriched in downregulated DEGs. Inositol phospholipid compounds (such as IP3 and DAG) on the cell membrane are important secondary messengers involved in signal transduction [54]. For example, many components of the phosphatidylinositol signaling system participate in the vacuolar diversification during pollen development and vesicle transport in pollen tube growth. A good regulation of phosphatidylinositol-4-phosphate and phosphatidylinositol 4,5-bisphosphate pools are necessary for polarized secretion in plants (reviewed in [55]). Gradients of these compounds have been observed in root hairs and pollen tubes where they are linked to polarized secretion [56,57]. Since the stigmatic papillae in faba bean are specialized in secreting the exudates, functions related with vesicle transport are expected to be found in this tissue.

3.2. DEGs within QTL intervals previously described for autofertility

From the significant markers found to be associated with autofertility traits by QTL analysis [32], two transcripts matched, in chromosome VI, with the significant QTL marker for NORMALQ, related with the quantity of pollen with a regular size. This QTL was associated with a LARP6c protein which is exclusively expressed in pollen, being necessary from male gametogenesis to the end of the fertilization process. LARP6c participates through spatio-temporal control of translation to male fertilization [58]. The detection of pollen genes in our stigma/style transcriptome data was unsurprising considering that the flower sampling was performed previous the anthesis and some pollen grains might have been already attached to the stigma.

Of particular interest in this study were the genes involved in the receptivity of the stigmas and the autofertility. Successful pollination, fertilization and seed set depend upon receptivity of stigmas during the few days following anthesis. Several DEGs identified in the transcriptome co-mapped in the intervals of QTLs related with the rupture of the stigmatic cuticle [32]: %RUPT (percentage of rupture of stigma cuticle), %RUPTAREA (percentage of ruptured area) and RUPTL (rupture length) (see [31] for further details about these measures), which co-localized as well in chromosome VI. For these QTLs, four DEGs matched with two genomic sequences within the QTLs intervals. The genomic sequence Vf6g026360 was identified as a SMALL AUXIN UP-REGULATED RNA 54-like protein (SAUR-like proteins). SAURs are one of the largest early auxine-responsive gene families and specially in Fabaceae lineages seems to be one of the largest (e.g. 141 and 236 different SAUR proteins found for Medicago truncatula and Glycine max, respectively). SAUR proteins play important roles in developmental regulation and environmental responses and different clades in the gene family show different expression profiles and functions [59]. In M. truncatula, it was found that some SAUR genes showed evident tissue-specific expression patterns. For example, MtSAUR18, MtSAUR57 and MtSAUR91 displayed unusually high expression levels in flowers, suggesting that they may be involved in flower development, whereas MtSAUR13 exhibited root meristem and root tip-specific expression patterns [60]. Auxin related proteins were also preferentially expressed in Nicotiana tabacum stigmas [43], corroborating the importance of this hormone in pistil development.

The second genomic sequence associated with a DEGs was Vf6g026920, identified as an ATP-binding cassette (ABC) transporter G family member 28. ABC transporters in plants are more numerous than in other organisms, and are classified into eight subfamilies: A-G and I. They are composed of nucleotide-binding domains (highly conserved) and transmembrane domains, being the latter very variable, allowing the transport of different substrates. Full size ABC proteins can work like transporters themselves, but half-size transporters can form complexes to perform their functions. Many full size ABCG transporters are implicated in defense against biotic stresses (e.g. see [61]). Two half-size ABCG transporters of M. truncatula that are present in peri-arbuscular membranes are implicated in arbuscule development in mycorrhizal symbiosis [62]. Another two half-size ABCG transporters are implicated in the stigma exertion in Medicago [63]. AtABCG28 is a critical half-size transporter of A. thaliana that establishes the correct level of reactive oxygen species (ROS) at the pollen tube and root tips. AtABCG28, is specifically localized to the membranes of secretory vesicles and expressed in mature pollen and growing pollen tubes. It is involved in sequestering polyamines (source of ROS) into the vesicles that move and fuse to the growing tip [64]. Since these QTLs are implicated in the rupture of the stigmatic cuticle, which is also related with the presence of exudate and receptivity of the stigma, high levels of ROS are expected in this tissue. Therefore regulation of ROS levels and transport of these substances are important to maintain the correct cellular functions and prevent cell damage.

4. Materials and Methods

4.1. Plant materials and sample collection

The recombinant inbred line (RIL) faba bean population of 124 individuals derived from the cross between lines Vf6 and Vf27 has been previously used for the localization of QTLs related to autofertility, dehiscence, flowering time and other yield related traits [31,65,66]. The parental line Vf6 is a highly autosterile and asynaptic line whereas Vf27 is considered highly autofertile. The materials selected in this study were six genotypes from this RIL population: the two parental lines (Vf6 and Vf27), two highly autosterile RILs (RIL19 and RIL96) and two highly autofertile RILs (RIL14 and RIL44).

Plants were grown in 5L pots under controlled conditions (22 ºC, 14 hours day- 10 hours dark). At the peak of flowering production for each line, flowers previous to anthesis were collected over several days and dissected to extract the flower style. Style samples were immediately frozen in liquid N₂ and stored at -80ºC until the RNA extraction was performed.

4.2. RNA extraction, sequencing and De novo assembly

Total RNA of approximately 100 styles per sample was extracted using TRIZOL reagent (St. Louis, Missouri, USA) with the Zymo Research Direct-zolTM RNA MiniPrep Kit according to the manufacturer’s instructions. A total of 18 samples were finally prepared consisting of three replicates for each of the six genotypes.

Samples were sent to STABVIDA (Caparica, Portugal) for quality control, library construction (with a Stranded mRNA Library Preparation Kit), sequencing (Illumina Novaseq, 150bp paired-end reads) and assembly. Raw sequences were trimmed to generate high quality reads. For each original read, the following parameters were applied: a quality trimming based on a quality score of 0.01 (error probability), a limit of the length of ambiguity of 2 nt and a minimum read length of 30 nt. Sample Vf27.18 was sequenced at a higher depth and the high quality sequence reads were used for the de novo assembly in Trinity 2.8.4 [67]. The assembled transcriptome of sample Vf27.18 was used as the reference sequence for the expression analysis. Raw reads of this study have been deposited into the NCBI Sequence Read Archive (SRA) database under the accession number PRJNA1044928.

4.3. Annotation and differential expression analysis

Contigs obtained from the assembly of sample Vf27.18 were annotated with TRAPID [68], a web application for taxonomic and functional analysis, using PLAZA 4.5 dicots [69] as database and clade Papilionoideae as similarity search database with a threshold of 10-5. GO graphs were summarized according to GO slims categories for plants.

The high quality reads from each sample were mapped against the de novo assembled transcriptome reference. A minimum similarity and length fractions of 0.8 were used as parameters to consider a correct mapped read. The differential expression analysis was performed with edgeR package [70] in R v. 4.2.1 [71]. The identified differentially expressed genes (DEGs) were filtered using a fold change value of > 2 or < -2 and a FDR (False Discovery Rate) p-value < 0.05 as thresholds.

In order to identify significant metabolic pathways correlated with putative autofertility genes, we focused on the subset of DEGs existing between all the autosterile and all the autofertile samples. A KEGG pathway enrichment analysis of the DEGs was performed using the KEGG pathway database in KOBAS-i, the more recent KEGG Orthology Based Annotation System [72] using Medicago truncatula as the reference database. A p-value < 0.05 was considered to indicate significant over-representation of a certain KEGG pathway. We also performed a Gene Ontology (GO) enrichment analysis in TRAPID, which determines the over-representation of a certain GO term, compared to the background frequency (i.e. Papilionoideae dataset). The Benjamini & Hochberg correction was further applied to control multiple testing and decrease the FDR (q-value <0.05 was established as threshold to determine if the GO term was enriched in the dataset).

4.4. Search of DEGs within QTL intervals previously described for autofertility

To identify candidate genes controlling autofertility we combined the results of a previous QTL mapping [32] with the transcriptome data. Molecular markers falling within the QTLs intervals were selected and the corresponding nucleotide sequences were extracted from ‘Vfaba_v2’ 60k SNP Array [73,74]. Marker sequences were first aligned against the faba bean genome [75] and the genome sequences were then aligned against the transcriptome sequences of the DEGs. Those DEGs falling within the QTL intervals were identified by BLASTx. Regarding the QTL related with the rupture of the stigmatic cuticle in chromosome VI, we selected all the genomic sequences within the QTL intervals (not only the ones corresponding to an associated marker in the genetic map). The 45 sequences included in the interval between Vf6g024800 and Vf6g123480 were aligned against the DEGs from the transcriptome and were further identified by BLASTx.

5. Conclusions

In this study we used RNA-sequencing to check for the first time the differential expression of gene transcripts between faba bean lines differing in autofertility. DEGs were overlaid onto QTLs detected in a recent high-density genetic map to find candidate genes associated with autofertility. Initial challenge in the current study was due to the lack of annotated stigma datasets. DEGs up- and downregulated were identified and some of them predicted to be related with the traits under study. One DEG (including two transcripts) matched with the significant marker associated with one QTL related with the quantity of pollen with a regular size and others DEGs co-mapped in the intervals of QTLs related with the rupture of the stigmatic cuticle. RNAseq combined with QTL mapping are a powerful approach for identifying candidate genes and the results derived from this work provides an important transcriptomic reference for style-stigma processes to aid understanding of the molecular mechanisms involved in faba bean fertilization. The new available transcriptomic data and the RIL population used will facilitate the fine mapping of the responsible genes and will provide targets for future study and improvement of autofertility traits in this crop.