For A. formosa, the most substantial expression difference occurred between DS3 and DS4, although the number of genes DE between DS4 and DS5 is also relatively high . This pattern more or less mirrors our sampling strategy in which finer sampling intervals were used between DS1 and DS4 with a larger developmental gap between DS4 and DS5. Given that there are relatively few changes across the earlier stages, for our first global analysis we focused on the end points of our sampling, bracketing Phase I of development, and identified the genes that increase or decrease in expression between DS1 and DS5. Within each species, between 5568 and 6933 genes were differentially regulated between DS1 and DS5. Counts of genes commonly up- and down-regulated between species are presented in Fig. 3a. A set of 1262 genes was commonly up-regulated across all species through development while a set of 1094 genes was commonly down-regulated across all species through development . Between 498 and 644 genes were uniquely up-regulated in just a single species across development while between 459 and 683 genes were uniquely down-regulated in a single species across development . Gene Ontology enrichment analyses were conducted for the genes commonly up- and down-regulated between DS1 and DS5 in all taxa. Genes related to mitotic activity, including DNA replication, mitotic chromosome condensation, and microtubule ontologies are enriched early in development, while genes belonging to oxidation reduction and fatty acid biosynthetic ontologies are enriched late in development . The sets of genes up- or down-regulated in all of the spurred taxa but not A. ecalcarata were also identified,indoor vertical farming with 521 genes commonly up-regulated in only the spurred taxa between DS1 and DS5 and 318 genes commonly down-regulated between these stages .
Although the petal is generally not considered a photosynthetic organ, GO enrichment analyses on the gene sets DE in only the spurred taxa showed an enrichment of genes with GO terms related to photosynthesis were up-regulated between DS1 and DS5 . Looking across all four species, 11,258 genes showed differential expression between DS1 and DS5 in at least one species, representing a full third of predicted genes in the Aquilegia genome . Conducting a principal component analysis of the expression levels of these genes across all samples showed that most of the variance can be explained by developmental time point . According to PC1 , the comparative developmental staging done across each species is fairly consistent with the expression data across taxa with a few exceptions. The A. ecalcarata samples assigned to DS1 have slightly higher PC1 values than the DS1 samples from the other taxa, suggesting that based on gene transcription, they may be more developmentally similar to the DS2 samples from the other taxa. Judging by PC1, the A. chrysantha samples assigned to DS4 and DS5 may be more transcriptionally similar to DS3 and DS4, respectively, of the other three taxa.In order to identify genes potentially involved in nectar spur development, differentially expressed genes between spurless A. ecalcarata and each spurred taxon were identified at each developmental stage. 15,588 genes are differentially expressed between A. ecalcarata and at least one of the spurred species during at least one developmental time point . Across the first three developmental time points, a greater number of the genes identified as DE between A. ecalcarata and each spurred taxon are expressed at a higher level in A. ecalcarata than the spurred species . By the fourth and fifth developmental stages, an approximately equal number of DE genes are expressed more highly in A. ecalcarata and each spurred species. DE gene sets for each of the spurred species versus A. ecalcarata were then compared to each other to identify genes that are commonly up- or down regulated in the spurred taxa versus A. ecalcarata at each developmental stage .
At each developmental stage, more genes are commonly up-regulated in A. ecalcarata relative to all of the spurred species than up-regulated in all of the spurred species relative to A. ecalcarata. The number of genes commonly up-regulated in spurred taxa gradually increases throughout development . The number of DE genes up-regulated in A. ecalcarata also increases across development, but the greatest increases occur between DS1 and DS2 , and between DS4 and DS5 . In total, 237 genes are commonly up-regulated in the spurred taxa relative to A. ecalcarata at all five developmental stages and 453 genes are commonly up-regulated in A. ecalcarata versus the spurred taxa at all five developmental stages . As a comparison, we also looked at the number of loci that were commonly up- or down-regulated in any one of the spurred species relative to the other three taxa . There was a common pattern in that there were more genes up-regulated in the focal species of the comparison than down-regulated , however, more genes were commonly DE when A. ecalcarata was the focal species. Conducting PCA on the set of genes that were commonly DE between all spurred taxa and A. ecalcarata at any developmental time point shows, as expected, that PC1 captures the difference between the spurred species and A. ecalcarata while PC2 captures common developmental differences across species . PC3 largely captures differences between A. sibirica and the other taxa while PC4 largely captures differences between the two North American species, A. formosa and A. chrysantha. While data points for PCs 3 and 4 generally cluster by species, interestingly, the A. ecalcarata DS5 samples cluster with A. formosa in PC3 and with A. chrysantha in PC4, although we have no hypothesis for why this may be. At each developmental stage, we tested for GO enrichment in the genes identified as either having higher expression in A. ecalcarata or commonly having higher expression in the spurred taxa versus a set of genes expressed in any species at that developmental stage . During the earlier developmental stages , genes related to heme/iron binding and oxidoreductase activity are over-represented in the A. ecalcarata datasets. By DS5, many GO categories related to mitosis are enriched in the set of genes identified as significantly up-regulated in the spurred taxa.
A prior study sought to identify genes important for nectar spur development in the horticultural variety A. coerulea ‘Origami’ by comparing gene expression between two regions of developing petals: the distal tip of the developing nectar spur and the petal blade. A. coerulea ‘Origami’ has relatively long nectar spurs, similar to those of hawkmoth-pollinated species, and these comparisons between developing spur and blade were made at two developmental times points,hydroponic vertical farming when petals were 1 mm and 3 mm long. As this prior study was conducted using an earlier Aquilegia genome annotation, we re-analyzed the data using the A. coerulea “Golsdmith” v3.1 genome annotation . Counts of genes DE either developmentally or tissue specifically are summarized in Supplemental Tables S11 and S12 . Focusing on DE genes between the blade and spur cup tissue, at the 1mm stage, 490 genes were more highly expressed in the blade and 280 genes were more highly expressed in the spur cup . At the 3 mm stage, 1178 genes were more highly expressed in the blade and 767 genes were more highly expressed in the spur cup . The results of this reanalysis are quite similar to the original analysis and exhibit a pattern common to that seen in the comparison of A. ecalcarata to the spurred taxa at early stages : more genes were commonly up-regulated at both the 1 mm and 3 mm stages in the spurless tissue sample than in the sample with spur tissue . These DE gene sets were compared to the lists of genes commonly DE between the spurred taxa and A. ecalcarata at roughly comparable developmental stages . Since A. ecalcarata petals consist primarily of blade tissue, we hypothesized that there would be more genes commonly up-regulated between whole A. ecalcarata petals and the A. coerulea ‘Origami’ blade tissue , and between entire petals of the spurred taxa and the A. coerulea ‘Origami’ spur tissue , compared to the opposite combinations . This is indeed the case, however, there are some genes that are more highly expressed in the blade tissue and the spurred taxa or in the spur tissue and in A. ecalcarata . For the purposes of identifying genes that make up a potential core module for spur development, we were particularly interested in genes that fall into the ‘blade’ and ‘spur’ classes defined above. Comparing genes that are more highly expressed in A. ecalcarata than the spurred taxa at all developmental stages to those that are expressed more highly in the blade than the spur at both the 1mm and 3mm time point reveals only 27 common genes . Looking at the intersect of genes that are more highly expressed in all of the spurred taxa at the stages we assessed with those that are expressed higher in the developing spur at both the 1 mm and 3 mm stages reveals an even smaller set of only 8 common genes . Of the 8 genes more highly expressed in ‘spur’ class, four show a general trend of decreasing through development in the spurred taxa,especially between DS4 and DS5 . These genes encode a dynein light chain protein, a myb/SANT-like transcription factor, a cytochrome P450 protein, and an oxysterol-binding protein. Two genes, encoding a leucine rich repeat protein and a heat shock factor transcription factor, show a pattern of increasing expression in the early stages assessed followed by a decrease, while another gene encoding a C2H2-type zinc finger transcription factor has more steady expression throughout development. One gene, encoding a xyloglucan endotransglycosylase protein, has relatively low expression in the early developmental stages with an uptick in expression moving into DS5 in the spurred taxa, especially in A. sibirica and A. formosa.
Interestingly, the STY homologs, which were among the most strongly differentially expressed loci in the blade to cup comparison and which were later found to be critical to nectary development, were not among the genes overlapping with the spurred to unspurred comparison. During later stages of petal development, however, there is some divergence between the STY homolog expression among the species sampled here, with A. ecalcarata showing lower expression than the other three . Among the 27 genes that are more highly expressed in the ‘blade’ class, several transcription factors stand out. These include a gene in the NAC transcription factor family with similarity to ANAC034 , an OVATE Family Protein with similarity to AtOFP4, and a gene with homology to the cryptochrome 2-interacting basic helixloop-helix transcription factor AtCIB1 .In Arabidopsis, many genetic factors have been identified that influence petal shape. These often function by regulating the timing of the transition from cell proliferation to cell differentiation . Genes identified as playing a role in controlling this transition can broadly be broken down into those that promote proliferation versus those that repress proliferation. The cell proliferation activators start off broadly expressed across the primordium but they become down-regulated in a basipetal wave that is complementary to an opposing pattern of up-regulation of the cell proliferation repressors. As the early phase of nectar spur development in Aquilegia involves prolonging localized cell divisions in the spur, we explored the expression of the Aquilegia homologs to these Arabidopsis candidates in our dataset to see if their expression patterns are generally consistent with a conserved role in Aquilegia petal development or differ between A. ecalcarata and the spurred taxa.Negative regulators of cell proliferation in A. thaliana petals include BIG BROTHER , DA1, the TCP genes TCP4 and TCP5, and the cyclin-dependent kinase inhibitor genes KIP RELATED PROTEIN 4and KRP2 . Although it may be expected that these genes would show the opposite pattern of expression from the cell-proliferation promoters, only a few of the Aquilegia homologs showed a pattern of increasing expression approaching DS5 . Two KRP-like genes, AqKRP2-like and AqKRP-like showed a strong increase in expression, especially between DS4 and DS5. Both of theses genes were also expressed more highly in A. ecalcarata than the spurred taxa. A third KRP-like gene with homology to KRP4 actually sharply decreased in expression during later developmental stages in all taxa assessed. Aside from the two KRP-like genes previously mentioned, AqDA1 is the only other candidate cell division repressor that appears to be increasing in expression toward DS5.