Study system
Loasaceae subfam. Loasoideae is a monophyletic and mostly Neotropical angiosperm subfamily. South Andean loasas diverged from the core Loasoideae during the mid-Eocene, around 45 Ma. Within South-Andean loasas, Blumenbachia Schrad., Loasa Adans., Pinnasa Weigend & R.H and Scyphanthus Sweet are entomophile, bee-pollinated genera. Divergence of these genera took place between the late Eocene, around 40 Ma., and the middle Oligocene, around 26 Ma. The South-Andean loasas genus Caiophora C. Presl. originated during the middle Oligocene, and several speciation events took place during the last 12 Ma. [6]. Remarkable speciation in the genus has been attributed to different Andean uplift events [6, 36] and to the acquisition of hummingbird (ornithophily) and rodent (therophily) pollination at high-elevation habitats [36]. Loasoideae flowers are protandric and their shape is highly complex. They present a corolla with separate pouch-shaped petals protecting the stamens and a whorl of androecium derived nectar scales [4]. Flowers are small and pendulous in bee-pollinated species and require the pollinator to land and hold onto the flower by grappling the nectar scales. Bee-pollinated flowers present open corollas, sometimes with highly reflexed petals, which make the nectar scales visible and easy to grasp (Fig. 1). Petal development in bee-pollinated species involves progressive expansion of the petal base [35, 37]. This shape change results in an open corolla with visible nectar scales. Expansion of the petal base is not completed in hummingbird-pollinated species and results in a narrow corolla [35, 37] (Fig. 1). Narrow corollas in hummingbird-pollinated species ensure that the body of the hovering hummingbird contacts the fertile flower structures [34]. Reversion from hummingbird- to bee-pollination took place at least once in Caiophora, and was paralleled by the evolution of small and open corollas with an expanded petal base [36] (Fig. 1).
This study includes the sampling of two species (Fig. 1). L. heterophylla originated around 15 Ma., is bee-pollinated and presents a small and open corolla composed of remarkably reflexed petals; C. hibiscifolia originated during the last 5 Ma., is pollinated by hummingbirds, and presents narrow and large corollas [6, 34]. The lineages giving rise to these two species diverged from a common ancestor around 40 Ma. [6]. We assumed that the pollination strategy and the corolla morphology of L. heterophylla represent the ancestral conditions in Loasoideae, based on maximum likelihood and Bayesian ancestral character state reconstructions (Additional file 1). L. heterophylla is distributed in low-elevation habitats in southern Chile; C. hibiscifolia is distributed in middle to high elevation habitats in the NW of Argentina [36].
Sampling
Sampling was conducted at the Botanical Gardens of Bonn University during the summer of 2013, except for two mature flower samples and a flower bud sample of L. heterophylla, which were collected during the summer of 2012. Our choice of species was constrained by availability of enough individual plants at the gardens. We sampled three outdoors individuals of L. heterophylla and three outdoors individuals of C. hibiscifolia. One mature male-phase flower and one flower bud with a diameter of 5 mm were sampled from each individual (this makes a total of 12 samples, including mature flowers and flower buds). One petal was collected from each flower and flower bud and stored in liquid air for RNA extraction. RNA sampling included an additional mature flower of L. heterophylla (making a total of 13 samples). An additional petal was collected from the mature flower and from the flower bud of one individual of each species and stored in liquid air for electron microscope imaging (a total of four samples).
SEM imaging
Scanning electron microscope (SEM) photos were obtained from the inner (concave) side of the petal, as the outer surface of the petal is covered with trichomes and the epidermal cells are not clearly visible. Before imaging, petals were critical-point-dried after fixation in 70% ethanol + 4% formaldehyde for at least 24 h and dehydrated with ethanol and acetone. After critical-point-drying (CP drying) petals were mounted as flat as possible on SEM holders. Dried samples were sputter-coated with a thin layer (< 30 nm) of palladium (Junker Edelmetalle, Waldbüttelbrunn, Germany). Scanning electron microscopy was performed with a LEO 1450 SEM (Zeiss, Jena, Germany). Images were recorded with a digital image acquisition system DISS 5 (Point Electronic, Halle, Germany). Petals were cut along the middle to reveal epidermal cells from the region adjacent to the petal midrib, and 8–9 pictures were obtained between the apical and basal end ends of the petal with a magnification of 100×. The first and the last picture in each sequence always corresponded to the apical and basal ends. Positions of the images along the petals were standardized across samples, so that a position at 0 corresponded to the basal end and a position of 1 corresponded to apical end. The SEM stage was tilted with the samples 30° to obtain good contrast, and the resulting distortion was compensated with the SEM-function ‘Tilt correction’.
Petal and cell geometry measurements
We randomly chose 10 cells from each SEM image and measured cell area (CA) in squared micrometers, and cell lobeyness. Cell lobeyness was calculated using a solidity index S ranging between 0 and 1. Cells with regular walls have S values that are closer to 1; cells with lobed walls have S values that are closer to 0. Cell lobeyness (CL) was calculated as 1 − S, so that cells with regular walls have CL values that are closer to 0, while cells with lobed walls have CL values closer to 1. We calculated the cell length-to-width ratio (CLWR) in the same longitudinal orientation of the petal, by measuring the maximum cell length in the longitudinal orientation of the petal and the maximum cell width in the perpendicular direction, and dividing the first by the second.
Pictures of the two halves of the cut petals used for SEM images were made to measure total petal area. Total petal area in squared micrometers was obtained by adding up the areas of the two cut petal halves. We obtained an estimate of the number of epidermal cells in the petals of the flower bud and the mature flower of each species by calculating the area of each SEM image, counting the number of epidermal cells in each SEM image, obtaining the average number of cells per squared micrometer in each SEM image and multiplying it by total petal area. We, therefore, have 8–9 calculations of cell number per petal, each corresponding to a different SEM image. We obtained a final average cell number per petal by averaging these calculations and also obtained the standard error of this estimated value. Cells were measured and counted with Image J [31].
Statistical analyses of CL, CA and CLWR along the middle rib of the petal
We used a linear model regression, including a quadratic term, on the relationship between cell lobeyness (CL) and the position along the petal midrib (PAMR) of the SEM image. We also used a linear model regression, including a quadratic term, on the relationship between the logarithm of cell area (CA), from here on log(CA), and PAMR. These models were used to test predictions about cell morphology along the petal midrib in L. heterophylla and C. hibiscifolia, under a scenario of delayed differentiation of epidermal cells in the later (Fig. 2). A linear model including a quadratic term was also used on the relationship between cell length-to-width ratio (CLWR) and PAMR. This last model was used to test whether differences in cell elongation at the petal base are related to interspecific differences in the expansion of the petal base. Since the relationship between CL, log(CA) and CLWR and PAMR was expected to differ between developmental stages and species, and we also expected the effect of the interaction between PAMR and developmental stage to depend on the species (Fig. 2), our models included a triple interaction among PAMR, developmental stage and species.
RNA preparation and sequencing
Frozen petal samples were pulverized with a mortar and pestle and total RNA was extracted and purified using the RNeasy Plant Mini Kit (Qiagen). Total RNA was checked for integrity using a BioAnalyzer with an Agilent RNA 6000 Nano Chip, following the manufacturer’s instructions (Agilent, https://www.agilent.com). We provided 5 μg of total cDNA from each sample to the Max Plank Genome Center Cologne for sequencing. cDNA was sequenced as 100 bp, paired-end reads using HiSeq2500 (Illumina). Libraries obtained for two out of the three L. heterophylla bud samples were sequenced later, at 150 bp. All data are available through NCBI’s Short Read Archive (BioProject accession: PRJNA763894).
De novo transcriptome assembly, read mapping and fragment counting
FastQ files for each species were quality trimmed with Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic), and fed to Trinity-v2.4.0 (https://github.com/trinityrnaseq/trinityrnaseq/wiki) to assemble species-specific reference transcriptomes. Completeness of each assembly was assessed using BUSCO (https://busco.ezlab.org) v4.1.3 with the embryophyta_odb10 lineage data set. Mapping of reads from each species to its assembled reference transcriptome and transcript abundance estimation was performed using the align_and_estimate_abundance.pl and abundance_estimates_to_matrix.pl perl scripts included in Trinity-v2.4.0 [16], using Bowtie2 [21] as alignment method and RSEM [22] as abundance estimation method.
Reference genome gene model annotation and transcriptome-to-genome mapping
Since there is no available genome from any species within the family Loasaceae, available genomes from the closest possible families were considered as a source for annotated “common” reference. The genome from the Chinese Happy Tree, Camptotheca acuminata (Nyssaceae) was chosen for this purpose (https://doi.org/10.1093/gigascience/gix065). The gene models were annotated using Trinotate (https://github.com/Trinotate/Trinotate.github.io/wiki). Details about gene model annotation can be found in the R Report of the RNAseq analysis (Additional file 2). L. heterophylla and C. hibiscifolia transcripts were mapped to C. acuminata gene models using BLASTx searches. Based on these mappings, a mapping table assigning a C. acuminata gene model to isoforms of each species that had a hit with an e value of 10e−3 or lower was built; all other isoforms were discarded.
Merging of count tables and exploratory analyses
Using the above mapping tables, a C. acuminata gene model was assigned to each isoform (row) of the read count table of each the two sequenced species. Reads for multiple isoforms matching a single C. acuminata gene model were collapsed by summing, while isoforms without a match were discarded. The resulting tables were then merged into a single table using C. acuminata gene model id as the key field. Using an inner join, only collections of transcripts with hits to C. acuminata gene models in both species were retained. Count data was normalized using a weighted trimmed mean of the log expression ratios (TMM, [29]), which normalizes by effective library size, but not by feature length. Then, counts were normalized to counts-per-million (cpm) using the cpm function (edgeR library). A filter was applied to remove any genes not having at least 1 cpm in at least three samples. We performed clustering and PCA exploratory analyses to check that samples clustered according to our sampling design, based on their RNA expression profile.
Differential gene expression analysis
Differentially expressed genes (DEGs) were detected using the DESeq2 (https://doi.org/10.1186/s13059-014-0550-8) approach, which estimates variance-mean dependence in count data and tests for DEGs using a model based on the negative binomial distribution. Four contrasts were created for DEGs: (1) L. heterophylla 5 mm bud vs. L. heterophylla mature flower; (2) C. hibiscifolia 5 mm bud vs. C. hibiscifolia mature flower; (3) L. heterophylla 5 mm bud vs. C. hibiscifolia 5 mm bud; (4) L. heterophylla mature flower vs. C. hibiscifolia 5 mm mature flower. The false discovery rate (FDR) threshold for significance was set to 0.05. Heatmaps of DEGs for L. heterophylla and C. hibiscifolia samples can be consulted in the R report (Additional file 2).
Gene set enrichment analysis (GSEA)
We conducted gene set enrichment analyses (GSEAs) for genes that are related to cell wall lobeyness in each of the four contrasts described above. Each GSEA requires a ranked list of the DEGs in the contrast. Were, therefore, ranked the differentially expressed genes in our four contrasts, based on their log2(fold-change). GSEA also requires a signature or an array of signatures that define a set of focal genes, cell lobeyness and cell elongation genes in our case. We created two signatures of cell wall lobeyness. The first includes genes that affect cell wall lobeyness trough interaction between turgor pressure and heterogeneous deposition of cellulose fibres and microtubule bundles on anticlinal cell walls, plus lobe outgrowth driven by actin filaments [27, 30]. We call this the ‘turgor pressure-cell wall interaction’ (TCWI) signature. The second signature includes genes that affect cell wall lobeyness through processes that are intrinsic to the cell wall and involve heterogeneous pectin demethylesterification of anticlinal cell walls [17]. We call this the ‘intrinsic cell wall property’ (ICWP) signature. The TCWI signature includes genes that code for kinesin-like proteins, actin-related proteins, Rac-like GTP-binding proteins, Rho of plants (ROP) proteins, CRIB domain-containing proteins, a cellulose synthase and a CLIP-associated protein. It also includes genes that code for the auxin-related proteins ABP1 and PIN. The ICWP signature includes genes that code for pectin methylesterases and pectinesterase inhibitors, as well as galacturonosyltransferases. We created an additional signature for genes that are related to differentiation and maturation of flower organs in the literature. These include genes coding for MADS-box proteins, for BEL1-like homeodomain proteins, for the Zinc finger protein JAGGED, for Auxin response factor ETTIN, for DELLA protein RGA, for protein SPOROCYTELESS and for TCP transcription factors [7, 10]. We call this the ‘flower differentiation-maturation’ (FDM) signature. We created a final signature for genes that are related to cell elongation in the literature. This included genes coding for DELLA proteins [41] and aquaporin genes [25, 41]. We call this the ‘cell elongation’ (CE) signature. To create the signatures, we first searched among the GO terms of the annotated C. acuminata reference transcriptome for protein and gene names in the protein families described above, as they appear in UniProt (https://www.uniprot.org/). The gene IDs of the filtered genes were kept and used in the GSEA. We did four GSEA analyses, one for each contrast. In each GSEA we looked for matches between the gene IDs of the TCWI, the ICWP, the FDM and the CE signatures and the gene IDs of the ranked DEGs list of the contrast. GSEA was done using the GSEA function of the clusterProfiler R library [43]. This function calculates a sum statistic for each signature and an enrichment score. It also permutes the rows of the ranked DEGs table, to calculate a null distribution and a P value for the enrichment score.