Gene regulatory network architecture in different developmental contexts influences the genetic basis of morphological evolution

Convergent phenotypic evolution is often caused by recurrent changes at particular nodes in the underlying gene regulatory networks (GRNs). The genes at such evolutionary ‘hotspots’ are thought to maximally affect the phenotype with minimal pleiotropic consequences. This has led to the suggestion that if a GRN is understood in sufficient detail, the path of evolution may be predictable. The repeated evolutionary loss of larval trichomes among Drosophila species is caused by the loss of shavenbaby (svb) expression. svb is also required for development of leg trichomes, but the evolutionary gain of trichomes in the ‘naked valley’ on T2 femurs in Drosophila melanogaster is caused by reduced microRNA-92a (miR-92a) expression rather than changes in svb. We compared the expression and function of components between the larval and leg trichome GRNs to investigate why the genetic basis of trichome pattern evolution differs in these developmental contexts. We found key differences between the two networks in both the genes employed, and in the regulation and function of common genes. These differences in the GRNs reveal why mutations in svb are unlikely to contribute to leg trichome evolution and how instead miR-92a represents the key evolutionary switch in this context. Our work shows that variability in GRNs across different developmental contexts, as well as whether a morphological feature is lost versus gained, influence the nodes at which a GRN evolves to cause morphological change. Therefore, our findings have important implications for understanding the pathways and predictability of evolution.


Introduction
A major challenge in biology is to understand the relationship between genotype and phenotype, and how genetic changes modify development to generate phenotypic diversification. The genetic basis of many phenotypic differences within and among species have been identified [e.g. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15], and these findings support the generally accepted hypothesis that morphological evolution is predominantly caused by mutations affecting cis-regulatory modules of developmental genes [16]. Moreover, it has been found that changes in the same genes commonly underlie the convergent evolution of traits [reviewed in 17]. This suggests that there are evolutionary 'hotspots' in GRNs: changes at particular nodes are repeatedly used during evolution because of the role and position of the gene in the GRN, and hence the limited pleiotropic effect of the change [18][19][20][21].
The regulation of trichome patterning is an excellent system for studying the genetic basis of morphological evolution [22]. Trichomes are actin protrusions from epidermal cells that are overlaid by cuticle and form short, non-sensory, hair-like structures. They can be found on various parts of insect bodies during different life stages, and are thought to be involved in, for example, thermo-regulation, aerodynamics, oxygen retention in semi-aquatic insects, grooming, and larval locomotion [23][24][25][26][27] (Fig 1).
Regions of dorso-lateral larval trichomes have been independently lost at least four times among Drosophila species [39,40]. Recombination mapping and functional studies have shown that in all cases analysed, this phenotypic change is caused by changes in several svb enhancers, resulting in a loss of svb expression [6,10,[39][40][41][42]. The modular enhancers of svb are thought to allow the accumulation of mutations that facilitate the loss of larval trichomes in certain regions without deleterious pleiotropic consequences. It is thought that evolutionary changes in larval trichome patterns cannot be achieved by mutations in genes upstream of svb because of deleterious pleiotropic effects, while changes in individual svb target genes would when over expressed [38,43]. Therefore, one could predict that changes in adult trichome patterns are similarly achieved through changes in svb enhancers [20,21].
The trichome pattern on femurs of second legs also varies within and between Drosophila species [1,44] (Fig 1). In D. melanogaster, an area of trichome-free cuticle or 'naked valley' varies in size among strains from small to larger naked regions. Other species of the D. melanogaster species subgroup only exhibit larger naked valleys [1,44]. Therefore, trichomes have been gained at the expense of naked cuticle in some strains of D. melanogaster. Differences in naked valley size between species have previously been associated with differences in the expression of Ultrabithorax (Ubx), which represses the formation of leg trichomes [44]. However, genetic mapping experiments and expression analysis have shown that naked valley size variation among populations of D. melanogaster is caused by cis-regulatory changes in miR-92a [1]. This microRNA represses trichome formation by repressing the svb target gene shavenoid (sha), and D. melanogaster strains with small and large naked valleys exhibit weaker or stronger miR-92a expression, respectively, in developing femurs [1,45]. Therefore, while svb is thought to be a hotspot for the evolutionary loss of patches of larval trichomes, it does not appear to underlie the evolutionary gain of leg trichomes in D. melanogaster.
Differences in GRN architecture among developmental contexts may affect which nodes can evolve to facilitate phenotypic change in different tissues or developmental stages. In addition, an evolutionary gain or loss of a phenotype may also result from changes at different nodes in the underlying GRN, i.e. alteration of a particular gene may allow the loss of a trait but changes in the same gene may not necessarily result in the gain of the same trait. Therefore, a better understanding of the genetic basis of phenotypic change and evaluation of the predictability of evolution require characterising the expression and function of GRN components in different developmental contexts, and studying how the loss versus the gain of a trait is achieved.
Here we report our comparison of the regulation of trichome development in legs versus embryos. Our results reveal differences in expression and function of key components of the GRN between these two developmental contexts. These differences indicate that svb is likely unable to act as a switch for the gain of leg trichomes because it is already expressed throughout the legs in both naked and trichome-producing cells. Instead, regulation of sha by miR-92a appears to act as the switch between naked and trichome-producing cells in the leg. This shows that differences in GRNs between different developmental contexts can affect the pathway used by evolution to generate phenotypic change.

Differences in gene expression between leg and larval trichome development
The embryonic expression, regulation, and function of many genes involved in larval trichome formation is well understood [for example, see 29,[32][33][34][35][36][37][38]42,46] (Fig 1A). To characterise the regulation of leg trichome development better we first carried out RNA-Seq of T2 pupal legs between 20 and 28 hours after puparium formation (hAPF): the window when leg trichomes are specified [44] (S1-S6 Files). We tested if genes known to play a role during larval trichome formation are also expressed in our samples and used a cut-off of 1 fragment per kilobase per million (FPKM) reads mapped to determine if a gene is most likely expressed or not. Note that we chose not to compare the actual expression levels of trichome genes in our leg data sets with those of previously published expression data for embryos because it is difficult to interpret what any quantitative differences in overall expression levels between these two heterogeneous mixtures of cells might mean with respect to trichome development. We found that key genes known to be involved in larval trichome formation are expressed in legs. These include Ubx, SoxN, tal, svb, and sha, as well as key components of the Delta-Notch, Wnt and EGF signalling pathways (Fig 1, and S1 Table). However, expression of several genes known to regulate larval trichome development [29,33,36] is barely detectable in legs (i.e. below or around 1 FPKM). These include Dichaete, Arrowhead, and abrupt, which are also known to regulate svb expression during larval trichome development [34,42] (Fig 1 and S1 Table). Furthermore, the expression of 28 of the 163 known targets of svb in embryos [29,33] is barely detectable in our dataset (FPKM at or below 1) (S2 Table). In addition, 12 out of the 43 genes thought to be involved in larval trichome formation independently of svb [33,36] are also expressed at levels of less than 1 FPKM in legs (S3 Table). Therefore, our RNA-Seq data evidence key differences in both upstream and downstream components of the leg trichome GRN when comparing it to what is known for the embryonic GRN that specifies larval trichomes.
Our leg RNA-Seq data also allowed us to compare expression between strains of D. melanogaster with different sizes of naked valley: Oregon R (OreR) which has a small naked valley and ebony 4 , white ocelli 1 , rough 1 (eworo) which has a large naked valley (Fig 1). The size of the naked valley in these two strains is caused by differential expression of miR-92a [1]. We found that none of the known regulators of svb are differentially expressed between these two strains. In addition, we did not detect any significant differences in the expression of svb itself or most of its target genes including sha (S1, S2 and S3 Tables). However, we did find a trend towards higher expression of jing interacting gene regulatory 1 (jigr1) in the large naked valley strain eworo, although this difference is not significant after p value correction for false discovery rate (FDR) (S1 Table). Interestingly, miR-92a is co-expressed with jigr1 during neuroblast selfrenewal [47] and it is located in one of its introns. Therefore higher expression of miR-92a may be indirectly detectable in eworo (S1 Table). These results are consistent with miR-92amediated post-transcriptional regulation causing differences in naked valley size, and since this only occurs in a small proportion of leg cells, the effect on transcripts is likely to be difficult to detect using RNA-Seq.

miR-92a is sufficient to repress leg trichomes and acts downstream of Ubx
We next further examined the function of specific genes during leg trichome development compared to their roles in the formation of larval trichomes. It was previously shown that mutants of miR-92a have small naked valleys [48], which is consistent with the evolution of this locus underlying natural variation in naked valley size [1]. We confirmed these findings using a double mutant for miR-92a and its paralogue miR-92b [47], which exhibits an even smaller naked valley (Fig 2). Note that we did not detect any changes to the larval trichome pattern in these mutants compared to heterozygotes. We examined the morphology of the proximal leg trichomes gained from the loss of miR-92a compared to the trichomes found more distally. We found that the trichomes gained were indistinguishable from the other leg trichomes (S1 Fig). This suggests that all of the genes required to generate leg trichomes are already transcribed in naked valley cells, but that miR-92a must be sufficient to block their translation. Indeed, we found that the extra trichomes that develop in the naked valley in the absence of miR-92a are dependent on svb because in a svb mutant background no trichomes are gained after loss of miR-92a (Fig 2). Furthermore, these results also show that trichome repression by Ubx in the naked valley requires miR-92a because trichomes in the miR-92a mutant develop in the region where Ubx is expressed [44]. Thus, our data confirm that Ubx plays opposite roles in the larval and leg trichome GRNs: in embryos Ubx activates svb to generate larval trichomes [46], while we show that Ubx-mediated repression of leg trichomes [44,49] depends on miR-92a (Figs 1 and 2).

Regulation of svb during leg trichome patterning
The results above suggest that svb is expressed in the naked valley but is unable to induce the formation of trichomes because of the presence of miR-92a. To test this further we examined the expression of svb transcripts in pupal T2 legs using in situ hybridization. However, this method produced inconsistent results among legs and it was difficult to distinguish between signal and background in the femur. Therefore we examined the expression of a nuclear localising GFP inserted into a BAC containing the entire svb cis-regulatory region, which was previously shown to reliably capture the expression of this gene [43]. We detected GFP throughout T2 legs at 24 hAPF including in the proximal region of the posterior femur (S2 Fig). This indicates that svb is expressed in naked valley cells that do not produce trichomes as well as in more distal trichome-producing cells.
We next investigated the regulatory sequences responsible for svb expression in T2 legs. To do this we carried out ATAC-Seq [50,51] on chromatin from T2 legs during the window of 20 to 28 hAPF when leg trichomes are specified [44]. Embryonic expression of svb underlying larval trichomes is regulated by several enhancers spanning a region of approximately 90 kb upstream of the transcription start site of this gene [5,10] (Fig 3). Several of these larval enhancers also drive reporter gene expression during pupal development [43]. We observed that the embryonic enhancers DG3, E and 7 contained regions of open chromatin according to our T2 leg ATAC-Seq data. However, we found additional accessible chromatin regions that do not overlap with known svb embryonic enhancers (Fig 3).
Deletion of a region including the embryonic enhancers DG2 and DG3 [Df(X)svb 108 ] results in a reduction in the number of dorso-lateral larval trichomes when in a sensitized genetic background or at extreme temperatures [5]. Moreover, Preger-Ben Noon and colleagues [43] recently showed that this deletion, as well as a larger deletion that also removes embryonic enhancer A ([Df(X)svb 106 ], see Fig 3), results in the loss of trichomes on abdominal segment A5, specifically in males. We found several peaks of open chromatin in the regions covered by these two deficiencies in our ATAC-seq dataset (Fig 3) and therefore tested the effect of Df(X)svb 106 on leg trichome development. We found that deletion of this region and Evolution of gene regulatory networks consequently enhancers DG2, DG3, Z and A did not affect the size of the naked valley or the density of trichomes on the femur or other leg segments of flies raised at 17˚C, 25˚C, or 29˚C (compared to the parental lines) (S3 Fig). This suggests that while this region may contribute to svb expression in legs, its removal does not perturb the robustness of leg trichome patterning.
Next, to try to identify enhancer(s) responsible for leg expression, we employed all available GAL4 reporter lines for cis-regulatory regions of svb (S4 Table)  To further test whether the expression of any of these regions is consistent with a role in trichome formation, we used them to drive expression of the trichome repressor miR-92a and the trichome activator sha-ΔUTR [1]. Intriguingly, driving miR-92a under control of only one of the fragments (VT057077) caused the repression of trichomes on all legs ( Driving sha-ΔUTR with VT057077 is sufficient to induce trichome formation in the naked valley (Fig 3) and on the posterior T3 femur (S5 Fig). Driving sha-ΔUTR under control of any of the other nine regions did not produce any ectopic trichomes in the naked valley on T2 or on any other legs. These results indicate that a single enhancer, VT057077, is able to drive svb expression throughout the second leg in both regions which normally produce trichomes and in naked areas.

svb and sha differ in their capacities to induce trichomes in larvae and legs
It was previously shown that miR-92a inhibits leg trichome formation by repressing translation of the svb target sha [1]. However, sha mutants are still able to develop trichomes in larvae, albeit with abnormal morphology [29]. These data suggest that there are differences in the functionality of svb and sha in larval versus leg trichome formation, and therefore we next verified and tested the capacity of svb and sha to produce larval and leg trichomes.
As previously shown [38], ectopic expression of svb is sufficient to induce trichome formation on normally naked larval cuticle (Fig 4). However, we found that ectopic expression of sha in the same cells does not lead to the production of trichomes (Fig 4). svb is also required for posterior leg trichome production [43] (Fig 2 and S6 Fig), but over expression of svb in the naked valley does not produce ectopic trichomes (Fig 4). Over expression of sha on the other hand is sufficient to induce trichome development in the naked valley [1] (Fig 4). These results show that svb and sha differ in their capacities to generate trichomes in larvae versus legs.
Interestingly, we observed that the ectopic trichomes produced by expression of sha-ΔUTR in the naked valley are significantly shorter than those on the rest of the leg (S1 Fig). This suggests that although sha is able to induce trichome formation in these cells, other genes are also required for their normal morphology. We observed that another characterised svb-target gene, CG14395 [33], is also a high-ranking predicted target of miR-92a [52]: its 3'UTR contains two conserved complete 8-mers corresponding to the binding site for this microRNA. We found that CG14395 is also expressed in pupal second legs according to our leg RNA-Seq data (S2 Table) and furthermore that RNAi against this gene resulted in shorter leg trichomes (S7 Fig). Therefore it appears that miR-92a also represses CG14395 and potentially other svb target genes in addition to sha to block trichome formation.

Over expression of tal or ovoB can induce trichomes
Svb acts as a transcriptional repressor and requires cleavage by the proteasome to become a transcriptional activator. This cleavage is induced by small proteins encoded by the tal locus [30][31][32]. We therefore tested if svb is unable to promote trichome development in the naked valley because it is not activated in these cells. We found that expressing the constitutively active form ovoB or tal in naked leg cells is sufficient to induce trichome formation (Fig 4), which is consistent with loss of trichomes in tal mutant clones of leg cells (S6 Fig). Furthermore, it appears that tal, like svb, is expressed throughout the leg (S6 Fig). It follows that svb and tal are expressed in naked cells but are unable to induce trichome formation under normal conditions because of repression of sha, CG14395 and possibly other genes by miR-92a. We Evolution of gene regulatory networks hypothesise that over expression of tal on the other hand must be able to produce enough active Svb to result in an increase of sha transcription to overwhelm miR-92a repression.

The GRNs for larval and leg trichome patterning differ in composition and evolution
The causative genes and even nucleotide changes that underlie the evolution of an increasing number and range of phenotypic traits have been identified [17]. An important theme that has emerged from these studies is that the convergent evolution of traits is often explained by changes in the same genes-so called evolutionary 'hotspots' [17,53]. This suggests that the architecture of GRNs may influence or bias the genetic changes that underlie phenotypic changes [18,19,21]. However, relatively little is known about the genetic basis of changes in traits in different developmental contexts and when features are gained versus lost [18].
It was shown previously that changes in the enhancers of svb alone underlie the convergent evolution of the loss of larval trichomes, while the gain of leg trichomes in D. melanogaster is instead mainly explained by evolutionary changes in cis-regulatory regions of miR-92a [1,6,10,[39][40][41]. We investigated this further by comparing the GRNs involved in both developmental contexts and by examining the regulation and function of key genes.
Our results show that there are differences between the GRNs underlying the formation of larval and leg trichomes in terms of the expression of components and their functionality. These changes are found both in upstream genes of the GRN that help to determine where trichomes are made, and in downstream genes whose products are directly involved in trichome formation (Fig 1). The latter may also determine the differences in the fine-scale morphology of these structures on larval and leg cuticle (Fig 1) [29].
Furthermore, while the key evolutionary switch in embryos, the gene svb, is also necessary for trichome production on the posterior leg, over expression of this gene is not sufficient to produce leg trichomes in the naked proximal region of the T2 femur. This is because the leg trichome GRN employs miR-92a, which inhibits trichome production by blocking the translation of the svb target gene sha and probably other target genes including CG14395. In the legs of D. melanogaster, miR-92a therefore acts as the evolutionary switch for trichome production, and consequently the size of the naked valley depends on the expression of this gene (Fig 5) [1].
Interestingly, we observed that the ectopic trichomes produced by over expression of sha-ΔUTR in the naked valley are significantly shorter than those on the rest of the leg (S1 Fig). Therefore, while sha is able to induce trichome formation in these cells, other genes, including CG14395, are also required for normal trichome morphology. This suggests that GRNs may be able to co-opt regulators, in this case possibly miR-92a, that can act in trans to regulate existing components. Such changes can facilitate phenotypic evolution by phenocopying the effects of 'hotspot' genes in contexts where their evolution may be constrained. While trichomes can be lost as a result of the loss of svb expression but not loss of sha alone, interestingly, over expression of miR-92a is also able to suppress trichomes on other structures, including wings [1,45], presumably through repression of sha and other genes like CG14395.

Other genetic bases for the evolution of leg trichome patterns?
In contrast to larvae, it is unlikely that mutations in svb can lead to evolutionary changes in legs to gain trichomes and decrease the size of the naked valley. This is because this gene (and likely all the other genes necessary for trichome production) is already transcribed in naked valley cells. In addition, a single svb enhancer is able to drive expression throughout the legs including the naked valley. Although other enhancer regions of this gene are able to drive some expression in patches of leg cells, none of these is naked valley-specific. This suggests that evolutionary changes to svb enhancers would be unlikely to only affect expression of this gene in the naked valley. It remains possible that binding sites could evolve in this global leg enhancer to increase the Svb concentration specifically in naked valley cells. This could overcome miR-92a-mediated repression of trichomes similar to our experiments where tal and ovoB are over expressed in these cells, or when molecular sponges are used to phenocopy the loss of microRNAs [54]. However, this does not seem to have been the preferred evolutionary route in D. melanogaster [1] (Fig 5).
Our study also corroborates that Ubx represses leg trichomes [44] whereas it promotes larval trichome development through activation of svb [46]. Moreover, our results indicate that Evolution of gene regulatory networks Ubx acts upstream of miR-92a in legs because it is unable to repress leg trichomes in the absence of this microRNA. It is possible that Ubx even directly activates miR-92a since ChIPchip data indicate that there are Ubx binding sites within the jigr1/miR-92a locus [55]. Intriguingly, there is no naked valley in D. virilis, and Ubx does not appear to be expressed in the second legs of this species during trichome development [44] (Fig 5). However naked valleys are evident in other species of the virilis and montana groups and it would be interesting to determine if these differences were caused by changes in Ubx, miR-92a, or even other loci (Fig 5).

Evolutionary hotspots and developmental context
To the best of our knowledge, our study is the first to directly compare the expression and function of components of the GRNs underlying the formation of similar structures that have evolved in different developmental contexts. Our results show that the GRNs for trichome production in larval versus leg contexts retain a core set of genes but also exhibit differences in the components used and in their wiring. These differences likely reflect changes that accumulate in GRNs during processes such as co-option [e.g. 56] and developmental systems drift [57][58][59], although it remains possible that the changes have been selected for unknown reasons.
Importantly, we show that the differences in these GRNs may help to explain why they have evolved at different nodes to lead to the gain or loss of trichomes. This supports the suggestion that GRN architecture can influence the pathway of evolution and lead to hotspots for the convergent evolution of traits [17][18][19]21]. Indeed, such hotspots can also underlie phenotypic changes in different developmental contexts. For example, yellow underlies differences in abdominal pigmentation and wing spot pigmentation among Drosophila species [7,11,60,61]. However, we demonstrate that it cannot be assumed that evolutionary hotspots in one development context represent the nodes of evolution in a different context as a consequence of differences in GRN architecture.
Our findings also highlight that the genes that underlie the loss of features might not have the capacity to lead to the gain of the same feature. Therefore, while evolution may be predictable in a particular context, it is very important to consider different developmental contexts and whether a trait is lost versus gained. Indeed, even when we map the genetic basis of phenotypic change to the causative genes it is important to understand the changes in the context of the wider GRN to fully appreciate how the developmental program functions and evolves. Since evolution is thought to favour changes with low pleiotropy [19,[62][63][64][65], the effects of genetic changes underlying phenotypic change should be tested more widely during development. Such an approach recently revealed that svb enhancers underlying differences in larval trichomes are actually also used in other contexts [43]. Interestingly, miR-92a is employed in several roles, including self-renewal of neuroblasts [47], germline specification [48], and circadian rhythms [66]. It remains to be seen if the evolutionary changes in this microRNA underlying naked valley differences also have pleiotropic consequences, and therefore if natural variation in naked valley size is actually a pleiotropic outcome of selection on another aspect of miR-92a function.

Fly strains, husbandry and crosses
Fly strains used in this study are listed in S4 Table. GAL4 lines for analysis of svb expression and RNAi lines for analysis of CG14395 were obtained from the Vienna Drosophila Resource Center (VDRC) [67,68]. Flies were reared on standard food at 25˚C unless otherwise indicated.
Replacement of the P{lacW}l(3)S011041 element, which is inserted 5' of the tal gene, by a P {GaWB} transposable element was carried out by mobilization in omb-GAL4; +/CyO Δ2-3; l(3) S011041/TM3,Sb flies as described in [31]. Replacements were screened by following UAS-GFP expression in the progeny. The P{GaWB} element is inserted in the same nucleotide position as P{lacW}S011041. Clonal analysis of tal S18.1 and svbR9 alleles were performed as previously described [69].
A transgenic line that contains the cis-regulatory region of svb upstream of a GFP reporter (svbBAC-GFP) [43] was used to monitor svb expression. Legs of pupae were dissected 24 h hAPF, fixed and stained following the protocol of Halachmi et al. [70], using a chicken anti-GFP as primary antibody (Aves Labs, 1:250) and an anti-chicken as secondary (AlexaFluor 488, 1:400). Images were obtained on a confocal microscope with a 60X objective. SUM projections of the z-stacks were generated after background subtraction. A filter median implemented in ImageJ software [71] was applied. The proximal femur image was reconstructed from two SUM projections using Adobe Photoshop.

Measurement of trichome length
For trichome length measurements, T2 legs were dissected, mounted in Hoyer's medium/lactic acid 1:1 and imaged under a Zeiss Axioplan microscope using ProgRes MF cool camera (Jenaoptik, Germany). Trichomes on distal and proximal femurs were measured and analysed using ImageJ software [71]. Statistical analyses were done in R version 3.4.2 [72].

RNA-Seq
Pupae were collected within 1 hAPF and allowed to develop for another 20 to 28 h at 25˚C. Second legs were dissected in PBS from approximately 80 pupae per replicate and kept in RNAlater. RNA was isolated using phenol-chloroform extraction. This was done in three replicates for two different strains (e 4 ,wo 1 ,ro 1 and OregonR). Library preparation and sequencing (75 bp paired end) were carried out by Edinburgh Genomics. Reads were aligned to D. melanogaster genome version 6.12 [73] using TopHat 2.1.1. [74]. Transcripts were quantified using Cufflinks 2.2.1 and differential expression analysis conducted using Cuffdiff [75] (S1-S7 Files). Genes expressed below or around 1 FPKM were considered not expressed. Raw sequencing reads are deposited in the Gene Expression Omnibus with accession number GSE113240.
The reads have been deposited in the Gene Expression Omnibus with accession number GSE113240. CG14395-RNAi males) and siblings without knockdown effect were used as controls (Ctrl). (A) Trichomes developing after knockdown of CG14395 in the proximal femur are significantly shorter around the naked valley area than on the remaining femur (distal part) and on femurs of the controls (p < 0.001). Data are normally distributed (Shapiro-Wilk normality test). Tukey's multiple comparison test was used to test for significance. (B) After knockdown of CG14395 in the whole leg, trichomes are significantly shorter both around the naked valley area and on the remaining femur (p < 0.001 and p < 0.01). Note that some controls show significantly different trichome lengths. Data are not normally distributed (Shapiro-Wilk normality test). Kruskal-Wallis and pairwise comparisons using Wilcoxon rank sum test were used to test for significance. (PDF) S1 Table. Comparison of RNA expression levels of upstream trichome network genes for T2 legs at 24 hAPF from two strains with different naked valley size (e 4 , wo 1