Co-expression clustering across flower development identifies modules for diverse floral forms in Achimenes (Gesneriaceae)

Plant materials

Ten of the 26 currently recognized species in Achimenes were sampled: five from Clade 1 and five from Clade 2 (Fig. 2). These species were chosen based on their hypothesized close relationships and their diversity in floral form (Fig. 1; Roalson, Skog & Zimmer, 2003). On the basis of previous molecular studies (Roalson & Roberts, 2016), two species were chosen as outgroups to represent related lineages: Eucodonia verticillata and Gesneria cuneifolia (Figs. 1 and 2). All selected species were grown in standard greenhouse conditions at Washington State University, under 16-h days, 24–27 °C and 80–85% humidity.

Whole flower buds from two stages of flower development were sampled for each species: an immature bud stage (Bud stage; Fig. 1) and an intermediate stage before anthesis (D stage; Fig. 1). Our definitions of the floral developmental stages in Achimenes were adapted from the stages determined for Antirrhinum majus, which uses a morphological and temporal framework (Vincent & Coen, 2004). The Bud stage represents the smallest reproductive bud distinguishable from vegetative buds, having lengths between 1 and 2 mm (Fig. 1). The appearance is marked by the emergence of the floral primordia, no accumulation of anthocyanin pigments and no elongation or growth of corolla tissue cells. The D stage represents a halfway point between the Bud and a pre-anthesis flower, where the length is about half the length of a fully developed flower (Fig. 1). Its appearance is marked by the accumulation of anthocyanins in the corolla tissue, the cells have elongated in the corolla tube and trichomes have emerged. Corolla spurs are found only in A. grandiflora and A. patens, and the spur begins development during the D stage as an outward growth of the corolla in the opposite direction from the tube opening (Fig. 1).

Each stage was sampled with two or three biological replicates, contributing four or six samples for each species. The bud stage was chosen to represent a baseline level of gene expression when floral phenotypes are largely similar between species that could then be compared to the intermediate stage where floral phenotypes diverge and many important pathways involved in flower development and pigmentation show increased expression levels (Roberts & Roalson, 2017). All tissues were immediately frozen into liquid nitrogen and stored at −80 °C.

RNA extraction and RNA-seq library construction

Total RNA was extracted from frozen tissue using the RNEasy Plant Kit (Qiagen, Hilden, Germany) and ribosomal depleted RNA samples were prepared using the RiboMinus Plant Kit (Thermo Fisher Scientific, Waltham, MA, USA). Stranded libraries were prepared using the NEBNext Ultra Directional RNA Library Kit (New England Biolabs, Ipswich, MA, USA), barcoded, and pooled based on nanomolar concentrations. Library quality and quantity were assessed using the Qubit (dsDNA HS assay; Thermo Fisher Scientific, Waltham, MA, USA) and BioAnalyzer 2100 (Agilent, Santa Clara, CA, USA). The library pool was sequenced across four lanes of an Illumina HiSeq2500 for paired-end 101 bp reads at the Genomics Core Lab at Washington State University, Spokane. All new sequence data were deposited in the NCBI Sequence Read Archive (BioProject: PRJNA435759).

Single replicate libraries for four of the species in the current study (A. cettoana, A. erecta, A. misera and A. patens) were previously sequenced for the same two timepoints (Bud and D stages) (Roberts & Roalson, 2017). These samples were combined with the newly generated data to assembly de novo transcriptomes and perform co-expression clustering. This data was deposited in the NCBI Sequence Read Archive (BioProject: PRJNA340450).

Transcriptome assembly

Per base quality, adapter contamination, and per base sequence content of raw reads was assessed with FastQC tools (http://bioinformatics.babraham.ac.uk/projects/fastqc). One library (AT6) returned poor yields and low-quality reads and was excluded from further analyses. Trimmomatic v.0.36 (Bolger, Lohse & Usadel, 2014) was then used to remove contaminating adapter sequences and dynamically trim low-quality bases (ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 HEADCROP:13 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50). Trimmed reads derived from an individual species were combined and de novo assembled into contigs using Trinity (--SS_lib_type RF; Grabherr et al., 2011; Haas et al., 2013), generating a reference transcriptome for each species. Sequence redundancy was reduced with CDHIT-EST (-c 0.99 -n 5; Li & Godzik, 2006). Open reading frames were then predicted for each transcriptome assembly with TransDecoder (Haas et al., 2013) using a BLASTp search against the SwissProt database (www.uniprot.org) to maximize sensitivity. Finally, transcriptome assembly completeness was assessed using BUSCO (Simão et al., 2015) against a set of single-copy plant orthologs (embryophyte_odb10 dataset).

Orthogroup inference

To help facilitate orthogroup inference, transcriptome assemblies from flowers in Sinningia eumorpha and S. magnifica were downloaded and combined with our assemblies (https://doi.org/10.5061/dryad.4r5p1; Serrano-Serrano et al., 2017a). Sinningia is closely related to our focal lineage in the Neotropical Gesneriaceae (Roalson & Roberts, 2016). To identify orthogroup sequences across all 14 species, we used a modified version of the approach described in Yang & Smith (2014). All peptide sequences were clustered using an all-by-all BLASTp search performed using DIAMOND (--evalue 1e−6; Buchfink, Xie & Huson, 2015) and results were clustered using MCL v14-137 (-I 1.4; Enright, Van Dongen & Ouzounis, 2002). Homolog clusters were aligned using MAFFT v7.271 (--genafpair --maxiterate 1000; Katoh & Standley, 2013) and alignments were trimmed using Phyutility v2.7.1 (-clean 0.1; Smith & Dunn, 2008). Codon alignments were then produced for each cluster using PAL2NAL v14 (Suyama, Torrents & Bork, 2006) and maximum-likelihood trees were constructed using FastTree v2.1.8 (Price, Dehal & Arkin, 2010) using the GTR model. Trees were then trimmed to exclude branches >10 times longer than its sister or longer than 0.2 substitutions per site. Final homolog group alignments were created from the trimmed trees using MAFFT and used for final homolog tree inference using RAxML v8.2.9 (Stamatakis, 2014) and the GTRGAMMA model. Orthogroups were inferred from the final homolog trees using the rooted ingroup (RT) method of Yang & Smith (2014), using S. eumorpha and S. magnifica as outgroups. Homolog trees were trimmed to include at least 2 ingroup species, resulting in 83,595 orthogroups.

Quantifying gene expression

Trimmed reads for each sample were aligned to the corresponding transcriptome assembly to quantify expression using default parameters in Kallisto v0.43.0 (Bray et al., 2016). The resulting counts for each contig were matched with the inferred orthogroups derived above, creating two separate matrices for the Bud and D stage samples. In the resulting gene expression matrices, rows corresponded to orthogroups and columns corresponded to individual samples. Counts for the Bud and D stage samples were separately transformed using variance-stabilizing transformation implemented in the Bioconductor package DESeq2 (Love, Huber & Anders, 2014) and quantile normalized using the Bioconductor package preprocessCore (Bolstad, 2018). Lastly, counts for both stages were separately corrected for confounding effects, such as batch or species effects, using principal component regression with the “sva_network” function in the Bioconductor package sva (Parsana et al., 2019).

Co-expression networks and identification of modules

The Bud and D stage RNAseq data was clustered into gene co-expression networks using the R package WGCNA (Langfelder & Horvath, 2008). Prior to network construction, orthogroups were filtered to remove any with >50% missing values or <0.3 expression variance. In total, 9,503 orthogroups common to both Bud and D stage data were used to construct separate networks for each stage. Parameters for the Bud stage network were as follows: power = 18, networkType = “signed”, corType = “bicor”, maxPOutliers = 0.05, TOMType = “signed”, deepSplit = 2, mergeCutHeight = 0.25. For the D stage network, the following parameters were used: power = 16, networkType = “signed”, corType = “bicor”, maxPOutliers = 0.05, TOMType = “signed”, deepSplit = 2, mergeCutHeight = 0.25. Soft-thresholding powers for each network (“power” parameter) were chosen as the lowest value such that the scale free topology model fit (R²) was ≥0.9. Signed networks consider positively correlated nodes connected, while unsigned networks consider both positively and negatively correlated nodes connected (Van Dam et al., 2018). We used signed networks here (“networkType” and “TOMType” parameters) because this type is considered more robust to biological functions with more specific expression patterns (Mason et al., 2009; Song, Langfelder & Horvath, 2012). Unsigned networks can capture only strong correlations and many negative regulatory relationships in biological systems are weak or moderate (Ritchie et al., 2009). Biweight midcorrelation (“corType” parameter) was used as this method is often more robust than Pearson correlation and more powerful than Spearman correlation (Song, Langfelder & Horvath, 2012). All other WGCNA parameters were kept at their default values.

Module eigengenes and orthogroup connectivity were calculated in each network separately using the “moduleEigengenes” and “intramodularConnectivity” functions in WGCNA, respectively. The module eigengene is defined as the first principal component of a module and represents the gene expression profile of a sample within a module (Langfelder & Horvath, 2008). Connectivity refers to the sum of connection strengths a node has to other network nodes.

Module hub genes (defined as the most highly connected nodes within a module) were identified in both networks based on module membership (kME) scores >0.9, calculated using the “signedKME” function of WGCNA. kME is defined as the correlation between the expression levels of a gene and the module eigengene (Horvath & Dong, 2008). Additionally, using connectivity measures, we defined peripheral genes in each network as the lowest 10% of connected nodes.

We tested whether Bud stage modules were preserved in the D stage and vice versa using module preservation statistics in WGCNA (Langfelder et al., 2011). Median Rank and Zsummary statistics were computed using the “modulePreservation” function of WGCNA, using 200 permutations (Langfelder et al., 2011). Median rank statistics were compared to the “gold” module, which consists of 1,000 randomly selected orthogroups. Modules with a lower median rank exhibit stronger preservation than modules with a higher median rank. Zsummary scores >10 indicate strong evidence of preservation and scores <2 indicate no evidence of preservation.

Phylogeny

We previously inferred a phylogeny for the 12 sampled species using 1,306 single-copy orthologs identified from the same transcriptome dataset used here (Roberts & Roalson, 2018). For comparative analyses of module-trait relationships, we randomly sampled 50 single-copy ortholog gene trees and rescaled branch lengths to be proportional to time (ultrametric) using the “chronos” function in the R package ape (Paradis, Claude & Strimmer, 2004). Bootstrap support was 100 for nearly every branch in the Roberts & Roalson (2018) phylogeny, therefore we chose to use randomly sampled ortholog trees to account for phylogenetic uncertainty.

Module and trait relationship

We correlated external traits (e.g., red flower color) with the module eigengenes using a phylogenetic Markov Chain Monte Carlo method, implemented in the R packages MCMCglmm (Hadfield, 2010) and mulTree (Guillerme & Healy, 2014). These analyses were performed individually for each module in both the Bud and D stage networks. Three floral traits and pollination syndrome were coded: primary flower color (purple, red, white and yellow), flower shape (funnelform, salverform and tubular), corolla spur (absent and present), and syndrome (bee, butterfly and hummingbird). We fit a multivariate mixed model with a non-informative prior, where the shape and scale parameters equal 0.001, and residual variance was fixed. Analyses were performed over the 50 randomly selected ortholog trees to account for phylogenetic uncertainty. A phylogenetic model was used to account for any interspecific non-dependency in the dataset by using a random effects structure that incorporated a phylogenetic tree model. Floral traits were set as categorical response variables and the module eigengenes were set as the predictor variable. For each tree, two MCMC chains were run for 250 k generations and discarding the first 50 k as burn-in. This process was run individually over each module in both the Bud and D stages. We checked for convergence between model chains using the Gelman–Rubin statistic (Gelman & Rubin, 1992). Potential scale reduction values were all less than 1.1 and effective sample sizes for all fixed effects were greater than 400. We considered fixed effects to be statistically significant when the probabilities in the 95% credible region did not include zero (Figs. S3–S6).

Evolutionary rate analysis

Codon alignments of the filtered orthogroups (n = 9,503) were produced using PAL2NAL v14 (Suyama, Torrents & Bork, 2006). Subsequent alignments and corresponding gene trees were used as input to the codeml package in PAML v4.9 (Yang, 2007) to estimate d_N/d_S (omega, ω). d_N/d_S for each orthogroup was estimated using a one-ratio model (model = 0, NSsites = 0), providing a single estimate for each orthogroup to match other single orthogroup metrics, such as connectivity. If d_N/d_S values were exceptionally large (=999) because of zero synonymous differences, those values were removed from the dataset.

First, we performed linear regression (LM) to test the effect of orthogroup connectivity, average expression levels, and their interactions (explanatory variables) on the estimated d_N/d_S values (main effect), LM = d_N/d_S ~ connectivity * expression. This was performed using the lm function in R with 10,000 bootstrap pseudoreplicates. Second, we ran two-sample t-tests with 10,000 permutations to test whether hub nodes and periphery nodes in each network had lower (for hubs) or higher (for periphery) d_N/d_S values than all background nodes. Third, we fit a LM with 10,000 bootstrap pseudoreplicates to test whether modules associated with a pollination syndrome (bee, butterfly, or hummingbird; see “Results” below) had increased d_N/d_S values, LM = d_N/d_S ~ syndrome. For all statistical analyses, d_N/d_S, connectivity, and expression values were log transformed to reach a normal distribution before being processed. Analyses were performed separately for each network using the separately estimated connectivity and expression levels.

Functional annotation and gene ontology enrichment analysis

Filtered orthogroups (n = 9,503) were annotated based on a search against all green plant proteins in the SwissProt database (release 2018_07; The UniProt Consortium, 2017) using Diamond BLASTp (—evalue 1e−6). Results were used to assign annotations and gene ontology (GO) terms. GO enrichment analyses were performed using the R package TopGO and the default “weight01” algorithm with the recommended cutoff of p < 0.05 (Alexa, Rahnenführer & Lengauer, 2006). Enriched GO term redundancy was removed and summarized using the semantic similarity measures implemented in the REVIGO web server (Supek et al., 2011).

Data and code availability

Raw sequencing reads are available from the NCBI Sequence Read Archive under BioProjects: PRJNA435759 and PRJNA340450. Data and code for all analyses can found at Zenodo: DOI 10.5281/zenodo.3517231.

Transcriptome assembly

Twelve transcriptomes were sequenced (72 libraries) and assembled in 12 species after sampling across two floral developmental stages (Table S1). The transcriptomes had an average assembly length of 244.9 Mb (±33.9 Mb), an average of 217,510 (±47 k) transcripts, and an average N50 length of 1,988 bp (±306 bp). Between 22,942 and 33,501 putative genes were detected with 2.34 (±0.35) putative isoforms. BUSCO identified an average of 87% of the complete, conserved single-copy plant orthologs (Table S1). Our focus in our comparative study was on the presence of conserved orthologs, and not on the presence of taxonomically restricted genes, which are not likely preserved across species.

Constructing the co-expression network for 12 species

We identified 83,595 orthologous clusters (orthogroups), each containing at least two species. After removing 74,092 orthogroups with too many missing values (>50% missing) or low expression variance (<0.3 variance) and keeping common orthogroups shared between the Bud and D stages, 9,503 orthogroups were used to construct two co-expression networks using the weighted gene correlation network analysis (WGCNA) approach. Modules of co-expressed orthogroups are inferred using the expression profiles of each sample regardless of species. Clustering identified 65 and 62 co-expression modules in the Bud and D stage networks, respectively (Fig. 3; Table S2). Module size in the Bud stage ranged from 22 (module ME63) to 592 (ME1) orthogroups (mean = 146) and ranged from 21 (ME61) to 704 (ME1) orthogroups (mean = 153) in the D stage. Out of the 9,503 orthogroups, 7,845 (83%) and 8,735 (92%) were assigned to modules in each network.

The remaining unassigned genes from each network that were not placed into a well-defined module were assigned to the ME0 module by WGCNA. Despite being unplaced, these orthogroups may still be functionally relevant within each network. Unassigned genes in the ME0 module were enriched for many processes, including transcription (Bud, n = 173 and D, n = 84) and organ development (D, n = 55), among other functions (Tables S3 and S4).

We tested whether modules in the Bud stage were preserved in the D stage, and vice versa. Overall, preservation of modules between the Bud and D stage networks was high. Eight of the 65 Bud stage modules (ME46, ME62, ME44, ME57, ME42, ME41, ME55 and ME40) were not strongly preserved in the D stage, while only 1 of the 62 D stage modules (ME59) was not strongly preserved in the Bud stage (Table S5). These patterns suggest that these non-preserved modules may play unique roles within their network. Among the eight Bud stage-specific modules were many orthogroups related to primary cell growth, meristem identity and differentiation, stamen development and flowering (Table S3). In the single D stage-specific module were many orthogroups related to methylation, flowering, and carotenoid biosynthesis (Table S4). Despite the relatively high preservation, correspondence between module membership in the Bud and D stages was moderate (Fig. 4). Orthogroup module membership tended to shift between each network. For example, most members of module ME15 in the Bud stage are transferred to modules ME3 and ME7 in the D stage (Fig. 4).

We defined network hubs (the most highly connected nodes within a module) with module membership (kME) scores >0.9. Using this definition, all modules in each network contained at least one hub, with the percentage of hubs ranging from 1% to 24% of orthogroups in a module (Figs. 5 and 6; Table S6). Of the 632 and 691 hubs identified in the Bud and D stages, a significant portion (110) were shared between networks (Fisher’s Exact Test, odds ratio = 3.02, P < 0.001). Functional enrichment of hubs in the Bud stage revealed that many were involved in the regulation of development and metabolic processes (n = 40), transport and protein localization (n = 39), biosynthetic processes (n = 19), and DNA methylation (n = 7), consistent with the notion that hubs usually play important roles in integrating other genes during network evolution (Ravasz et al., 2002) (Table 2). In the D stage, hubs were enriched for functions related to translation (n = 33), transport (n = 9), ethylene and abscisic acid signaling pathways (n = 17), carotenoid and anthocyanin biosynthetic processes (n = 11), and the regulation of gene expression (n = 20), among others (Table 2).

Nodes in the periphery of each network were loosely defined by considering them to be the lowest 10% connected (connectivity_Bud < 2.824, connectivity_D < 4.651). These definitions identified 951 periphery nodes (Table S7), of which 301 were peripheral in both networks (Fisher’s Exact Test, odds ratio = 5.67, P < 0.001). These periphery nodes in the Bud stage were enriched for functions related to reproductive structure development (n = 75), transport and localization (n = 107), and signaling (n = 389), among others (Table 3). In the D stage, the periphery was enriched for functions related to the regulation of flower development (n = 64), the regulation of flowering time (n = 18), the regulation of flavonoid biosynthesis (n = 7), and cell fate specification (n = 12), among others (Table 3). Greater than expected numbers of peripheral nodes in the D stage were related to regulation (Fisher’s exact test, odds ratio = 1.191, p = 0.02752) and transcription (Fisher’s exact test, odds ratio = 1.370, p = 0.001549). In contrast, there were no differences in the expected number of peripheral nodes related to regulation (Fisher’s exact test, odds ratio = 0.955, p = 0.7044) or transcription (Fisher’s exact test, odds ratio = 1.047, p = 0.3549) in the Bud stage.

Modules associated with floral traits

We calculated module eigengenes, single values which represent the gene expression profiles of a sample within a module (Figs. S1 and S2). The extent of module involvement in various biological processes can be tested by correlating eigengenes with external traits, such as flower color. We tested the relationship between the module eigengenes and the four floral traits (pollination syndrome (bee, butterfly, hummingbird), flower color (purple, red, yellow and white), flower shape (funnelform, salverform and tubular), and corolla spurs (absence and presence)), while controlling for any phylogenetic bias in the dataset using a phylogenetic mixed model (Figs. 5 and 6; Figs. S3–S6). We also examined whether eigengene values (>0.1 or <−0.1) in single species contributed to a strong correlation or whether eigengene values (>0.1 or <−0.1) in multiple species contributed to a strong correlation (Figs. 5 and 6; Figs. S3–S6). We considered there to be strong evidence for trait association when a significant correlation from MCMCglmm was coupled with elevated eigengene (>1.0) or decreased eigengene (<−1.0) expression in multiple species that shared any of the significant traits. Likewise, weak evidence for trait association was considered when the MCMCglmm results were significant, but only a single species or no species sharing the significant traits had elevated (>1.0) or decreased (<−1.0) eigengene expression. WGCNA does not use trait information when inferring modules, so modules that correlated strongly with any floral trait were inferred without a priori knowledge.

Evolutionary rates correlate to network location

We ran several analyses to explore how several network variables may affect evolutionary rates of the orthogroups within each network (as measured by d_N/d_S). First, we tested the effects of orthogroup connectivity, expression levels, and their interaction on the d_N/d_S values (Table S7). Results indicated a collective effect between connectivity and expression on d_N/d_S in both the Bud stage (LM, F_{(3, 9,548)} = 92.58, p = 0, ${\eta }^2$ = 0.03) and the D stage networks (LM, F_{(3, 9,548)} = 78.36, p = 0, ${\eta }^2$ = 0.02). We found that orthogroup connectivity was negatively correlated with d_N/d_S in the Bud stage network (LM, t = −16.31, p = 0) and the D stage network (LM, t = −15.13, p = 0) (Fig. 7), while expression was positively correlated with d_N/d_S in the Bud stage (LM, t = 2.12, p = 0.0344) and not in the D stage (LM, t = 1.51, p = 0.1324) (Fig. S7). The interaction between connectivity and expression also indicated a marginal dependency of both predictors on d_N/d_S in the Bud stage (LM, t = −2.022, p = 0.0433) but not in the D stage network (LM, t = −1.823, p = 0.0683). These patterns have been observed in model species and may be universal to all eukaryotes.

Second, we tested whether the nodes we defined as hubs in each network had lower d_N/d_S values, suggestive of increased evolutionary constraint (Table S9). Hubs are considered to functionally significant and may control large portions of the network. We found that the hubs in the Bud stage network had lower d_N/d_S than all background genes (Two-sample t-test, t_(679.31) = 7.0934, p = 1.65e−12, Cohen’s d = 0.383; Permutation t-test. mean diff. = 0.266, p = 0.0001) (Fig. S8). Similarly, hubs in the D stage network also had lower d_N/d_S than all background genes (Two-sample t-test, t_(755.91) = 4.1918, p = 1.55e−05, Cohen’s d = 0.206; Permutation t-test, mean diff. = 0.143, p = 0.0001) (Fig. S8).

Third, we tested whether the most peripheral nodes (the lowest 10% connected nodes) had higher d_N/d_S values, which may suggest relaxed evolutionary constraint (Table S9). These nodes are loosely connected within the network and their roles may be more apt to fluctuate during evolution and development. These peripheral genes had higher d_N/d_S than all background genes in the Bud stage network (Two-sample t-test, t_(1,487.70) = −10.272, p = 0, Cohen’s d = −0.253; Permutation t-test, mean diff. = −0.176, p = 0.0001) and the D stage network (Two-sample, t_(1,512.30) = −10.387, p = 0, Cohen’s d = −0.252; Permutation t-test, mean diff. = −0.175, p = 0.0001) (Fig. S8).

Lastly, we tested whether modules correlated with bee, butterfly, or hummingbird pollination syndromes showed increased evolutionary rates over the non-associated modules. Results indicated a very weak effect of syndrome association on d_N/d_S in the Bud network (LM, F_{(3, 9,548)} = 5.186, p = 0.001406, ${\eta }^2$ = 0.0016) and the D network (LM, F_{(3, 9,548)} = 4.277, p = 0.005038, ${\eta }^2$ = 0.0013) (Table S9). Only bee syndrome associated modules had marginally lower d_N/d_S values in both networks (Bud stage, LM, t = −3.58, p = 0.000346; D stage, LM, t = −2.508, p = 0.0122). Butterfly associated modules had marginally lower d_N/d_S only in the D stage (LM, t = −2.548, p = 0.0108), while hummingbird associated modules did not have any effect in either network.

Hub nodes contain important candidates for floral form

To date, a large number of studies conducted in diverse angiosperm systems, have found that many genes are associated with flower development (Quattrocchio et al., 1999; Howarth et al., 2011). Despite evidence for extensive gene/genome duplications that could expand the genetic repertoire (Tang et al., 2008), many recent studies have suggested that the generation of novel genes is rare (Paterson et al., 2009; De Smet et al., 2013). Thus, the evolution and development of diverse floral forms may predominantly occur through the co-option of existing pathways (Preston, Hileman & Cubas, 2011). To our knowledge, gene co-expression network analyses have not been previously used across multiple species in a non-model lineage to understand flower development and evolution.

Hub nodes (genes) are considered to be conserved, highly connected nodes that are central to network architecture and may directly orchestrate numerous biological processes (Ravasz et al., 2002; Hu et al., 2016; Mähler et al., 2017). Our analyses identified over 600 gene nodes as hubs within each network, the majority of which were only hubs within a specific floral stage. These hub nodes tended to have lower d_N/d_S estimates, suggesting that they not prone to rapid evolutionary change (Table S9). Unsurprisingly, hub nodes in the Bud stage were enriched for many processes related to early flower development, including cell growth, leaf formation, and photosystem assembly (Table 2). Similarly, hub nodes in the D stage were enriched for genes that may be related the development of distinct floral forms, including carotenoid biosynthesis, anthocyanin biosynthesis, and cell morphogenesis (Table 2). Within this context of flower development, we identified a number of hubs that we consider important candidates for involvement in the production of flower color and the development of floral form. Many of these candidates were considered previously (Roberts & Roalson, 2017) as potentially important for the diversification of floral form in Achimenes. Undoubtedly, this list of hubs contains many genes that may important for other biological traits that we do not consider here, such as scent, nectar, and pollen production.

Peripheral genes are enriched for transcriptional regulators in later stages of development

Peripheral nodes have low connectivity within co-expression networks (Table S7), tend to have higher rates of molecular evolution (Masalia, Bewick & Burke, 2017), and have been hypothesized to contribute more to evolutionary innovations than network hubs (Ichihashi et al., 2014). Of the 951 nodes we defined as peripheral in both Bud and D stage networks, roughly a third were consistently peripheral under our definitions (lowest 10% connected; Table S7). After identifying the nodes with homologies to transcriptional regulators, we found there was no difference in the expected number in the early Bud stage, while there were greater than expected numbers of transcriptional regulators in the periphery of the D stage (see “Results”). No such patterns existed when we tested these transcriptional regulators as hub genes. While Mähler et al. (2017) found an enrichment for transcription factors among hub genes in a single species network of Populus tremula, our results in a multi-species network find the opposite pattern. Shifting the expression patterns of various genes and pathways can lead to phenotypic divergence among closely related species (West-Eberhard, 1989). This can be seen during later stages of flower development when distinct traits, such as red or purple flower color, can appear in different species due to shifting expression patterns within the anthocyanin biosynthetic pathway (Smith & Rausher, 2011). The connection between the apparent overabundance of transcriptional regulators in later stages of flower development and the rapid phenotypic divergence in floral form will need to be explored further.

Several notable peripheral nodes were identified for their potential involvement in important floral processes (Table S7). Three homologs were identified in the periphery of the Bud stage with potential involvement in the determination of zygmorphic (bilateral) floral symmetry: two homologs of divaricata (DIV; cluster13684_1rr.inclade1.ortho1.tre, cluster8585_1rr.inclade1.ortho1.tre; ME26 and ME0) and one homolog of dichotoma (DICH; cluster13245_1rr.inclade1.ortho1.tre; ME0). Bilateral symmetry has been proposed as a mechanism that facilitates pollen transfer by insects (Galliot, Stuurman & Kuhlemeier, 2006). Both DIV and DICH are TCP transcription factors that were first characterized in Antirrhinum (Almeida, Rocheta & Galego, 1997; Luo et al., 1996), and show localized expression to the ventral and dorsal portions of the petal, respectively. Achimenes flowers at the Bud stage are lacking or only just beginning to establish floral asymmetry and the relegation of both DIV and DICH to the network periphery may reflect this small role. The role of these genes has not been studied in the Gesneriaceae, but have been examined more recently in other non-model systems for their role in the evolution of pollination syndromes (Zhang, Kramer & Davis, 2010; Preston, Martinez & Hileman, 2011). Additionally, a homolog of the homeobox wuschel gene (WUS; cluster9273_1rr.inclade1.ortho1.tre; module ME0) was identified from the periphery of the D stage network. WUS was first characterized in Arabidopsis for its role in meristem cell maintenance (Laux et al., 1996) and its expression is localized only to the meristem (Mayer et al., 1998). Once the floral meristem has been initiated, WUS starts to be repressed by AG (hub gene, described above) (Sun et al., 2009; Lohmann et al., 2001). The relegation of WUS to the periphery in our networks at the D stage might then be expected given that it plays such a prominent role during early floral development and the vegetative to reproductive transition.

Evolutionary rates are correlated to network connectivity

Increased protein evolutionary rates during rapid diversification has been suggested in other radiations of animals and plants (Kapralov, Votintseva & Filatov, 2013; Brawand et al., 2014; Pease et al., 2016). Our study is among the first to examine patterns of network evolution across multiple closely-related species of plants with multiple transitions in flower type and pollination syndrome. Hummingbirds are hypothesized to have had an important influence on diversification in Neotropical Gesneriaceae (Roalson & Roberts, 2016; Serrano-Serrano et al., 2017b), while the role of butterflies has not been explored. Modules correlated with hummingbird or butterfly pollination syndromes did not exhibit clear increases or decreases in evolutionary rates compared to other modules (Table S9). Detecting selection is subject to many factors, such as a gene’s transcriptional abundance and its importance within the protein interaction network (Lemos et al., 2005). For instance, we found that higher network connectivity was strongly associated with lower evolutionary rates across both networks in our study (Fig. 7; Table S9). This corroborates the recent results found in other eukaryotic systems (Morandin et al., 2016; Masalia, Bewick & Burke, 2017; Josephs et al., 2017; Mähler et al., 2017) while also showing the apparent strength of correlation between connectivity and evolutionary rate differs between studies and systems.

Network effects (such as the type of network analysis employed) and the expression levels of highly abundant genes are known to affect evolutionary rates, potentially confounding analysis of selection pressure (Krylov et al., 2003). We found that d_N/d_S in early development may be dependent on both connectivity and expression levels, while this interactive effect was gone during later development (Table S9). As there are few highly connected genes (hubs, which are important determinants of the observed network structure), a random mutation would be less likely to affect such a gene. Changing the amino acid sequence of a hub gene could have multiple associated effects, some deleterious (Hahn & Kern, 2005; Luisi et al., 2015). In comparison, a random mutation would be more likely to affect a gene with lower connectivity, while also resulting in variants with fewer negative consequences. The hub genes have more direct interactions and are more likely to be essential and involved in more biological processes, placing them under greater constraint.

Numerous network modules correlated to floral form

One primary use for gene co-expression network analyses is to identify modules of co-expressed genes whose overall expression (as measured by the eigengene) may correlate with different traits. These traits can be both qualitative or quantitative, such as social caste (Morandin et al., 2016), seed oil production (Hu et al., 2016), or anther development (Hollender et al., 2014). While we tested only qualitative traits here, future studies could also correlate eigengene expression to quantitative variables such as nectary size or sucrose concentration, since these appear to be another important component of pollination syndrome delimitation (Perret et al., 2001; Katzer, Wessinger & Hileman, 2019; Vandelook et al., 2019). Since our dataset was composed of samples from 12 species, we performed correlation analyses using a phylogenetic model to test whether module expression correlated to different flower colors, flower shapes, corolla spurs and pollination syndromes. Nearly half of our modules in both networks had evidence for an association with these floral traits (Figs. 5 and 6). There were no modules whose correlation to a pollination syndrome had increased eigengene expression in all species sharing that trait (Figs. 5 and 6). Instead many modules with correlation to a trait showed a phylogenetic pattern with a shared increased eigengene expression among the most closely related species. For example, module ME23 in the Bud stage network had a correlation to butterfly pollination and increased eigengene expression in A. cettoana and A. longiflora, both members of Clade 1 (Fig. 2). Our aim was to identify sets of co-expressed genes that would provide candidates for involvement in the development of these traits.

For instance, many of the candidate hubs we identified (see above) were found in modules that were correlated to one or multiple traits to which they might contribute, such as phantastica (PHAN; cluster11220_1rr.inclade1.ortho3.tre; module ME4, correlated to corolla spurs) and DFR (module ME48, correlated to hummingbird pollination and purple flowers) (Table S2) (Waites et al., 1998). Additionally, modules correlated to various floral traits were enriched for multiple genes that might be involved in the development of those traits (Table S8). Module ME53 was correlated to the corolla spurs found in A. grandiflora and A. patens in the D stage network and was enriched for many genes involved in cell division and cell growth (Fig. 6; Table S4). Other times, modules were correlated to a specific trait but the enriched GO terms did not provide a clear idea of how these co-expressed orthogroups may contribute (Tables S3 and S4). In our analyses, we considered relatively few traits and only focused on those that have traditionally been considered the most important for the differentiation of floral form in Achimenes, namely flower color, flower shape, and corolla spurs. Undoubtedly there are a myriad number of molecular and biochemical pathways that we have not considered here that underlie the development of these complex floral phenotypes. With the idea that networks represent complex systems of gene and protein interactions, it may be that these seemingly simple and generalized floral phenotypes are produced through much more complex means than many have previously considered.

Strengths and limitations of network analyses in non-model organisms

Our analysis covers genes with orthologs found in most species, representing a set of evolutionarily conserved genes. While no reference genome currently exists for Neotropical Gesneriaceae, we constructed, annotated, and presented de novo transcriptomes for 10 Achimenes species as well as two outgroup species, E. verticillata and G. cuneifolia. Our approach enables functional genomic studies in this non-model system, but the potential for misassemblies may bias our downstream results (Hornett & Wheat, 2012). Using a gene tree-based approach to carefully identify orthogroups, we focused our attention on conserved gene regulatory machinery and not on taxonomically restricted genes. The latter may be involved in species-specific functions. Nevertheless, it is worth noting that many recent studies have indicated that transitions in floral form may be due to the shifting expression of genes found in conserved genetic pathways (Wessinger & Rausher, 2015). Taxonomically restricted genes are likely interacting with existing regulatory pathways, and the manner in which they integrate into the co-expression modules identified here will be a fascinating topic for future research. With extended sampling within Achimenes we could infer species-specific co-expression networks and perform more detailed cross-species comparisons. We expect that these taxonomically restricted genes may be poorly connected with the network, which may allow for increased diversification.

We characterized gene expression patterns by sequencing RNA across two timepoints during flower development. Previous gene expression experiments across three timepoints (early Bud, intermediate D, and late Pre-anthesis stages) in four Achimenes species suggested that the majority of differential gene expression occurred between the early Bud stage and the intermediate D stage (Roberts & Roalson, 2017). We believe our approach allows us to capture dynamic transcriptional patterns during flower development for comparison across multiple species. We found that the co-expression networks for the Bud and D stages were largely preserved, despite there being a few stage specific modules (Table S5). However, this approach may greatly over-simplify the nature of transcriptional regulation in these species. Organ- and structure-specific expression studies may have greater power to detect floral form-specific differences, both in terms of the number of differentially expressed genes, and co-expression network structure (Hollender et al., 2014; Suzuki et al., 2017; Shahan et al., 2018). As a result, our data likely represents an underestimate of the true number of modules. Future studies could focus on the comparative analysis of expression patterns in specific organs, such as the petals, stigma, and stamens, and additional timepoints, which would provide greater functional insight into the evolution, development, and diversification of these organs.

Network analysis represents a more complex approach than differential expression analyses because it can capture system-level properties (Langfelder & Horvath, 2008). Most phenotypes involve interactions of proteins from diverse biochemical pathways. Although inferred modules do not necessarily correspond to entire biochemical pathways, or other components of cellular organization, the approach performs well in reconstructing the complexity of protein-protein interaction networks (Allen et al., 2012). In many systems this approach has succeeded in identifying candidate genes for various biological traits, for example floral organ development in strawberry (Hollender et al., 2014; Shahan et al., 2018) or carotenoid accumulation in carrots (Iorizzo et al., 2016). Our results demonstrate the difficulty in identifying gene clusters that underlie highly correlated floral traits, especially in non-model systems, as seen by the overlapping correlation of modules associated with flower shape and color (Figs. 5 and 6). It has been shown in many plant systems that floral traits tend to be selected together by pollinators (Cuartas-Domínguez & Medel, 2010; Sletvold, Grindeland & Ågren, 2010), making it challenging to separate correlated modules from causal modules. Co-expression networks are foundationally about correlations between genes (Van Dam et al., 2018), making them useful in identifying novel sets of genes that may be functionally relevant for trait of interest. Other approaches, such as quantitative trait locus (QTL) analyses may provide more promise in identifying distinct genes that are involved in the divergence of floral traits. QTL studies have already shown promise in identifying genetic loci that underlie pollination syndrome divergence in Mimulus (Yuan et al., 2013), Penstemon (Wessinger, Hileman & Rausher, 2014) and Rhytidophyllum (Alexandre et al., 2015).

We were able to identify numerous candidate orthogroups involved in the development of floral form in our non-model system. While no functional genomic data exists yet for Achimenes, our annotations are based on the manually curated SwissProt database. De novo assemblies provide important data for non-model organisms, but there can be biases introduced during assembly and annotation due to potentially missing genes (Hornett & Wheat, 2012). We attempted to limit this bias by focusing only on shared orthogroups found in the majority (at least six) of our sampled Achimenes. Applying this network-based approach to the diversification of flowers across multiple species in Achimenes, we recovered distinct modules of co-expressed genes that may underlie a range of floral phenotypic traits tied to different pollination syndromes. Although our data come from whole flowers during development, and do not allow us to test specific hypotheses regarding whether particular genes have shifted expression location in the developing flower, they do suggest that there are conserved regulatory modules that may have been co-opted in convergent flower types multiple times. It will be interesting to apply network analyses to study floral diversification in additional lineages and examine how conserved or divergent the patterns are across angiosperms.