Different expression pattern of flowering pathway genes contribute to male or female organ development during floral transition in the monoecious weed Ambrosia artemisiifolia L. (Asteraceae)

Plant materials and cultivation conditions

The RNA-seq experiments were based on common ragweed samples from four bulked individuals. Seeds were collected from Cserszegtomaj, Zala county, Hungary (GPS: 46.79528, 17.26005) and were surface sterilized for 1 min in 70% (v/v) ethanol and washed with deep water for 23–30 min, followed by dipping in a 7% solution of calcium hypochlorite with TWEEN 20 for 20 min, and rinsing three times with sterile distilled water. The sterilized seeds were then germinated on MS medium (Murashige & Skoog, 1962) with 0.3 mg/l metatopolin hormon supplemented with 2% sucrose/litre and solidified with 0.8% phytoagar. The pH was adjusted to 5.8 prior to autoclaving at 120 °C for 20 min. All the cultures were maintained at 25 °C under a 16/8 light/dark photoperiodic regime with a light intensity of about 2,400–2,800 lux. Node segments (1.5–2 cm long) were excised from 21-day-old plants and transferred to above described MS medium.

Nine developmental stages of female flowers were collected from 12 in vitro cultivated individuals of common ragweed. The classification of the different phenological phases was based on visual observation (Fig. 1). For RNA-seq two categories of nine phenological phases (1F and 2F) were established and combined such as early 1F (Figs. 1A–1C) and late 2F (Figs. 1G–1I) developmental stages. Male flower buds were collected from genotypes used in F, M, L transcriptoms (Virág et al., 2016). For gene expression comparision we used 1M and 2M libraries as indicated in the Figs. 1J–1M, 1P, 1Q and 2M). In order to equal distribution, we determined the weight of the samples from each developmental stages (Table S1).

Along RT-qPCR analysis we used four additionally wild growing genotypes (collected in Keszthely, Zala county, Hungary, GPS: 46.7654716, 17.2479554).

RNA-seq library construction and sequencing

Samples were frozen in liquid nitrogen immediately after collection. RNA extraction and on-column DNase digestion were performed from each categories using TaKaRa Plant RNA Extraction Kit (Takara Bio Inc, Kusatsu, Japan) following the manufacturer’s instruction. RNA integrity and quantity were assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Total RNA of 2.0 μg from 1F, 1M, and 2F, 2M samples were used as starting material for cDNA synthesis and library construction. Enrichment of mRNA, cDNA synthesis and library preparation for Illumina NextSeq paired-end sequencing were carried out using TruSeq™ RNA sample preparation kit (Low-Throughput protocol) using an oligo(dT)₁₈ as described by Virág et al. (2016). The RNA-Seq was performed using Illumina NextSeq5000 system. Raw sequence data in FASTQ format were deposited in the NCBI sequence read archive (SRA) database under the accession SRR5965731 (1F) and SRR5965732 (2F). We used 1M and 2M for gene expression comparison in the case of investigated genes. 1M and 2M transcriptom datasets and their analysis will be published elsewhere.

Sequence reads of male inflorescence (M), leaves (L), and female (F) inflorescence containing all nine developmental stages were used from earlier reported and deposited data in NCBI SRA under the accession SRP08007 (Virág et al., 2016).

Construction of Ambrosia artemisiifolia transcript datasets, identification of differentially expressed sequences and filtering pistillate specific transcripts

During preprocessing raw reads quality were examined using the FastQC quality control software (Andrews, 2010). Based on FastQC report, sequences found to be represented more than 0.1% of the total and low-quality bases (corresponding to a 0.1% sequencing error rate) were removed and trimmed using a self-developed application “GenoUtils” written in Visual Studio integrated development environment in C#. Conversion of .fastq to .fasta files was also performed using this in house application.

A sample-specific multiple assembly of cleaned reads from five libraries (M, L, F, 1F, and 2F) were performed. De novo assembly of each sample type were performed by using Trinity (Haas et al., 2013). For reference guided alignments Bowtie2 short read aligner (Langmead et al., 2009) was applied. The assembly strategy was performed in a way with taking into account the refinement of pistillate tissue specific sequences in both developmental classes as the final outcomes. Assembly strategy is summarized in the flow chart (Fig. 2). In detaile. First, a combined read set was assembled from the three sample libraries (M, F, and L) to generate a “reference” de novo transcriptome assembly that was deposited in the NCBI TSA database under the accession GEZL00000000 and reported in our prior study (Virág et al., 2016). During this processes, reads with certain length of overlap were combined to form contigs (kmer size, K = 25). This combined dataset was used as reference shotgun assembly in the further alignments. In order to perform sample-specific expression analysis, we screened the DE sequences by aligning the original sample reads to the reference followed by abundance estimation using RSEM (Li & Dewey, 2011). The resulting DE transcripts were further clastered according to their expression patterns by applying Microsoft SQL Server Management Studio. Protein coding regions were extracted from the reference assembly using TransDecoder and further characterized according to likely functions based on sequence homology or domain content using BLAST+ (Altschul et al., 1990). In this way coding sequences (CDS) as unique transcripts of each five samples were screened and compared. Separately, de novo assemblies of 1F and 2F libraries were performed using Trinity and deposited in the NCBI TSA database under the accession GFWB00000000 (1F) and GFWS00000000 (2F). Using these reference contig sets 1F and 2F unique transcripts were realigned applying BLASTn with E value less than 10⁻⁵. The resulting narrowed and specified sequences were annotated using NCBI nr protein database.

Selection of genes for floral transcriptome characterization

Genes representing different floral regulatory pathways—stimulated by environmental and endogenous signals—were selected for characterize genetic events during male and female floral development. Architecture formation were analized by determining ABC(E) genes expression pattern at generative growth (differences in FOIGs). Genes that enable and promote floral meristem initiation during vegetative growth (FMIGs) such as photoperiodic, vernalization, GA pathway, and flower development genes (Figs. 3 and 4) were also investigated in both wild growing and in vitro samples. Normalized expression of selected genes were performed digitally and expression values were compared.

Digital gene expression analysis

The obtained reads from each library were mapped to the selected gene CDSs by using Bowtie2. The mapped reads were used to estimate the transcriptome level by the reads per kilobase per million mapped reads (RPKM) method (Mortazavi et al., 2008): $$\mathrm{R}\mathrm{P}\mathrm{K}\mathrm{M}={\displaystyle \frac{\mathrm{C}\mathrm{D}\mathrm{S}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}\ast {10}^9}{\mathrm{C}\mathrm{D}\mathrm{S}\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{h}\ast \mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}\mathrm{m}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{d}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}}}$$In this equation CDS represent the CDS of the investigated gene. Read count means the mapped reads from the total. CDS length is the nucleotid sequence length of coding region of the investigated gene and total mapped reads are the total raw reads of given cleaned library. For automatic calculation of RPKM values of genes a self developed pipeline as modul of GenoUtils were applied.

Gene expression analysis with RT-qPCR amplification

Reverse transcription was performed starting from one µg total RNA using Maxima H Minus First Strand cDNA Synthesis Kit with dsDNase (Thermo Scientific, Waltham, MA, USA) according to the manufacturer’s protocol, using an oligo(dT)₁₈ and random hexamer primers (Thermo Scientific, Waltham, MA, USA), final volume was 20 µl. cDNA (one µl) was used for real-time PCR amplification on a Bio-Rad CFX96 System. qPCR analysis (and efficiency) was performed with one μl of cDNA on a Bio-Rad CFX96 System using Xceed qPCR SG Mix (Institute of Applied Biotechnologies, Praha, Czech Republic). The relative gene expression was calculated with ΔΔCt method using Bio-Rad CFX Manager™ Software v3.1. Primer efficiency was analyzed with CFX Manager™ Software v3.1 (Bio-Rad, Hercules, CA, USA). PCR was performed as follows: an initial activation of the polymerase enzyme at 95 °C for 2 min was followed by 46 cycles at 95 °C 5 s, 59 °C for 30 s, ended with a 5 s melt analisys ranging between 65 and 95 °C with 0.5 increments. Primers (Table S2) were designed with Primer3 (Koressaar & Remm, 2007; Untergrasser et al., 2012) based on in silico sequence prediction of CDS.

Gene expression analysis was established based on three technical and biological replicates, and normalized with the reference gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (Hodgins et al., 2013). In a pilot experiment GAPDH was selected as reference gene. In this investigation based on the literature CDS of 15 housekeeping genes were determined. Subsequently the relative expression levels were determined in three sample types (L, M, F) were six reference genes showed the same level, respectively. Among the equally expressed genes three reference genes (GAPDH, TUA, TUB) were selected for RPKM value calculation (Table S3) and PCR experiments. GAPDH was the most characteristic along these investigations. Based on RPKM values GAPDH expression referd on equivalent level order of magnitude level in each liberies ranging from 2,480 to 4,914.

Annotation of DE transcripts

Whereas Ambrosia artemisiifolia floral transcriptome was a newly targeted transcriptome the annotation procedure was performed as described by Haas et al. (2013). Potential coding regions within reconstructed transcripts that were insufficiently represented by detectable homologies to known proteins were predicted based on metrics tied to sequence composition by applying TransDecoder include with Trinity. Running this application on the Trinity-reconstructed transcripts the candidate protein-coding regions could be identified based on nucleotide composition, open reading frame (ORF) length, and (optional) Pfam domain content. To predict their functions the latest non-redundant protein database 10/03/2016 (ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz) were used.

For further annotation of unigenes using various bioinformatics approaches, the unigenes were firstly searched against the Nr database and the Swiss-Prot protein database using local BLASTx (Gish & States, 1993) with E value cutoff of 10⁻⁵. With Nr annotation, Blast2GO (Götz et al., 2008) was used to get gene ontology (GO) annotation according to molecular function, biological process and cellular component ontologies (http://www.geneontology.org).

RNA sequencing and assembly of Ambrosia artemisiifolia floral libraries

The M, F, and L libraries were obtained through HiSeq2000 platform. Library construction and assembly was described detailed in our prior study (Virág et al., 2016). In brief, after cleaning we obtained 18,472,374 (M); 15,290,200 (F) and 17,435,976 (L) from total raw reads of 24,110,256; 23,264,636 and 24,330,693 (2 × 100 bp). Assembly of combined read sets resulted in 229,116 transcripts and 162,494 trinity genes (unigenes) and used as a reference (Table S4).

A total of 39,664,366 (1F) and 37,127,852 (2F) 2 × 80 bp reads were generated by the Illumina NextSeq500 system. High quality (> Q20) bases were more than 91% in (paired-end) reads of both sample. Percentage of unresolved bases (Ns) and overrepresented sequences was observed to be minimal 0.197% and 0.0004% in average. In total of 92% were high quality sequences from raw data: numerically 36,491,216 (1F) and 34,157,623 (2F). These high quality, processed paired-end reads were used to assemble into contigs and transcripts. Their de novo assembly resulted in 109,452 (1F) and 97,239 (2F) contigs. In comparision the reference guided assembly led to 40% more contigs 147,457 (1F) and 141,58 (2F), respectively. However, N50 length showed a minor 5–11% dispersion between the ranges 626.9–697.5 in both methods. Sequencing statistics are summarized in Table S4.

Selection of Ambrosia artemisiifolia unique and differentially expressed genes

In order to identify exclusively expressed transcripts in the sex-specific floral tissues we used Bowtie2 alignments of F, M, L, 1F and 2F libraries for the Ambrosia reference transcriptome GEZL00000000. In order to prove the 1F and 2F exclusive transcripts in female sample (F), contigs were realigned using BLASTn algorithm to the 1F and 2F de novo Trinity transcripts (GFWB00000000, GFWS00000000). Filtering method is summarized in Fig. 2. The transcript datasets were extracted from Trinotate heat map and postprocessed using Microsoft SQL Server Management Studio (Dataset S1). The queries for unique organ specific sequences resulted in 5,659 (M), 1,691 (F), and 4,267 (L). The number of exclusive and shared transcripts in wild growing male and female flowers and leaf are visualized in Venn diagram. For visualization we used FunRich (V3) software (Pathan et al., 2015) (Fig. 5). The number of annotable unique transcripts were 11,300 in M, F, and L libraries. Among which the most important organ specific genes are summarized in Tables 3 and 4 and their sequences are in Dataset S2. The presented genes were within the upper limit (>80%) of the total expression level value.

Annotation and functional classification of Ambrosia artemisiifolia unigenes

Functional annotation of the transcriptome sequences and analysis of the annotation of three-specific RNA pools (M, F, L) were performed using BLASTx search (Gish & States, 1993) and Blast2GO (Conesa et al., 2005). This process resulted in 81.56% translated contigs in total. Based on GO annotation the cell part, organic cyclic compound binding and organic substance metabolic process were the most abundant GO slims in the categories of cellular component, molecular function, and biological process, respectively (Fig. S1). The high physiological activity of Ambrosia reproductive tissues was suspected from the high representation of primary metabolic process and cellular metabolic process groups within the biological process category.

Representation of the significant sub-ontologies of cellular component, molecular function, and biological process are summarized in Fig. S1. GO enrichment analysis were performed with the investigated floral genes wich resaults are summarized in Figs. S2–S4. Statistics and background information can be found in Datasets S4 and S5.

Validation of the Ambrosia artemisiifolia floral transcriptome

The identification of differences between the male and female floral transcriptomes was based on in silico analysis. We performed a local blast alignment by using the determined 80 CDS fasta sequences (MK098047–MK098126) against male (M) and female (F) transcriptoms. No qualitative differences were found among the sequences of commonly expressed genes in two floral organs. The blast output files are in Dataset S3. The observed expression differences were validated investigating the transcript expression levels of homolog genes responsible for pollen and embryo formation in model organism Arabidopsis thaliana. The reliability of male transcriptome was characterized by overexpression of the homologs of AMS, LAP3, LAP5, LAP6, MSE1, PME, PG1, PG2. The female transcriptome validity was justified by overexpression of YAB1, YAB4, and ANT transcription factors. Additionally, we introduced SUPERMAN (SUP, MK098096, GO:0003676), Transcription repressor MYB5 (MYB5, MK098099, GO:0003677), Transcription repressor MYB61 (MYB61, MK098101, GO:0003677) and Protein BRANCHED2 (TCP12, MK098098, GO:0003700) as female candidate genes based on our filtering by Microsoft SQL Server Management Studio. These regulators are involved in plant and floral tissue development, but their role during sex determination is still unclear. Expression differences and function of organ specific genes are summarized in Tables 1 and 2.

RT-qPCR validation of tissue specific genes

For the RT-qPCR validation we selected the genes showing the largest differences (Tables 1 and 2) and most interesting from functional aspects. TCP12, SUP, and LAP6 were selected to justify their exclusive expression in female and male flowers. STIG1 were chosen according to its opposite expression in male flowers. Additionally we selected also PI (MK098083, GO:0000977) as ABC transcript exclusively present in staminate tissues.

Gene expressions during floral meristem initiation (vegetative growth)

Functional and structural analysis of Arabidopsis thaliana gene homologs were investigated to elucidate floral regulation in the common ragweed. Expression pattern of photoperiodic, vernalization, GA, and flower development pathways comparing with floral architecture formation are summarized in Figs. 3 and 4.

Gene expressions during floral architecture formation (generative growth)

The ABC(E) genes responsible for the formation of flower organs were determined on the basis of various compositae species (Gerbera sp, Helianthus sp, Tagetes sp, Chrysanthemum sp) found in public databases (NCBI, UNIProt). The CDSs and their expression patterns of the following 13 gene homologs were determined: AP1 (MK098082), AP2 (MK098081), AP3 1 (MK098086), AP3 2 (MK098085), AP3 3 (MK098084), PI, AG (MK098087), SEP1 (MK098094), SEP2 1 (MK098092), SEP2 2 (MK098091), SEP3 1 (MK098090), SEP3 2 (MK098089) and SEP4 (MK098088) (Figs. 3 and 4). All of investigated floral organ genes except AP2 (GO:0003677) have the same GO ID: GO:0000977.

Transcription factors differentially expressed in F, M, and L libraries

Based on analysis of the different TF groups that are sex-specific and DE in each library the following factors were identified in our samples according to Rocheta et al. (2014). Gene expression values of all investigated libraries and functional descriptions are summarized in Tables S5 and S6.

Expression of female flower specific genes

In order to identify exclusively expressed genes in female flowers, the combined assembly of five (M, L, F, 1F, and 2F) non-normalized libraries were analyzed. In the analysis, 1,691 unique transcripts were found in the wild growing F library, representing the nine developmental stages. Transcripts longer than 700 bp were further investigated. After this filtering, we found 60 transcripts showing evaluable ORF, among which nine genes were well annotable using NCBI nr and Swiss-Prot databases. Based on this, the following protein coding and TFs were identified: Protein kinase 1A (PBL9, MK098095, GO:0004674), SUP, Alpha carbonic anhydrase-7 (ACA7, MK098097, GO:0004089), Laccase-2 (LAC2, MK098102, GO:0005507), TCP12, MYB5, MIB61, Zinc finger protein WIP2 (MK098100, GO:0003677), and Protein Brassinosteroid enhanced expression 1 (BEE1, MK098103, GO:0006355) (see Table 3). In order to find genes that might correspond to early and late flower development, exact CDSs of these genes were determined in silico and validated with sanger sequencing from in vitro plant materials. RPKM values were calculated in F, L, M, 1F, 2F, 1M, and 2M libraries using validated CDSs.

Expression of male flower specific genes

The uniquely expressed transcript number in male tissues were 4,549 (Fig. 5) of which sequences longer than 700 bp were filtered out. In this way total CDSs of 41 genes were identified based on NCBI nt Blast and NCBI ORF finder databases. The 12 most characteristic male specific genes, RPKM values and functional properties are summarized in Table 4.

Validation of the Ambrosia artemisiifolia floral transcriptome

As expected, the investigated genes showed differential expression ratios with a meaningful results according that male tissues were characterized with an overexpression of genes (AMS, LAP3, LAP5, LAP6, STIG1, PG1, PG2, PME, MSE1) taking part in pollen structural development, fatty acid and exine constituent biosynthesis during the anther and pollen differentiation. In contrast, expression of these genes was completely absent in the pistillate tissues in which ovule polarity, initiation, structural development and embryogenesis influencing genetic elements (YAB class TFs) showed about 67% higher presence than in staminate samples. The non-exclusive occurrence of these sequences in pistillate flowers may be explained by the contribution to setting up the polarity in stamen maybe through their action in supporting the development of the floral meristem or/and regulating margin growth (Eshed et al., 2004) that was observed also in all lateral organs (Chen et al., 1999; Siegfried et al., 1999). Additionally, the validation method was confirmed through the very low incidence of candidate sequences in leaf transcriptomes.

RT-qPCR validation of tissue specific genes

To evaluate the results, genes were grouped according to their specificity. According to this, the female transcriptome was validated based on TCP12 and SUP expressional pattern. Female specificity of SUP and TCP12 was confirmed (Fig. 6B). The expression of TCP12 was three times higher than of the SUP in female flowers confirming the in silico RPKM results. Male transcriptome was validated according the expression pattern of LAP6, STIG1, and PI. The non exclusive expression of PI was observed indicating the ABC function and B class activity in both flower types. Since in female flower the petal differentiation is vestigial and this whorl is united into a single tubular flower, therefore PI homolog transcript are less present than in male flowers where petals are united into short anthers with long stigma and stamen structures represents a large amount.In silico expession of STIG1 homologe in male flowers (Table 1) were confirmed by RT-qPCR, indicating the modified pistil origin in the pistillodium structure. Howerver this transcript showed a minimal expression also in female flowers (Fig. 6A). This fenomena reveals the common origin of these flower types which were separated during evolutuion.

Gene expressions during floral meristem initiation (vegetative growth)

Plant development and architecture is regulated by meristems that initiate lateral organs on their flanks. Gene regulatory networks that control the transition of a vegetative shoot apical meristem into an inflorescence meristem (FMIGs), together with those necessary to specify floral meristem identity (FOIGs) have been explored in Arabidopsis thaliana and are highly complex and redundant. We have evaluated and report the most prominent cases in this study. Since the morfogenesis of the in vitro cultivated female inflorescence of common ragweed is more intensive than in wild-growing plants, therefore genetic regulation of the pistillate inflorescence can be more effectively investigated in these samples. The investigated pathways and related genes showed the following correllations. No interaction was found between LFY, FRIGIDA (FRI) and FLC transcription factors that determine the initiating steps during vegetative growth. Homologous transcript of this genes were not found in any of the investigated libraries, and were not expressed in the examined samples. This indicates, that the blossom phase has exceeded the initial signals in the bio-rythm of collected samples. The combined expression of LFY and FLC leads to the so-called “late-flowering phenotype” (Michaels & Amasino, 2001) absent in the 1F and 2F individuals confirming the complete lack of expression of these genes, that may be traced back to the altered hormonal conditions (cytokinin-supplemented culture medium) during in vitro cultivation.

Overexpression of photoperiodic pathway related transcripts of PHYB (MK098056, GO:0000155) in 1F, 2F, 1M and 2M libraries and COP1 (MK098059, GO:0009785), CO (MK098060, GO:0005634) and FT (MK098073, GO:0008429) in the L library were found. Significantly higher expression of PHYB than PHYA (MK098057, GO:000155) were observed in all samples independently of the cultivation conditions indicating the red (morning) illumination affects more intensively the flowering processes in Ambrosia artemisiifolia. These results confirm the findings of Deen, Swanton & Hunt (2001) at genetic level, according which Ambrosia artemisiifolia belongs to the short day plants. The transcription of FT is activated by CO in leaf, and the FT protein moves to the shoot apex to induce flowering, therefore we found these transcripts overexpressed only in leaf tissues in wild-growing samples. Since Ambrosia artemisiifolia is a short day plant, higher expression of COP1 in leaf may be explained by the function to repress CO degradation (Song, Ito & Imaizumi, 2010). Interestingly, FT-mRNA was only detected in early stage inflorescences (1F) suggesting that the initial flowering signal loses its regulatory function during female flower development, and TFs responsible for floral organ development come to the fore as control elements.

The vernalization pathway gene homolog VRN1 (MK098069, GO:0003677), showed high values in all developmental stages (1F, 2F, 1M, 2M), which function is to inhibit the flowering repressor FLC, thus maintain generative morphogenesis of the flowering meristem. Expression of VRN1 was increased during female flower maturation (Sung & Amasino, 2004). We found SOC1 (MK098072, GO:0000957) homolog in 1F and 2F libraries. SOC1 can integrate flowering signals from different pathways and interact with FT-mRNA to develop inflorescence. The SOC1-FT interaction, based on gene expression values, showed intensive regulation in the early stages of the common ragweed (1F).

SOC1 associating with AGL24 (MK098093, GO:0000977) mediates effect of GAs on flowering under short-day conditions, and regulates the expression of LFY, which links floral induction and floral development (Liu et al., 2008). The lack of SOC1, LFY (RPKM = 0) and overexpression of FT (RPKM = 1,487) in wild-growing samples indicated that flowering transition was regulated more intensively by the photoperiodic pathway in this samples. In contrast, hormonal pathway related genes (SOC1, GID1A (MK098047, GO:0016787), GAI (MK098050, GO:0003712), Ga2ox8 (MK098051, GO:0016491), CAL (MK098071, GO:0000977), AGL24) were more expressed in in vitro plants. Therefore we concluded, that both cultivation condition reached the AP1 expression that is the “first pass” of floral architecture formation, however, the vegetative pathway was altered under in vitro conditions. This phenomena has led to altered morphology and genesis of floral genders. The female inflorescences grow faster and in greater abundance in the case of in vitro individuals than in the wild-growing plants (see Fig. 7). Male inflorescences, however, were barely appeared in vitro. Under wild-growing conditions male flowers dominate because of the intense lightening mediated signals.

The effective integration of hormonal signaling in plant vegetative growth is essential in which GAs and cytokinins act antagonistically in leaf formation and meristem maintenance (Gan et al., 2007). Cytokinins may act antagonistically with GA preventing LFY induction (Li, Li & Smith, 2017) that was observed in 1F and 2F samples in our flowering system. However, their integrative role in flower formation is still contested. We demonstrated the importance of cytokinins as part of hormonal regulation during Ambrosia floral transition (overexpression of TSF in L library), however, no sex-specificity was observed (low RPKM values in 1F, 2F, 1M, 2M libraries). We demonstrated also that cytokinin signaling may act positively on GA-pathway genes during flowering therefore GA pathway genes (GID1A and GID2 (MK098049, GO:0005515)) were found overexpressed in female inflorescences.

Gene expressions during floral architecture formation (generative growth)

The AP2, AP3 2, and PI genes were not or very low expressed in the female samples, therefore we considered them as male specific genes of common ragweed. PI with AG determine the stamen and AP2 with AP3 determine the petal formation that are male flower elements in the common ragweed (Fig. 8A). The expression of the AP1 gene showed a significantly higher value in female flowers than in the male inflorescence, suggesting that the initiation of petals and stamens is stronger in the pistillate flowers. This overexpression may be traced back to the hormonal pathway domination and through SOC1 during vegetative phase in female flowers (Kaufmann et al., 2010). In the early stage of development of female inflorescences (1F) a high expression rate were found for SEP1, SEP2, SEP3, and SEP4. Supposedly these genes are responsible for the development of petal and carpel primordias. Hovewer, SEP3 2 and SEP2 2 showed high expression in the late phase (2F) that may be involved in the pistill prolongation (Fig. 8B).

Transcription factors differentially expressed in F, M and L libraries

Transcription factors play important roles in plant development and flower morphogenesis responding to the environment. These factors interact specifically with sequences located in the promoter regions of the genes they regulate. TFs are classified in families according to the structure of their DNA-binding domain.

Based on RPKM values, homologs of CONSTANS-LIKE zink-finger TF family was found to be characteristic in both female and male samples (COL4, MK098124, GO:0005634, COL5, MK098125, GO:0005634). Upregulation of COL4 was found in 1F and 2F libraries indicating a stronger interaction of female tissues with the circadian clock and light signals than in male flowers. COL4 null mutants flowered early under short or long days. In contrast, OsCOL4 activation-tagging mutants (OsCOL4-D) flowered late in either environment (Lee et al., 2010). Of these the late flower showed higher RPKM value in the correspondence of CO transcript level. Down regulation of COL9 (MK098126, GO:0005634) was observed in all wild-growing samples. Since overexpression of COL9 also delayed flowering by reducing expression of CO and FT (Cheng & Wang, 2005), our results confirm the natural circadian rhythm of wild growing samples. Over-expression of COL5 can induce flowering in short-day grown Arabidopsis (Hassidim et al., 2009). We observed the higher COL5 expression in female flowers, mainly in the case of in vitro plants where pistillate flowering showed intensive growth under short light conditions.

The next most characteristic transcript was an ILR3 (MK098117, GO:0046983) homolog of IAA-LEUCINE RESISTANT 3 protein representing the helix-loop-helix protein family. ILR3 is required to maintain Fe homeostasis in correlation with a large amount of Fe³⁺ observed in the sepals in A. thaliana (Zhang et al., 2015). Since this transcript was found to be expressed roughly equally in all investigated samples (1F, 2F, 1M, 2M, M, and L), approximately, we assumed that to keep at appropriate level the iron homeostasis is essential to all flower types in fact to maintain hormonal and redox balance for fertility (Sudre et al., 2013).

Among the hormone related factors the auxin responsive AUX/IAA protein family genes such as IAA27 (MK098121, GO:0005634) and IAA9 (MK098120, GO:0005634) were found to be characteristic in both reproductive tissues. These transcriptional factors are repressors of early auxin response genes at low auxin concentrations (Liscum & Reed, 2002). RT-qPCR and in situ hybridization have shown that tissue-specific gradient of IAA9 expression which was established during flower development, and the release of which, upon pollination and fertilization, triggers the initiation of fruit development (Wang et al., 2009). These results explain the higher expression rate in female tissues in both early and late phenophases in Ambrosia. The auxin efflux carrier homolog, PIN-FORMED 1 (PIN1, MK098122, GO:0009734) was also found to be upregulated in early female flowers (1F) indicating the female tissue determination might be under strong control of auxin. These results correspond to the study of Ceccato et al. (2013) where PIN1 expression and cellular localization in the ovule were essential for correct auxin efflux in early stages of female gametophyte development. In the ethylene response pathway the protein ETHYLENE INSENSITIVE3 (EIN3, MK098123, GO:0003677) regulates ethylene responsive genes as a positive regulator and results in flowering delay via repression of the floral meristem-identity genes LFY and SOC1 (Achard et al., 2007). In our flowering system it was upregulated in M, 1M, 2M libraries resulting no expression of SOC1 and LFY in male tissues (see Figs. 3 and 4). No LFY homolog was identified in any of the floral libraries in corellation of higher expression of this TF. Wild growing Ambrosia individuals grow male flowers intensively in late summer, which is a delayed state compared to the female blossom. However, SOC1 transcript was found only in the in vitro samples referring to the altered (earlier) floral regulation we assumed that this fenomena is based on changes in hormonal pathways.

NAC-domain containing TFs are essential for normal plant morphogenesis which representatives are the CUP-SHAPED COTYLEDON proteins encoded by CUC1 (MK098118, GO:0003677) and CUC2 (MK098119, GO:0003677) genes less expressed in our flower samples. Mutations in Arabidopsis CUCl and CUC2 caused fusion of cotyledons, sepals, and stamens suggesting that these factors are responsible for mechanisms that separates organ developing at adjacent positions (Aida et al., 1997). No expression of CUC1 was detected in wild growing pistillate tissues and its expression was lowest among the investigated TFs in staminate samples. In in vitro samples CUC2 was expressed only in early phenotypes (1F) in female flowers with a minimal expression rate. Since these factors are involved mainly in early development of cotyledonary primordia (Aida et al., 1997), the phenomena of downregulation of these genes may be explained by the flowering separation that has overtaken this critical state in our samples.

The homeodomain like MYB transcription factors are involved in the control of the cell cycle of plants. The majority of them are miRNA regulated and are involved in numerous cell biochemical pathways like hormonal regulation. The GA pathway related MYB33 (MK098116, GO:0003677) homolog was expressed only in male flower in wild growing Ambrosia samples in a correspondence with Rocheta et al. (2014) and Millar & Gubler (2005) in which studies the GAMYB-Like MYB33 was proved to facilitate anther development redundantly. In in vitro originated 1F and 2F libraries homolog transcripts of MYB33 were also found indicating repeatedly altered hormonal regulation of this system.

Expression of female flower specific genes

As expected, all of these investigated genes are not expressed in M, L, 1M, and 2M libraries (RPKM = 0) except of BEE1 in leaf samples with a very low value. In early female flower development the serine/threonine protein kinase PBL9 (MK098095, GO:0004674) showed an exceptionally high value. SUP, TCP12, and MYB5 were expressed two times greater than in matured flowers before the ovule formation. Two catalytic protein genes LAC2 and ACA7 were also found as unique genes in late 2F samples with low expression suggesting a not characteristic role, but taking part in necessary metabolic processes during floral maturation

PBL9 (alternatively APK1) is very weakly expressed in flower tissues. It was described as Arabidopsis thaliana serine/threonine protein kinase that phosphorilates tyrosine, serine and threonine (APK1) during signal transduction (Hirayama & Oka, 1992). The APK family is negatively regulated by the floral homeotic protein AG (AGAMOUS) (Ito et al., 1997) involved in the control of organ identity, which expression was equal in early and late female libraries, respectively. Unique expression of PBL9 in female flowers indicated that the development of staminate and pistillate tissues are regulated by different kinases at protein level during signal transduction mechanisms. Similarly, a receptor-like protein kinase was also reported as female specific gene in the monoecious Q. suber (Rocheta et al., 2014). Additionally, different expression of these gene showed altered protein phosphorilation during female flower formation, accurately it is more significant at the beginning of the pistillate flower development.

Despite of SUP is not a sex determination gene (not located on the sex—Y—chromosome) it was considered as female specific gene. Their function in female flower differentiation pathway was proved in dioceious plant Silene latifolia. During the development of the female flower in Silene latifolia, the expression of SUP is firstly detectable in whorls 2 and 3 when the normal expression pattern of the B-class flowering genes was already established and persisted in the stamen primordia until the ovule had matured (Bowman et al., 1992; Kazama et al., 2009). It is probably controlled by sex determination genes on the Y chromosome in dioecious plants.

The role of MYB transcription factor superfamily in plant reproductive development is unequivocal. MYB proteins include a conserved domain, the MYB DNA-binding domain regulating a variety of plant-specific processes and they may play roles in different plant species (Ambawat et al., 2013). For example R2R3 MYB, and MYB88 regulates female flower reproduction in A. thaliana (Makkena et al., 2012), MYB35 (MSE) regulates male specific expression in Asparagus officinalis (Murase et al., 2017). MYB1 and MYB16 regulates petal development in Petunia hybrida and A. thaliana (Baumann et al., 2007; Noda et al., 1994). Overexpression of MYB61 was demonstrated in several aspects of plant growth and development, such as xylem formation and xylem cell differentiation, and lateral root formation (Arsovski et al., 2009; Liang et al., 2005; Matías-Hernández et al., 2017; Romano et al., 2012). Importance in seed coat formation, trichome morphogenesis and mucilage secretion of MYB5 and MYB61 was also described by several studies (Arsovski et al., 2009; Gonzalez et al., 2009; Li et al., 2009; Liu, Jun & Dixon, 2014) indicating overexpression in seed coat and siliques. In our investigations MYB5 and MYB61 showed unique expression in female flowers, however, MYB5 (MK098099) MYB61 (MK098101) was found to be more characteristic in early and late phenophases, respectively.

TCP12 belongs to the plant-specific TF family TCPs. The founding members of this family have been characterized to be involved in growth, cell proliferation, and organ identity in plants (Yang et al., 2012). TCP1-like genes are involved in floral development with different regulatory and evolutionary mechanisms underlying the diverse forms of floral symmetry, such as dorsal/ventral dosage effect (Busch et al., 2011). We found that the female specific TCP12 was stronger expressed in earlier phenophases suggesting a TCP12-related zygomorphic regulation at the beginning of floral morphogenesis in this species.

The transcription factor BEE1 is a positive regulator of brassinosteroid signaling (BR) regulating plant development and physiology. The BR signal is transduced by a receptor kinase-mediated signal transduction pathway which ultimately results in altered expression of numerous genes. During flowering brassinosteroids stimulate floral transition through positively regulating the circadian clock pathway genes and inhibiting FLC repressor (Domagalska et al., 2007; Zhu, Sae-Seaw & Wang, 2013). Deficiency of this phytohormones lead to the delay of flowering and because of the disruption of the pollen tube, BR deficient plants are practically male-sterile (Clouse, 1996). Characteristic overexpression of the BR enchancer BEE1 uniquely in late-female flowers suggest on one hand an altered regulation of BR signal transduction than in male flower development, on the other hand it confirms an enhanced BR mediated floral transition related to circadian clock in Ambrosia artemisiifolia.

Transcription factor WIP2 is a zink finger protein and was found to be responsible for transmitting tract formation in the gynoecium, which is important for pollen tube growth (Crawford & Yanofsky, 2008). WIP2 plays role in replum development and regulation of AGL8/FUL (MK098076, GO:0000957), which is required for normal pattern of cell division, expansion and differentiation during morphogenesis of the silique (Marsch-Martínez et al., 2014). Expression of this gene was equal in early and late developmental flowers indicating an intensive cell division in the gynoecial tissues during fruit development.

LAC2 is a multicopper-containing glycoprotein belonging to the laccase family. The exact function of LAC genes is not yet cleared. Laccase genes were highly expressed in lignifying tissues and are known to be associated with lignin synthesis (Gavnholt & Larsen, 2002; Sato et al., 2001). In Arabidopsis array experiments, the majority of investigated laccase members were found primarily in root tissues and in flower and stem samples secondly. In a mutant screening experiment, LAC8 and LAC2 deficient mutants lead to early flowering and reduced root elongation phenotypes. LAC2 expression was also detected in flower tissues in the Arabidopsis array (TAIR database) (Gavnholt & Larsen, 2002). In our Ambrosia flowering system, LAC2 was found to have been expressed uniquely in late 2F flowers; however, it was present with a reduced value of order compared to the other female-specific genes in this library.

The role of ACA7 in flowering pathways is not yet reported. It belongs to the Zn metalloenzymes carbonic anhydrases (CAs) that catalyze the reaction of reversible hydration of CO₂ with exceptionally high efficiency. CAs have been involved in a broad range of biochemical processes that involve carboxylation or decarboxylation reactions including photosynthesis and respiration. CAs also participate in pH regulation, inorganic carbon transport, ion transport, water, and electrolyte balance (Henry, 1996; Smith & Ferry, 2000; Tashian, 1989). However, an interesting correlation revealed, that ACA7 location was mapped on the same chromosome region as EARLY FLOWERING 3 (ELF3) including the earliness per se locus Eps-Am1 in diploid wheat (Triticum monococcum). This locus affects the duration of early developmental phases and it is responsible for the optimal photoperiod sensing and vernalization during flowering time (Gavnholt & Larsen, 2002). Due to ACA7 location it may be considered that carbon-based catalytic processes are strongly related to ACA7 in female flower development in Ambrosia artemisiifolia. These results suggest that female flower differentiation may differ mainly in regulatory elements of the flowering pathway genes such as altered signal transduction and TF activity.

Expression of male flower specific genes

Of the annotable genes, the mostly expressed sequences showed homology with the followings:

The homolog of desiccation-related protein, PCC13-62 (MK098104, GO:no terms) showed the highest RPKM value, 1,126.37, which is activated during the drought stress in the pollen, seeds, and vegetative organs in succulent plants (Giarola et al., 2018). It was expressed exclusively in the stylopodium of the hermaphrodite flower in Mucona sempervivens, where the nectary is located (Zha et al., 2013). Since the male flower pistillodium is a modified pistil with a function to help dispensing the pollen from staminate flowers in Ambrosia artemisiifolia (Payne, 1963), the anatomical presence of this transcript may be explained. Functionally, however, this protein may serve to protect the cytoplasmic contents in the pollen, since the pollen grain is an extremely desiccated structure at maturity. Desiccation-related proteins were also reported in the A. thaliana columbia ecotype pollen (Sheoran et al., 2006).

The second most significantly expressed sequence was the CYP450 86B1 (MK098105, GO:0004497) homolog. This gene was found to be expressed in many parts of the flower, including epidermis, anther, and pistil. In young inflorescences, this gene effects the suberin and cutin biosynthesis playing an important role in defensive mechanisms during biotic and abiotic stresses and supports the maintenance of reproduction (Compagnon et al., 2009; Duan & Schuler, 2005).

The third most significantly expressed gene homolog was found from the lipolytic enzyme family, GDSL2 (MK098106, GO:0016021), which was expressed during petal differentiation and is responsible for lipid catabolic processes (Ji et al., 2017). Ji et al. identified a GDSL lipase gene in Brassica rapa playing an important role in the anther and pollen development. Arabidopsis GDSL lipase 2 plays role in pathogen defense via negative regulation of auxin signaling and a positive correlation between late flowering and resistance to Fusarium oxysporum in Arabidopsis thaliana natural ecotypes was observed related to GDSL lipases (Lyons et al., 2015).

The above described male characteristic genes with outstanding high RPKM values refer to the most intensive physiological processes that are related to pollen exposure and stability and biotic and abiotic defense during anther formation.

Because all of the most male characteristic genes are related to flowering, secondary metabolism, plant cell structure, biotic and abiotic stresses, these may be considered to be essential to male specificity in the common ragweed. Gene homologs discussed below were found to be important for male flower physiology based on RPKM values. TET8 (MK098111, GO:0016021), which encodes a transmembrane protein tetraspanin 8, member of tetraspanin family. TET8 play important role in plant development, reproduction and stress responses (Reimann, Kost & Dettmer, 2017). Expression of tetraspanins including tetraspanin 8 were described during the pollen development and in mature pollen in Arabidopsis thaliana (Boavida et al., 2013; Honys & Twell, 2004; Pina et al., 2005). pMADS2 was found expressed in flowers in the second whorl, however in vegetative organs were not found in Petunia (Van Der Krol et al., 1993). In Ambrosia, this gene was found to be homologous to PI that is predominantly expressed in petals and stamens, and less in carpels and sepals. The 3-oxoacyl-(acyl-carrier-protein) synthase coding OAS (MK098107, GO:0006633) is responsible for the condensation reaction of fatty acid synthesis and confers resistance to low temperatures by maintaining chloroplast membrane integrity (Dayan & Duke, 2014). This protein is also involved in the regulation of fatty acid ratio during seed metabolism. It is required for embryo development, especially at the transition from globular to heart stage (Hakozaki et al., 2008). Another gene family responsible for membrane structure, the syntaxins (SYN) showed also unique expression in male tissues. These proteins play role in fusion of transport vesicles to target membranes. Inactivation of genes from two syntaxin families has led to lethality of the male gametophyte (Sanderfoot et al., 2001). NIP 4 1 and NIP 4 2 are pollen specific aquaporins playing role in pollen germination and pollination in A. thaliana. These genes are essential to the growth of pollen tubes and were considered to be exclusive components of the reproductive apparatus of angiosperms with partially redundant roles in pollen development and pollination (Di Giorgio et al., 2016).

MYB proteins are a superfamily of TFs that play regulatory roles of gene expression controlling direct organ development and defense responses in plants. We found uniquely expressed MYB transcription factors related to flower development in male tissues. MYB80 (MK098109, GO:0003677) showed the highest expression in male tissues playing an essential role during anther development and production of tapetum and microspores (Higginson, Li & Parish, 2003) (Phan et al., 2011). MYB26 (MK098115, GO:0003677) regulates lignified secondary cell wall thickening of the anther endotechium, which is necessary for anther dehiscence. It may play role in specifying early endothecial cell development by regulating a number of genes linked to secondary thickening (Mitsuda et al., 2006; Yang et al., 2007). MYB35 (MK098110, GO:0003677) is required for anther development and early tapetal function during microspore maturation. It regulates callose dissolution and is required for microspores release from the tetrads (Zhu et al., 2011). The MYB44 (MK098112, GO:0003677) represses the expression of protein phosphatases 2C in response to abscisic acid (ABA). It confers resistance to abiotic stresses dependent of ABA. The overexpression of AtMYB44 enhanced stomatal closure to confer abiotic stress tolerance in transgenic Arabidopsis (Jung et al., 2008).