All plants used in this study were in the Columbia (Col-0) ecotype. For all analyses, Arabidopsis thaliana seeds were incubated for 2 days at 4°C after sowing. Plants were grown at 22°C under long day conditions (LD). The ces-4 mutant was obtained from NASC (SAIL_674_A01 lines). This line phenocopy the ces-3 (haf), ces-1 and ces-2 mutant alleles that have been previously characterized (Crawford and Yanofsky, 2011; Eremina et al., 2016). The stk, shp1, and shp2 mutants were provided by M. Yanofsky (Pinyopich et al., 2003), the bee1 bee2 bee3 triple mutant was provided by J. Chory (Friedrichsen et al., 2002). The other mutant combinations analyzed were obtained by crossing. The pCES::GUS line was provided by Prof. B. Poppenberger (Poppenberger et al., 2011).

Genotyping of the mutants was done by PCR analysis. All the primers used in this study are listed in Supplementary Table S11.

Arabidopsis flowers and siliques were fixed and embedded in paraffin as previously described in Dreni et al. (2011). Tissue sections (8 μm) were hybridized with a CES specific antisense RNA probe and observed with a Zeiss Axiophot D1 microscope. Images were captured on an Axiocam MRc5 camera (Zeiss) using the Axiovision program (version 4.1).

GUS staining of pCES::GUS was performed as described in Balanzà et al. (2016).

The percentage of unfertilized ovules per fruit was determined by opening at least 20 mature siliques for each genotype, counting the number of unfertilized ovules and the total number of developed seeds per silique under a stereomicroscope. The statistical significance of the differences in unfertilized ovules was analyzed using an Anova test followed by the Tukey HSD test (^∗∗p < 0.01). All the experiments were repeated four times. For each experiment 5 siliques were collected from 4 independent plant lines of the same genotype.

Biological samples were collected and fixed overnight at 4°C in FAA solution (3.7% formaldehyde, 5% acetic acid and 50% ethanol). Fixed tissues were washed with water and post-fixed with aqueous 2% osmium tetroxide for 2 h at room temperature. Tissues were rinsed several times in deionized water and dehydrated in a graded series of ethanol for 15 min per rinse. This step was followed by critical point drying with liquid CO₂ and sputter-coating with gold in a Nanotech sputter coater. Specimens were analyzed using a LEO 1430 Scanning Electron Microscope.

Experiments of pollen tube guidance were performed as previously described by Mendes et al. (2016) and images were captured on an Axiocam MRc5 camera (Zeiss) using the Axiovision program (version 4.1).

Coding Sequences (CDS) of CES and STK were amplified by PCR with the primers indicated in Supplementary Table S11 and first cloned into pDONOR207 (Life Technologies) and subsequently into the pYFPN43 and pYFPC43 vectors.¹ BiFC experiments were conducted as previously described in Balanzà et al. (2016). Essentially, the abaxial surfaces of tobacco leaves were observed 2–5 days after agroinfiltration and the interactions were monitored using a Leica TCS confocal microscope.

The two-hybrid assays were performed at 22°C in the yeast strain AH109 (Clontech), using the cotransformation technique (Egea-Cortines et al., 1999). The coding sequences of CESTA, STK, and BEE1 were cloned in the Gateway vector GAL4 system (pGADT7 and pGBKT7; Clontech) passing through pDONOR207 (Life Technologies). Yeast two-hybrid interaction assays were performed on selective yeast synthetic dropout medium lacking Leu, Trp, Adenine, and His supplemented with different concentrations of 3-aminotriazole (1, 2.5, 5, 10, and 15 mM of 3-AT).

For the qRT-PCR experiment of ces-4 (knockout validation) total RNA was extracted from whole inflorescences using the LiCl method as previously described (Verwoerd et al., 1989). First-strand cDNA was synthesized using an IMProm-II ^TM Reverse Transcription System (Promega). Enrichments fold were detected using a SYBER Green assay (Bio-Rad).² The qRT-PCR assay was performed in triplicate using a Bio-Rad iCycler iQ optical system (software version 3.0a). The qRT-PCR assay was repeated three times. For each experiment the RNA was extracted from three independent plant lines.

Total RNA was extracted from three biological replicates (1 gr, obtained from a pool of 5 plants for each replicate) from both wild-type and bee1 bee3 stk ces-4 mutant inflorescences till stage 12, using the Macherey Nagel “Nucleospin RNA Plant” according to the manufacturer’s instructions. RNA concentrations and integrity were determined using Qubit Fluorometer and the Qubit^TM RNA XR Assay Kit (Thermo Fisher Scientific). Sequencing libraries were prepared using the NEBNext Ultra II Directional RNA library Prep Kit for Illumina (NEB) according to the manufacture’s instruction and sequenced on the HiSeq Illumina platform. Reads were mapped on the reference Arabidopsis thaliana genome (TAIR, version 10) using the bowtie2 program (Langmead and Steven, 2012; Supplementary Table S4). Estimation of gene expression levels was performed using RSEM (Li and Dewey, 2011). Identification of differentially expressed genes was performed by the quasi-likelihood F-test as implemented by edgeR (Robinson et al., 2009). A False Discovery Rate (FDR) cut-off value of 0.05 was applied for the identification of significantly differentially expressed genes. Graphical representation of the data was performed by means of the gplots R library. To validate the data obtained from the RNAseq experiment, we extracted total RNA from wild type and bee1 bee3 stk ces-4 mutant inflorescences till stage 12, using the Macherey Nagel “Nucleospin RNA Plant” according to the manufacturer’s instructions. The cDNA and qRT-PCR were performed as described in the previous paragraph “RNA extraction and qRT-PCR analysis.” The qRT-PCR assay to validate RNAseq experiment was repeated three times. For each experiment the RNA was extracted from three independent plant lines that were used for RT-PCR analysis.

Enrichment analysis of Transcription Factor Binding Sites (TFBS) was performed by means of the Pscan tool, using the “non-redundant” set of Jaspar TFBS matrices (Fornes et al., 2019) and considering a region spanning from the TSS to 500 bp upstream (default settings). Archetype consensus sites for the W-box, MADS-box and G-box elements were established as the Jaspar TFBS matrix showing the highest enrichment in our dataset according to Pscan. To identify genes with high confidence W-box, MADS-box and G-box elements in their promoters, the complete set of 33239 Arabidopsis transcripts according to the TAIR10 annotation were scanned using the highly significantly enriched TFBS matrices as defined above. Distribution of scores were established for each element, and an Expectation-Maximization algorithm as implemented in the R EM-cluster package was used to fit a 2 Gaussian mixture model. Based on the EM model promoters were partitioned in “high” and “low” scoring. Only promoters that were confidently assigned to high scoring clusters (confidence > 95%) were considered to have a bona fide TFBS element. Similarly, to our previous analyses promoters were defined as regions spanning 500 bp upstream of the TSS. P-values for the enrichment of genes with W-box, MADS-box and G-box elements in their promoters were performed by means of the hypergeometric distribution test as implemented in the R stats package.

The GenBank/EMBL accession numbers for the genes shown in this study are the following: BEE1, At1g18400; BEE2, At4g36540; BEE3, At1g73830; CES, At1g25330; SHP1, At3g58780; SHP2, At2g42830; STK, At4g09960. The accession number of the genes identified in the transcriptomic analysis are indicated in Supplementary Tables S5, S6.

In order to identify new targets, partners and/or regulators of STK, a bioinformatic analysis was previously performed in the laboratory to discover genes co-expressed with STK (Mendes et al., 2016). This kind of analysis assumes that correlation between the patterns of gene expression in a large range of different experimental conditions could indicate a functional relationship (Aoki et al., 2007). As already described, the procedure to identify genes that have expression profiles that correlate with STK was proven successful since among the top listed genes SHP1 and SHP2 were found, which act redundantly with STK in the control of ovule identity (Pinyopich et al., 2003). Furthermore, top rated were VERDANDI (VDD) and VALKYRIE (VAL), targets of the STK-SEPALLATA3 (SEP3) complex, which both act in the control of the death of the synergid cell that is receptive of the pollen tube (Mendes et al., 2016). CES (At1g25330) rated in position 14 of the correlators, having an expression correlation with STK of p(LIN) 0.759 and p(LOG) 0.505 (Mendes et al., 2016).

Previously, it has been reported that CES is expressed from stage 8 to 15 of flower development, first in CMM and later in the stigma, septum and funiculi (Crawford and Yanofsky, 2011; Poppenberger et al., 2011). Since STK is expressed during female reproductive organ development (mainly septum, transmitting tract, replum, funiculus and ovule/seed) (Pinyopich et al., 2003; Mizzotti et al., 2014; Herrera-Ubaldo et al., 2019; Di Marzo et al., 2020), we investigated in more detail the spatio-temporal expression pattern of CES during pistil, ovule and seed development. We performed both RNA in situ hybridization and the analysis of the pCES::GUS reporter gene construct (Figure 1). The RNA in situ analysis experiments showed that CES was expressed in ovules during early stages of development until mature ovule formation and its transcripts were also present in the transmitting tract before fertilization (Figures 1A–C). In mature ovules, a high expression was detected in the funiculus (Figure 1C). After fertilization and during embryo development CES mRNA was present at high levels in the early globular stage, and at heart and torpedo stages mRNA levels were lower (Figures 1D–F). Analysis of the pCES::GUS plants confirmed the in situ data (Figures 1G–I) and these experiments showed that STK and CES partially overlap in their expression patterns during ovule and transmitting tract development (Pinyopich et al., 2003; Crawford and Yanofsky, 2011; Poppenberger et al., 2011; Herrera-Ubaldo et al., 2019).

To further analyze the role of CES, we characterized a CES loss-of-function allele that we obtained from the SAIL collection (Alonso and Stepanova, 2003). This allele, that we called ces-4, has a T-DNA insertion in the second exon, differently from ces-1, ces-2, and ces-3 (haf), previously published T-DNA knockout insertion alleles (Crawford and Yanofsky, 2011; Poppenberger et al., 2011; Eremina et al., 2016). The qRT-PCR analysis performed to check CES expression in the ces-4 mutant indicated that this is also a null allele (Supplementary Figure S1). Analysis of ces-4 homozygous plants did not reveal any general phenotype at the level of the entire plant, which was consistent with previous reports for the ces (haf) mutant alleles (Crawford and Yanofsky, 2011; Poppenberger et al., 2011; Eremina et al., 2016). We crossed the ces-4 mutant with stk to obtain the stk ces-4 double mutant. The double mutant was indistinguishable from wild type plants apart from a large number of unfertilized ovules at the lower part of the siliques. In wild type plants on average 1.0% of the ovules in a silique remained unfertilized. In the stk mutant we observed on average 6.2% unfertilized ovules, and in the ces-4 single mutant 11.6% (^∗∗p < 0.01 in respect to wild type; Figure 2 and Supplementary Table S1). In the stk ces-4 double mutant this phenotype was further increased to 28.1% of unfertilized ovules (Figure 2 and Supplementary Table S1). Like in the ces-4 single mutant, also in the stk ces-4 double mutant the unfertilized ovules were located always at the lower part of the siliques (Figure 2 and Supplementary Table S1).

Since SHP1 and SHP2 are closely related to STK and have shown to control together with STK ovule identity redundantly (Pinyopich et al., 2003), we investigated whether they also play a role together with CES in the control of ovule fertilization. Therefore, the shp1 shp2 double mutant was crossed with the ces-4 mutant to obtain the shp1 shp2 ces-4 triple mutant. The shp1 shp2 double mutant showed on average 7.0% unfertilized ovules (Figure 2 and Supplementary Table S1). The shp1 shp2 ces-4 triple mutant showed a clear increase to 16.2%, also in this case the unfertilized ovules were localized at the lower part of the siliques (^∗∗p < 0.01 in respect to wild type and shp1 shp2; Figure 2 and Supplementary Table S1). Unfortunately we could not analyze the stk shp1 shp2 triple mutant, because the triple mutant does not develop normal ovules (Pinyopich et al., 2003). These results suggest that the MADS-domain transcription factor STK, SHP1 and SHP2 have a role together with CES in the control of ovule fertilization.

The contribution of the BEE genes to this process was also investigated (Figure 3 and Supplementary Table S2). CES is phylogenetically related to BEE1 and BEE3, which have been described as redundant with BEE2 (Friedrichsen et al., 2002). Moreover, it has been shown that CES together with BEE1 and BEE3 redundantly control the development of the transmitting tract in the pistil (Crawford and Yanofsky, 2011; Poppenberger et al., 2011). The single mutants bee1 and bee3 and the multiple mutant combinations bee1 bee3 and bee1 bee2 bee3 did not show any obvious phenotype (Figure 3 and Supplementary Table S2). To investigate redundancy between BEE genes and STK, we obtained the bee1 bee2 bee3 stk quadruple mutant, that presented 9.2% of unfertilized ovules, with ^∗∗p < 0.01 when compared to wild type and the triple mutant bee1 bee2 bee3, which again were preferentially positioned at the lower part of the siliques (Figure 3 and Supplementary Table S2). This percentage is lower than the bee1 bee2 bee3 ces-4 quadruple mutant, which presented in average 28.7% of unfertilized ovules (Figure 3 and Supplementary Table S2). Despite the weaker phenotype, our results demonstrated that STK and BEEs also have a role together in ovule fertilization. To investigate if combining the ces-4 mutant could further enhance the unfertilized ovule phenotype, we generated the bee1 bee3 stk ces-4 quadruple and the bee1 bee2 bee3 stk ces-4 quintuple mutants. The results of this analysis showed that bee1 bee3 stk ces-4 had on average 46.6% unfertilized ovules, and a similar result was obtained in the bee1 bee2 bee3 stk ces-4 (44.0%) (Figure 3 and Supplementary Table S2), showing that the addition of the bee2 allele did not further enhance the phenotype. The unfertilized ovules in these cases are localized from the middle to the bottom part of the siliques. Altogether, these results indicated that the MADS-domain transcription factors function not only together with CES, but also with the BEEs in the control of the fertilization of all the ovules within the Arabidopsis fruit.

To better understand the defect that caused the presence of a high percentage of unfertilized ovules at the lower part of the siliques in the different mutant combinations, we performed reciprocal crosses between wild type and stk ces-4 or shp1 shp2 ces-4 mutant plants, to analyze if the female, male or both parts were defective. The outcome of the reciprocal crosses shown in Supplementary Figure S2, suggested that the unfertilized-ovules phenotype was only present when the mutant genotypes were used as mothers. This result indicates that the defects resides only in the female reproductive organ.

To further investigate the defect present in the female reproductive organs, which caused the high number of unfertilized ovules at the lower part of the siliques in the multiple mutant combinations, we analyzed the septum morphology by Scanning Electron Microscope (SEM), at stage 17-B of fruit development, when the siliques reach their maximum size (Roeder and Yanofsky, 2006).

The stk and ces-4 single mutants presented a septum similar to wild type, while the double mutant stk ces-4 disclosed a defect in septum morphology (Figure 4A). The double mutant was characterized by a partial failure of septum fusion in the central portion of the septum in conjunction with the medium portion of the ovary (Figure 4A). The quadruple mutant combinations bee1 bee2 bee3 stk and bee1 bee2 bee3 ces-4 did not display defects in septum morphology when compared to wild type (Figure 4A). However, the damage of the septum appears stronger in the quadruple mutant bee1 bee3 stk ces-4 and in the quintuple mutant bee1 bee2 bee3 stk ces-4 where a hole in the central portion of the septum was seen in conjunction with the medium region of the ovary (Figure 4A).

The septum morphological analysis at stage 17-B suggests a transmitting tract defect caused by a failure of the carpel fusion which happens from stage 9 to 12 (Roeder and Yanofsky, 2006). At these stages, fusion occurs within the two crests of the CMM leading to septum and transmitting tract formation (Bowman et al., 1999; Roeder and Yanofsky, 2006; Reyes-Olalde et al., 2013).

To confirm that the defects found in the septum and transmitting tract formation, described for stk ces-4, bee1 bee3 stk ces-4, and bee1 bee2 bee3 stk ces-4 caused a defective growth of the pollen tubes through the ovary transmitting tract, pistils were hand pollinated with wild type pollen and an aniline blue staining was performed after 36 h (Figure 4B). In the wild-type-hand-pollinated control pistils all the pollen tubes reached the bottom part of the ovary. However, pollen tube growth was compromised in the three mutant combinations. In the stk ces-4 double mutant the pollen tubes did not grow further than the middle region of ovary length (Figure 4B). A stronger phenotype was observed in bee1 bee3 stk ces-4 quadruple and in the bee1 bee2 bee3 stk ces-4 quintuple mutants, where the pollen tubes stopped growing at the upper part of the ovary (Figure 4B).

The differences in septum morphology and in pollen tube growth in the double mutant stk ces-4, in the quadruple mutant bee1 bee3 stk ces-4 and in the quintuple mutant bee1 bee2 bee3 stk ces-4 clearly indicate that the fusion of the carpels, that permits the correct development of the septum and consequently of the transmitting tract, was partially disturbed in these mutants and the cause that the ovules at the lower part of the siliques remained unfertilized.

The experiments described above showed a genetic interaction between STK and CES. To investigate if there is also a physical interaction between these proteins, Bimolecular Fluorescence Complementation (BiFC) was performed. As a positive control, cells were co-transformed with CES (YC-CES) and (YN-CES), as CES was shown to homodimerize in yeast (Poppenberger et al., 2011; Figure 5A). Reconstitution of YFP fluorescence in the nuclei of co-transformed cells confirmed that CES is able to interact with STK in vivo (Figure 5A), supporting our idea that CES and STK form a complex that controls transmitting tract formation. The interaction between STK and CES was also confirmed by Y2H experiments (Figure 5B). As positive control in these experiments we used the previously published interaction between BEE1 and CES (Poppenberger et al., 2011; Figure 5B). Furthermore, we observed that STK and BEE1 were not able to interact. The complete list of vector combinations, interactions and growth media used in this Y2H assay are described in Supplementary Table S3.

To obtain a snapshot of the regulatory network controlled by STK, CES, and BEE proteins, in the development of the septum/transmitting tract, transcriptional profiles of the quadruple mutant bee1 bee3 stk ces-4, which displayed the most severe phenotype, were studied by RNA sequencing. RNA was extracted from the quadruple mutant and wild type inflorescences until stage 12 of pistil development, before fertilization since the transmitting tract is fully developed at this stage (Bowman et al., 1999; Roeder and Yanofsky, 2006; Reyes-Olalde et al., 2013).

For the identification of significantly differentially expressed genes, a FDR cut-off value of 0.05 was applied and a total of 202 differentially expressed genes were identified. A total of 31 genes were upregulated in the mutant (Supplementary Table S5) while 171 genes were downregulated (Supplementary Table S6). Furthermore, functional enrichment analyses revealed that gene ontology categories such as glycosinolate, glycoside and carbohydrate biosynthetic and metabolic processes were enriched for the upregulated genes (Supplementary Table S7), while the downregulated genes were involved in the cell death, cell wall, plant-type cell wall biosynthesis components and response to hormone stimulus (Supplementary Table S8).

Remarkably, among the downregulated genes, three genes were previously described as direct targets of STK (Herrera-Ubaldo et al., 2019): AT3G26140, AT1G28710, and AT1G06080. The AT3G26140 gene which encodes a glycosyl hydrolase (GH5_11) is expressed in CMM during pistil development (Aspeborg et al., 2012; Herrera-Ubaldo et al., 2019). The GH5 enzymes are all mannan endo-beta-1,4-mannosidases which are known to be involved in cell wall remodeling processes (Aspeborg et al., 2012; Lombard et al., 2014; Herrera-Ubaldo et al., 2019). AT1G28710 encodes for a nucleotide diphospho-sugar transferase protein which is involving in the biosynthesis of polysaccharides, that are fundamental components of the ECM (Lennon et al., 1998; Tung et al., 2005; Herrera-Ubaldo and de Folter, 2018). AT1G06080 [DELTA 9 DESATURASE1 (ADS1)], encodes for a fatty acid desaturase which is involved in lipid biosynthesis, expressed in flowers (Fukuchi-Mizutani et al., 1998; Yao et al., 2003; Heilmann et al., 2004; Table 1, see also Figure 6 for heatmap description). Downregulation of these three genes was also confirmed by qRT-PCR (Supplementary Figure S3).

PCD is a key process in the transmitting tract for ovule fertilization (Crawford et al., 2007). Interestingly, several genes involved in the regulation of cell death were downregulated in the quadruple mutant (Table 1 and Figure 6). One of the most important is GRIM REAPER (GRI, AT1G53130) that encodes for a small protein with an extracellular localization that is able to positively regulate cell death (Wrzaczek et al., 2009, 2015; Table 1 and Figure 6). Moreover, we observe that ACCELERATED CELL DEATH 6 (ACD6, AT4G14400) another positive regulator of cell death is also downregulated (Rate et al., 2007; Table 1 and Figure 6).

Additionally, we notice the downregulation of three genes encoding for the Arabinogalactan proteins (AGPs) ARABINOGALACTAN PROTEIN 5 (AGP5, AT1G35230), AGP7 (AT5G65390), and AGP14 (AT5G56540). AGPs are structurally complex plasma membrane and cell wall proteoglycans that are implicated in diverse developmental processes, including plant sexual reproduction (Lennon et al., 1998; Coimbra et al., 2007, 2009). Importantly, these protein are able to regulate cell death in Arabidopsis cell culture (Gao and Showalter, 1999; Guan and Nothnagel, 2004). Intriguingly, we also detect a significant downregulation of ATWSCP (AT1G72290), which encodes for a Kunitz-type protease inhibitor, member of a small family of proteins, that are water-soluble chlorophyll proteins, and LIPID TRANSFER PROTEIN 2 (LTP2, AT2G38530) (Table 1 and Figure 6). ATWSCP is expressed in the septum and transmitting tract during its development, and after fertilization the protein is located in the medial domain of the silique where is able to regulate cell death and pollen tube growth (Bektas et al., 2012; Boex-Fontvieille et al., 2015). The LTP2 gene encodes for a lipid transfer protein necessary for the movement of hydrophobic wax components through the hydrophilic cell wall matrix to reach the cuticle (Kunst and Samuels, 2003). LTP2 mediates cell wall loosening in vitro (Nieuwland et al., 2005).

Numerous genes involved in cell wall biosynthesis processes are mis-regulated in the quadruple mutant, these include BETA-XYLOSIDASE 3 (BXL3, AT5G09730), which encodes for an enzyme involved in arabinan and xylan catabolic process (Minic et al., 2006), and PECTIN METHYLESTERASE 58 (PME58, AT5G49180), encoding for a pectin methylesterases that contributes to cell wall modification (Turbant et al., 2016; Table 1 and Figure 6); both genes were downregulated. Among the upregulated genes, there are three genes of the GDSL family of serine esterases/lipases, AT1G75890, AT4G28780, and AT4G16230 (Table 1 and Figure 6). These genes are involved in numerous aspect of plant development, but in particular the acetylxylan esterase activity of these enzymes could contribute in the degradation of cell walls in association with xylanases, cellulases, and mannanases (Akoh et al., 2004).

Moreover, some genes involved in the auxin signaling pathway were downregulated. These include genes encoding for Aux/IAAs proteins, INDOLE-3-ACETIC ACID INDUCIBLE 1 (IAA1) (AT4G14560), IAA2 (AT3G23030), IAA14 (AT4G14550), and IAA17 (AT1G04250) (Table 1 and Figure 6). The Aux/IAA proteins are the early auxin response proteins and participate in auxin signaling through interacting, as repressors, with AUXIN RESPONSE FACTOR (ARF) proteins (Lavy and Estelle, 2016). In addition, the SMALL AUXIN UPREGULATED RNA 9 (SAUR9) and SAUR12 were also downregulated. The SAURs are small transcripts induced by auxin, that are localized in elongated tissues (McClure and Guilfoyle, 1987; Knauss et al., 2003). SAUR9 and SAUR12 are part of the SAUR10-clade (van Mourik et al., 2017).

Finally, we decided to explore the enriched binding sites from our RNAseq data of the quadruple mutant in order to find possible direct target genes of the complex formed by STK, CES, BEE1, and BEE3. The binding sites used by the MADS-domain and bHLH transcription factors have been characterized: MADS factors bind to CArG boxes (Egea-Cortines et al., 1999; de Folter and Angenent, 2006) while bHLHs recognize E-boxes, being G-boxes the most common form (Quail, 2000; Toledo-Ortiz et al., 2003), in order to regulate the activity of their target genes. On the other hand, W-box is the binding site used by the WRKY transcription factors (Ciolkowski et al., 2008).

Transcription factor binding sites enrichment analyses, performed by means of the Pscan tool (Zambelli et al., 2009) suggest a strong enrichment of G-box and G-box like elements in the upregulated genes (Supplementary Table S9), while a significant enrichment of W-box and MADS-box elements is observed in the downregulated genes (Supplementary Table S10). Notably, according to our analyses G-box elements are identified in 35% (p = 4E-05) of the upregulated genes, while W-box elements and MADS box elements are identified in 17 and 19 downregulated genes (p-values 9.3E-04 and 0.035), respectively. Unsurprisingly, considering the limited number of genes that are differentially expressed in our experimental settings, we observe that only one gene has both a G-box and a CArG-box element in its promoter (AT3G49270, upregulated), only one gene has both a W-box and a MADS-box element (AT1G05450, downregulated), and none of the differentially expressed genes have a W-box and a G-box.

In Arabidopsis, the gynoecium is formed by the fusion of two carpels, which fuse vertically at their margins and give rise to the CMM. From this meristem the transmitting tract develops, an essential tissue to facilitate the growing pollen tubes to reach the ovules and fertilize them.

Recently, a role for STK in relation with NTT was shown during the development of transmitting tract (Herrera-Ubaldo et al., 2019). Here we reveal that STK may also act in the downstream pathway of NTT in a protein complex with CES. Our genetic analysis evidence the role and redundancies between the MADS-box genes STK, SHP1, and SHP2, the bHLH transcription factor gene CES and the closely related bHLH transcription factors BEE1 and BEE3 in the development of the transmitting tract in Arabidopsis thaliana. The CES gene was selected based on high correlation of its expression pattern with STK (Mendes et al., 2016) but also because this gene was shown to play a role in transmitting tract development (Crawford and Yanofsky, 2011; Poppenberger et al., 2011). CES, BEE1, and BEE3 were shown to be expressed with distinct but overlapping patterns within the reproductive tract (Crawford and Yanofsky, 2011; Poppenberger et al., 2011). In the ces (haf) bee1 bee3 triple mutant ECM biosynthesis and PCD fail to occur within this tissue (Crawford and Yanofsky, 2011).

The stk mutant did not present any type of reproductive defect (no statistical differences when the number of unfertilized ovules was compared with wild type), while the single mutant ces-4 is characterized by a small but significative percentage of unfertilized ovules (Figure 2). However, the stk ces-4 double mutant had 28.1% of unfertilized ovules, which were almost exclusively found at the basal part of the silique (Figure 2). This phenotype suggests, similar to the ces bee1 bee3 triple mutant (Crawford and Yanofsky, 2011; Poppenberger et al., 2011), a defect in transmitting tract development.

Interestingly, we also observed a genetic interaction between STK and BEE1 and BEE3 in transmitting tract development since the bee1 bee2 bee3 stk quadruple mutant showed 9.21% unfertilized ovules at the lower part of the silique, while stk had only 6.2% and the triple mutant bee1 bee2 bee3 displayed 1.8% of unfertilized ovules (Figure 3). In the light of the fact that this phenotype was weaker than that observed for stk ces-4 (Figure 2), we suggest that CES plays a more predominant role in transmitting tract formation together with STK. To understand a possible role of SHP1 and SHP2 in the formation of the transmitting tract, the triple mutant shp1 shp2 ces-4 was obtained. This triple mutant showed 16.2% of unfertilized ovules which is also a milder phenotype when compared to stk ces-4 (Figure 2). This suggests that STK is important for transmitting tract development, although is not able to cover for 100% SHPs function since the shp1 shp2 ces-4 triple mutant presents a considerable percentage of unfertilized ovules.

Recently, it has been reported that in Arabidopsis thaliana the correct accomplishment of ovule fertilization directly influences the final length of the fruit. In the ntt mutant, in which the seeds only develop in the upper part of the fruit (Crawford et al., 2007), the authors described by using a fruit live imaging platform (FLIP) that in the upper part of the ntt siliques (where fertilization happens) the cells composing the valves elongate properly like in the wild type control while in the bottom part (where no fertilization occurs) the cells do not elongate properly (Ripoll et al., 2019). The quadruple mutant bee1 bee3 stk ces-4, which, similarly, to the ntt mutant displayed the highest percentage of unfertilized ovules from the middle to the bottom part of the fruit (Crawford et al., 2007), was characterized by a clear decrease in fruit length as can be seen in Figure 3, confirming the importance of ovule fertilization for fruit size.

The different degree of unfertilized ovule phenotypes obtained in the different mutant combinations points to redundancy between MADS-box, bHLH, and the BEEs transcription factor encoding genes in a dosage-dependent manner. The unfertilized ovules phenotype was only observed when different mutant alleles were combined and was more pronounced in the quadruple and quintuple mutants, suggesting an additive effect. Such a dosage dependent effect was already described for other MADS-box genes, like FLOWERING LOCUS C (FLC) that acts, in a dosage-dependent manner, as a potent repressor of the floral transition (Bastow et al., 2004; Michaels and Amasino, 2007; Sheldon et al., 2007).

Y2H and BIFC experiments showed that STK is able to interact with CES, but not with BEE1 (Figure 5 and Supplementary Table S3). This result suggests that CES might bridge the interaction between MADS and BEE proteins to form a hypothetical multimeric complex that controls correct transmitting tract formation. Cooperative action of MADS-domain and bHLH transcription factors has been previously described in mice. Molkentin et al. (1995) showed that MEF2 factors, a family of MADS-domain proteins expressed in muscle cells and other cell types, acted as coregulators of myogenic bHLH proteins to activate muscle genes (Molkentin et al., 1995). Moreover, the authors show that these two types of transcription factors physically interact, and that either factor can interact with the other when one is bound to DNA (Molkentin et al., 1995). Interestingly, the authors present this as a general mechanism for the regulation of transcription in specific cell types (Molkentin et al., 1995). This might also be a mechanism for the combinatorial control of MADS-box (with STK as a main player) and bHLH (CES as principal component) transcription factors for the formation of the transmitting tract tissue in Arabidopsis. Studies to understand the nature of a multimeric complex formed by MADS-box and bHLH proteins would really be a very interesting field for future studies.

When we discovered, by backcrossing, that the male part was not compromised in the stk ces-4 and shp1 shp2 ces-4 multiple mutants (Supplementary Figure S2), we decided to perform a detailed analysis of the septum morphology of the mutants displaying the most stronger phenotype related to the percentage of unfertilized ovules. The single mutants stk and ces-4 possess a similar septum respect to the wild type (Figure 4). However, when STK was mutated together with CES the fruits showed a septum fusion defect in the lower middle part in concomitance with the medial region of the ovary. The defective septum was found in the double mutant stk ces-4, in the quadruple mutant bee1 bee3 stk ces-4 and in the quintuple mutant bee1 bee2 bee3 stk ces-4. Notably, no defective carpel fusion was revealed for the quadruple mutants bee1 bee2 bee3 stk and bee1 bee2 bee3 ces-4 (Figure 4) which have a significant percentage of unfertilized ovules (Figure 3), confirming the hypothesis that only when STK and CES are mutated simultaneously the fusion of the carpels was in part compromised. Gynoecium defects were previously described only in the ntt stk double mutant, where there is a failure of the fusion with a septum that has a hole (Herrera-Ubaldo et al., 2019). The organ fusion defects were not disclosed for the single mutants stk and ntt and for the triple mutant ces bee1 bee3 (Crawford and Yanofsky, 2011; Herrera-Ubaldo et al., 2019). In the ces bee1 bee3 and probably in our quadruple mutant bee1 bee2 bee3 ces-4 the unfertilized ovule phenotype is due to defects, as previously described by Crawford and Yanofsky (2011), in ECM components synthases and PCD failure in transmitting tract. In conclusion, only when the dimer STK-CES is not formed, the carpels fusion is compromised.

These phenotypes could be explained by mis-regulation of the genes involved in cuticle formation. The cuticle is a lipidic alteration of the cell wall (Yeats and Rose, 2013). The alteration in cuticle biosynthesis, deposition and wax component transports could cause organ fusion defects in plants, but in particular in flowers and carpels (Nawrath, 2006; Sieber et al., 2007; Panikashvili et al., 2010). Intriguingly transcriptional profiles of the quadruple mutant suggest a strong downregulation of LTP2 (Table 1 and Figure 6), a gene involved in lipid transport and/or wax transport to form the cuticle (Kunst and Samuels, 2003). Notably, the LTP5 gene, expressed in the transmitting tract, which belongs to the same gene family as LTP2, and which did not show significantly altered expression in the quadruple mutant according to our analysis, has previously been shown to be involved in pollen tube guidance (Kunst and Samuels, 2003; Chae et al., 2009). We believe that also LTP2 might be involved in the correct fertilization of the ovules influencing not just the pollen tube growth but also the correct fusion of the two carpels that permits to form the septum and transmitting tract. Moreover, in the bee1 bee3 stk ces-4 quadruple mutant we detected the overexpression of three genes encoding for GDSL esterase/lipase involved in lipid catabolism (Table 1 and Figure 6). The overexpression of a GDSL named, CUTICLE DESTRUCTING FACTOR 1 (CDEF1), causes abnormal organ fusion in leaves, stems and flowers (Takahashi et al., 2010). This gene codes for a protein which is considered a cutinase (Takahashi et al., 2010). It is possible that like for the CDEF1 overexpression, the excess of cutinase in the quadruple mutant respect to wild type caused abnormal organ fusion.

Cell death and tissue degeneration are required for pollen tube growth in the transmitting tract of Arabidopsis (Crawford et al., 2007). The phenotype described for combined ces and bee mutants showed that they fail to form the ECM and to accomplish PCD within the transmitting tract (Crawford and Yanofsky, 2011). Interestingly, STK has already been shown to control PCD in transmitting tract with NTT (Herrera-Ubaldo et al., 2019) but also in the receptive synergid cell through regulation of its direct target genes, VDD and VAL, two Reproductive Meristem (REM) transcription factors (Matias-Hernandez et al., 2010; Mendes et al., 2016).

As we show in the Table 1, in the quadruple mutant bee1 bee3 stk ces-4 several genes involved in cell death are downregulated. The ATWSCP gene encodes for a Kunitz-type protease inhibitor localized preferentially in the transmitting tract and septum, which is involved in ECM and PCD but also in the exit of pollen tubes from the septum to reach the ovules during the fertilization phase (Bektas et al., 2012; Boex-Fontvieille et al., 2015). Also, ACD6 is downregulated (Table 1 and Figure 6), this gene encodes for a protein that promotes cell death in Arabidopsis (Rate et al., 2007). Furthermore, three Arabinogalactan genes (AGPs) resulted downregulated in the quadruple mutant (Table 1 and Figure 6). The AGPs have been associated with ECM composition of the transmitting tract in different plant species (Hoggart and Clarke, 1984; Webb and Williams, 1988; Gane et al., 1995) and they have been shown to play a role in the regulation of PCD during pollen tube growth (Coimbra et al., 2009).

It was recently suggested that the auxin pathway plays an important role in the specification and formation of the tissues along the pistil axis (Larsson et al., 2013; Robert et al., 2015). However, the role of auxin signaling in transmitting tract development is still not fully understood. ARF6 and ARF8 regulate CES expression in the medial domain of the pistil during transmitting tract development (Crawford and Yanofsky, 2011). The double mutant arf6 arf8 displayed reduced alcian blue staining, that detects acidic polysaccharides, the major components of the ECM (Crawford and Yanofsky, 2011). According to our RNAseq data, four Aux/IAAs and two SAURs genes were downregulated in the quadruple mutant in respect to wild type (Table 1 and Figure 6). Aux/IAA proteins and SAURs might be activated in presence of brassinosteroids since evidences of a crosstalk between auxin and brassinosteroid have already been described (Nakamura et al., 2006). Moreover, SAUR9 expression is induced by a combination of auxin and brassinosteroids (van Mourik et al., 2017). Probably a STK-CES-BEE1-BEE3 complex is able to connect the auxin-brassinosteroid crosstalk in transmitting tract tissue, however, the precise role is of this crosstalk in the regulation of transmitting tract development remains unclear.

Finally, from the analysis of the enriched binding sites an interestingly candidate was found to be a possible direct target of the complex composed by the MADS and bHLH transcription factors described in this work. The AT3G49270 gene (upregulated, see Supplementary Table S5) has both CArG and G boxes binding sites. It encodes for an extensin-like protein that has been suggested to have a possible role in postpollination responses such as cell wall protein modification or degradation associated with pistil receptivity (Boavida et al., 2011).

In summary, the data presented here suggest that a MADS-box-bHLH protein complex mastered by the physical interaction of STK-CES, plays a role in the regulation of genes involved in the PCD, ECM and cell wall biosynthesis pathways and by that controlling correct carpel fusion and transmitting tract development to guarantee efficient ovule fertilization (Figure 7).