Linking gene expression to unilateral pollen-pistil reproductive barriers

Unilateral incompatibility (UI) is an asymmetric reproductive barrier that unidirectionally prevents gene flow between species and/or populations. UI is characterized by a compatible interaction between partners in one direction, but in the reciprocal cross fertilization fails, generally due to pollen tube rejection by the pistil. Although UI has long been observed in crosses between different species, the underlying molecular mechanisms are only beginning to be characterized. The wild tomato relative Solanum habrochaites provides a unique study system to investigate the molecular basis of this reproductive barrier, as populations within the species exhibit both interspecific and interpopulation UI. Here we used a transcriptomic approach to identify genes in both pollen and pistil tissues that may be probable key players in UI. We confirmed UI at the pollen-pistil level between a self-incompatible population and a self-compatible population of S. habrochaites. A comparison of gene expression between pollinated styles exhibiting the incompatibility response and unpollinated controls revealed only a small number of differentially expressed transcripts. Many more differences in transcript profiles were identified between UI-competent versus UI-compromised reproductive tissues. A number of intriguing candidate genes were highly differentially expressed, including a putative pollen arabinogalactan protein, a stylar Kunitz family protease inhibitor, and a stylar peptide hormone Rapid Alkalinization Factor. Our data also provide transcriptomic evidence that fundamental processes including reactive oxygen species signaling are likely key in UI pollen-pistil interactions between both populations and species. Our transcriptomic analysis highlighted specific genes, including those in ROS signaling pathways that warrant further study in investigations of UI. To our knowledge, this is the first report to identify candidate genes involved in unilateral barriers between populations of the same species.


Pollination type significantly influences a small number of genes
We expected that stylar transcriptomes would exhibit differences due to the type of 143 pollination and whether or not a pollination was compatible. We performed three separate 144 pairwise comparisons to look at differential gene expression in unpollinated (UP) LA1777 styles 145 compared to different pollination types (intrapopulation, self and interpopulation). Pairwise 146 comparisons between pollinated and UP LA1777 styles revealed only a small number of 147 differentially expressed genes, all of which were expressed at very low levels (<0.3 counts per 148 million (cpm)), Additional File 1, Table S1 (self vs UP); Additional File 1 Table S2 149 (intrapopulation vs UP) and Additional File 1, Table S3 (interpopulation vs UP)). In all 150 pollination treatments we identified pollination-induced differential expression of a number of 151 hypothetical proteins, RNA-directed DNA polymerases and retrovirus-related reverse 152 transcriptases. However, only three genes showed similar regulation among all pollination treatments: a hypothetical protein (Sopen05g010690) and two genes that have high homology to a Copia-like retrotransposon (Sopen01g019140 and Sopen02g008750). 155 Five genes were similarly differentially expressed in incompatible crosses with LA1777 156 as female (self and interpopulation crosses, Additional File 1, Table S1 and Table S3) versus 157 compatible crosses (intrapopulation, Additional File 1, Table S2), three of which encoded 158 hypothetical proteins. In addition, a cystathionine beta synthase (CBS) domain-containing 159 protein (Sopen06g019460) was downregulated 3-fold in incompatible crosses. It has been 160 proposed that CBS proteins are redox regulators involved in the modification of cell wall 161 composition, and in Arabidopsis, changes in CBS expression can reduce self-fertility [55]. A 162 homolog of Defective meristem silencing 3 (DMS3, Sopen03g019410), a gene that is involved in 163 silencing and epigenetic modification, was also downregulated 3-fold in incompatible crosses 164 [56]. 165 Thirty two genes that were differentially regulated between interpopulation pollination 166 and UP styles (Additional File 1, Table S3), but not in other treatments (self vs UP, Additional 167 File 1, Table S1 or intrapopulation vs UP, Additional File 1, Table S2).Interesting candidates 168 included a Ras-related GTPase (Sopen04g023040, Rab3), a pollen specific calcium binding 169 annexin (Sopen04g003390), and a K + transporter (Sopen11g006300) were upregulated ~3-fold in Intriguingly, in Pyrus pyrifolia and Papaver rhoeas, the disruption of ROS signaling in 177 incompatible (self) pollen tubes leads to depolymerization of the actin cytoskeleton and an 178 increase in nuclear DNA degradation [62,63]. 181 Although a number of potentially interesting candidates involved in UI were identified 182 from our pairwise comparisons, they were expressed at relatively low levels, often showed only 183 low fold-changes, and had relatively high p-values (p < 0.05). A principal components analysis 184 (PCA) of all stylar treatments showed that most of the variation between samples could be 185 explained by source population, and that there was no consistent grouping of styles due to 186 pollination treatment (Additional File 2, Fig. S1). The fact that few significant differences were 187 detected between UP styles and pollination treatments may reflect results of previous studies 188 indicating that the genes involved in UI-competence are expressed in styles regardless of 189 pollination status [31,37,54]. In other words, the UI-competent styles appear 'primed' to reject 190 incompatible pollen, as large changes in gene expression are not seen between UP styles and 191 those pollinated with compatible versus incompatible pollen parents. Because of these results, 192 we also performed an analysis in which all stylar treatments within a population (UP, pollinated 193 with self, intrapopulation and interpopulation pollen) were pooled to capture differential 194 expression between UI-competent (LA1777) versus UI-compromised (LA0407) styles. 195 We identified 179 genes that were significantly upregulated in UI-competent versus 196 compromised styles, and 179 that were significantly downregulated (Additional File 1, Table S4   197   and Table S5). However, we focused our interest predominantly on genes that showed a mean expression level over 3 normalized cpm and were upregulated over 10- compromised styles, we found ten that are involved in oxidation-reduction reactions (Table 1).

203
These included three Cytochrome P450 genes (Sopen02g037430, Sopen06g021780,  In UI-competent styles, we also identified three upregulated genes that are putatively 213 involved in defense response. These included an Endo-Chitinase involved in jasmonic acid and 214 ethylene signaling (Sopen02g027710) that was highly expressed (182 cpm) and upregulated over 215 20-fold, a Beta-Glucosidase (Sopen11g004500) and a Disease Resistance Protein 216 (Sopen10g025030) ( Table 1). These genes are of interest, as many molecular components 217 involved in plant-pathogen interactions show significant overlap with those involved in pollen 218 tube growth and guidance [64].

219
One gene of particular interest that was highly upregulated in UI-competent styles is the 220 hypothetical protein Sopen02g033850 (Table 1), which shows 80% amino acid identity to the peptide hormone Rapid ALkalinization Factor (RALF) from the wild potato, S. chacoense. 222 Small secreted peptides including RALFs are involved in a wide variety of plant functions 223 including development and immunity [65], and a pollen-specific RALF from S. lycopersicum 224 (SlPRALF) was found to negatively regulate pollen tube elongation in vitro [66]. We also 225 identified a protein involved in oligo-peptide transport (Sopen05g001950) which was 226 upregulated >28-fold in UI-competent styles and may be involved in the transport/secretion of 227 peptide hormones such as RALFs.

228
Another gene of interest that was highly upregulated in UI-competent versus 229 compromised styles was a Kunitz family protease inhibitor (Table 1). In Nicotiana, the Kunitz 230 family protease inhibitor, NaStep, is highly expressed in the pistils of SI species and is thought to 231 stabilize HT-proteins, although the mechanistic basis of this interaction remains unknown [67].

232
The NaStep protein is taken up by both compatible and incompatible pollen tubes, and the 233 transgenic suppression of NaStep in SI Nicotiana species compromises rejection of both self-and 234 some types of interspecific pollen tubes [67]. The reduced expression of this gene in UI-235 compromised LA0407 may reflect this population's lack of ability to reject some types of 236 interspecific pollen tubes [1,6].

237
Although we were most interested in genes showing the highest upregulation in UI-238 competent styles, we also analyzed genes that are highly downregulated. We found that three of 239 the top 25 most highly downregulated genes in UI-competent styles were involved in oxidation-240 reduction reactions, and one was putatively involved in defense response (Additional File 1, 241   Table S5). Notably, five of the top 25 genes downregulated in UI-competent styles are 242 putatively involved in cell wall modification (Additional File 1, Table S5). These include two (Sopen00g008620) and two Glycosyl Hydrolases (Sopen04g034210 and Sopen07g001240), all 245 of which were upregulated over 13-fold. Eight additional genes involved in cell wall 246 modification were downregulated in UI-competent styles, although to lesser extents (Additional 247 File 1, Table S5).

248
In addition to our analysis of genes that are highly up-or downregulated in UI-competent 249 versus compromised styles, we investigated the expression of 33 a priori candidates (Additional 250 File 1, Table S6) based on reports in the literature of stylar genes that may be involved in UI [26, 251 28, 31, 54]. These candidates included genes that were identified in a recent study comparing 252 stylar transcriptomes between UI-competent S. pennellii 0716 and the cultivated tomato S. 253 lycopersicum M82 which shows no UI response [54]. Only three of the 33 a priori candidates 254 showed over a 2-fold increase in styles of UI-competent LA1777 versus UI-compromised 255 LA0407 (Additional File 1, Table S6). These included a Pectin-Methylesterase Inhibitor 256 (Sopen04g027820, upregulated 2.4-fold) and two Glucosyltransferases (Sopen08g002330 and 257 Sopen08g002350, both upregulated over 8-fold); however all were expressed at low levels, < 1 258 cpm. HT-protein, a small asparagine-rich protein required for S-RNase-based UI [33] was 259 highly expressed in styles from both populations and was downregulated slightly in LA1777 260 (HT-A, 1.3-fold). Interestingly HT-B transcript was expressed in both populations, although 261 neither population accumulates functional protein due to an early stop codon in the transcript 262 [31]. Each individual of LA1777 harbored two unique S-RNase alleles as expected (data not 263 shown), whereas LA0407 did not express S-RNase transcript, as has been documented in 264 previous studies [31].

Differential gene expression related to UI-competence in pollen
As coordinated interactions between pollen and pistil are required for successful 268 fertilization, we also compared pollen transcriptomes between populations. Pollen tubes of 269 LA1777 reach ovaries of all SI and SC species within the tomato clade [1]. However, LA0407 270 pollen tubes are rejected by all SI species, including by SI S. habrochaites LA1777 [1]. The 271 inability of LA0407 pollen to traverse the styles of LA1777 and other SI Solanum species 272 suggests that it has lost a pollen factor(s) required for S-RNase resistance. For our analysis, we 273 therefore considered LA1777 pollen to be UI-competent, and LA0407 pollen to be UI-274 compromised. 275 We identified 90 genes that were upregulated in UI-competent pollen (Additional File 1, 276 Table S7) and 99 that were downregulated (Additional File 1, Table S8). Ten of the top 25 genes 277 upregulated in UI-competent versus UI-compromised pollen were annotated as hypothetical 278 proteins ( Table 2; Additional File 1, Table S7). Upon further analysis, the hypothetical gene  [68][69][70], and pistil AGPs are known to 282 stimulate pollen tube growth [69][70][71][72]. We also identified a Rab-GTPase (Rab4A, 283 Sopen01g033860) that was upregulated nearly 200-fold in UI-competent pollen. Pollen specific 284 proteins from this family have been found to promote pollen tube tip growth and play a role in 285 the ability of the pollen tube to sense directional cues [73,74]. 286 Our analysis of pollen also identified differentially expressed genes that may be involved 287 in protein degradation pathways ( Table 2; Additional File 1, Table S7 and Table S8). These pistil side factors such as S-RNases through the proteasomal degradation pathway [11, 15, 21, 291 26, 27, 29, 75]. We found two F-box/Skp-2 -like genes (Sopen05g003280 and 292 Sopen12g022970) that were upregulated over 150-fold in UI-competent pollen ( Table 2). 293 Another gene involved in protein degradation, an aspartyl protease (Sopen01g032940) was 294 highly upregulated 189-fold in UI-competent pollen ( Table 2). No previously identified SLFs 295 [28] were significantly differentially regulated, nor was the pollen UI factor Cullin1 [26] 296 (Additional File 1, Table S9). 297 We identified two protein kinases that were upregulated over 50-fold in UI-competent 298 pollen, one of which encodes a calcium binding serine/threonine wall-associated kinase 299 (Sopen10g028190), and the other a cysteine-rich receptor like kinase (RLK) (Sopen02g013470) 300 ( Table 2). Because RLKs have a proven role in pollen tube growth [76], these genes are of 301 potential interest. Further, pollen-expressed RLKs are involved in reception of peptide and 302 hormone signals from the pistil [77][78][79].

303
Genes involved in transcriptional regulation are likely to be important components of 304 pollen tube growth. We identified two types of transcription factors known to play key roles in 305 stress response [80,81] that were highly upregulated in UI-competent pollen: a NAC-domain 306 containing protein (Sopen03g020850) and a bZIP family protein (Sopen12g006170) that has 307 previously been localized to pollen ( Table 2; Additional File 1, Table S7).

308
An analysis of genes that were highly downregulated in UI-competent pollen identified 309 an F-box Protein (Sopen11g004020) and a Cysteine-rich RLK (Sopen05g014070) that were 310 downregulated over 45-and 100-fold respectively (Supplemental Table S8). Genes showing 311 over 100-fold decreases in UI-competent pollen also included two Histone Deacetylases 312 (HDACs, Sopen11g004040 and Sopen11g004050; Additional File 1, Table S8). The proper function of HDACs has been linked to successful pollen tube germination and tip growth in 314 Picea willsoni [82]. 317 Our initial analyses comparing transcriptomes of LA1777 and LA0407 reproductive 318 tissues led to some intriguing candidates that might be integral to UI. In an attempt to further 319 narrow down candidate gene and to expand our analysis to a broader range of S. habrochaites 320 populations, we performed a second RNA-seq experiment that included UP styles and pollen 321 from S. habrochaites populations with selected phenotypes (Table 3; Additional File 1, Table   322 S10). 323 First we confirmed that the UI phenotypes of the individuals tested reflected previous 324 results [6] (phenotypes summarized in Table 3). A PCA of all samples shows that much of the 325 sample variation (1 st PC) can be explained by tissue type, whereas a smaller percentage of 326 variation (2 nd PC) is explained by source population and potentially mating system (Fig. 3).

327
Interestingly, the additional populations selected for transcriptome analysis, which lie between 328 LA1777 and LA0407 geographically, cluster between these two populations in the PCA.  Table 3). Because styles of LA2119 are unique in that they accept 334 interpopulation LA0407 pollen while rejecting interspecific pollen, we interpret the information in Table 4 to reflect the variation in UI phenotype of LA2119 styles (i.e., genes that are up in 336 LA2119 are top candidates in interspecific UI, whereas genes that are down in LA2119 are top 337 candidates in interpopulation UI).

338
The top candidates for interspecific UI-competence in styles identified in the linear 339 discriminant analysis (LDA) included the hypothetical protein Sopen01g011020, that is highly 340 expressed and upregulated >100-fold (Table 4). Six genes potentially involved in oxidation- Oxidoreductase, Sopen12g006800).

346
Top candidates for interpopulation UI-competence in styles identified in the LDA 347 included two potentially involved in ROS signaling (Table 4). A DELLA-like transcription 348 factor (Sopen01g0311700) which is responsive to gibberellins [83] and participates in ROS 349 signaling by regulating ROS accumulation [84], was upregulated 11-fold. A Hemoglobin 350 Protein (Sopen08g021970) containing redox active transition metals was upregulated 3-fold.  The LDA of pollen genes was more straightforward in that LA2119 and LA1264 were 356 expected to be UI-competent (i.e., similar to LA1777) whereas LA1223 is UI-compromised and 357 may be missing pollen factor(s) required to traverse SI styles. As shown in Table 5, using the LDA we identified 22 genes that were upregulated in UI-competent pollen. Two of these encode  transporter was upregulated over 16-fold only in incompatible interpopulation pollinations.

395
These transporters generally pump reactive H 2 O 2 into the apoplast, an acidic environment with 396 low numbers of ROS scavengers, resulting in oxidative stress [85]. A K + channel and a Ca 2+ -397 binding annexin protein, both of which were annotated as being pollen-expressed were also 398 upregulated in interpopulation interactions, as was a Rab-GTPase. All three of these proteins can 399 play important roles in ROS generation and signaling. In pollen tubes, annexins may provide an 400 important link between Ca 2+ , the membrane and the cytoskeleton [91]. Further, many annexins 401 form Ca 2+ channels which are vital not only for the oscillating Ca 2+ influx associated with pollen 402 tube tip growth, but also for facilitating ROS signaling, cell elongation and cell wall remodeling 403 [85,92]. Interestingly, SKOR K + channels like the one identified in our analysis can also act as ROS-activated Ca 2+ channels [92]. Small GTPases, including Rabs, increase NADP(H)-oxidase 405 activity in a Ca 2+ -dependent manner [58,93], and their proper function is required for pollen 406 tube growth [58,60,73,74,94]. For example, the overexpression of both active and mutant 407 forms of the RAB11 protein leads to the inhibition of pollen tube growth in Nicotiana [74] 408 suggesting that the correct balance of multiple Rabs is required for effective pollen tube growth.

409
In UI-competent styles, ten of the top 25 most highly upregulated genes are putatively 410 involved in ROS generation and/or signaling, including an NADP(H)-oxidase and multiple 411 cytochrome P450s (Table 1). We also identified two highly upregulated gene candidates in the 412 flavonoid pathway, which could produce flavonoid compounds to act as pro-or anti-oxidants 413 ( Table 1). In our subsequent analysis using transcriptome data from additional S. habrochaites 414 populations, we found that six of the nine candidate stylar genes that may be involved in 415 interspecific UI were ROS pathway genes (Table 4). Surprisingly, only two ROS-linked genes 416 were upregulated in interpopulation UI-competent styles, one of which encodes a DELLA-like 417 transcription factor that inhibits ROS accumulation and restrains cell expansion [83,84]. In sum, 418 these results suggest a dynamic interplay between pollen and pistil that must be held in a tight 419 balance for pollen tubes to successfully grow through styles to reach the ovary.

420
The generation of ROS is linked to cell expansion, growth, cell wall cross linking and 421 callose deposition [95]. Pollen tube walls consist of a number of polymers (callose, cellulose 422 and pectin) that are highly cross-linked to each other; however the mechanisms by which cell 423 wall modification is regulated remains unclear [96][97][98]. Studies using microarray analysis have 424 found an upregulation of cell-wall modification genes in pollinated versus unpollinated styles 425 [51]. However, few have specifically investigated specificity to the SI or UI response (but see callose deposition between SI and UI crosses in S. peruvianum wherein SI crosses showed large 428 levels of callose deposition at the pollen tube tips, but interspecific crosses did not [99]. Our 429 pairwise comparisons in UP versus pollinated LA1777 styles did not reveal any genes involved 430 in cell wall modification. However, we did identify a number of differentially expressed genes 431 involved in cell wall modification in our larger analysis of stylar tissue (Additional File l, Table   432 S4 and Table S5), most of which were highly downregulated in UI-competent styles.

433
One of the most intriguing gene candidates from our analysis of stylar tissue included a 434 putative RALF peptide hormone (Sopen02g033850) that was upregulated over 22-fold in UI-435 competent styles. Peptide hormone signaling is involved in numerous processes during pollen-436 pistil interactions, from pollen hydration to fertilization [77]. A tomato RALF has been shown to 437 reduce pollen tube growth during specific windows of development [66], and therefore a stylar-  In UI-competent pollen, we identified a highly upregulated putative AGP 448 (Sopen12g014190) ( Table 2) that contains a signal peptide, a hydrophobic C-terminal domain, 449 eight dipeptides that are found in known AGPs and is composed of > 35% Pro/Ala/Ser/Thr (PAST) amino acids: a defining characteristic of AGPs (Additional File 2, Fig. S2; [70,101]). In 451 Arabidopsis, pollen AGPs are required for pollen tube growth, and may be involved in complex 452 signaling cascades [68][69][70], and pollen AGPs have been localized to pollen tube tips in some 453 species [102]. Another gene that warrants further study is a Rab-GTPase (Sopen01g033860) that 454 was upregulated nearly 200-fold in UI-competent pollen. It will be interesting to see if 455 increasing the expression of these genes in UI-compromised pollen is able to increase rates of 456 pollen tube growth through SI styles. demonstrates variability in mating system [41,44], as well as interspecific [31,35] and  Table 3 for more information.  Pollinations were covered with mesh bags to prevent unintended pollen deposition by 491 pollinators. Pollen tube growth was assessed using fluorescence microscopy as previously 492 described [31], and the length of styles and the point in the style at which no more than three pollen tubes passed were measured using ImageJ 1.47v (http://rsb.info.nih.gov/ij/; [104]). 496 The primary RNA-seq experiment was performed to identify genes involved in 497 interpopulation pollen tube rejection observed in crosses between LA1777 females and LA0407 498 males. Samples consisted of unpollinated styles, pollinated styles (self-pollination, 499 intrapopulation pollination, or interpopulation pollination), and pollen ( Fig. 1; Additional File 1, (including stigmas) were harvested directly into RNALater (Qiagen) and stored at 4º C for one 513 week, after which styles were blotted dry and immediately frozen at -80º C until processing.   528 Prior to mapping and assembly, reads were trimmed and filtered using the SHEAR 529 program (http://www.github.com/jbpease/shear; [54]). Briefly, SHEAR first uses the Scythe 530 algorithm (https://github.com/vsbuffalo/scythe; [105]) to remove adapters from the 3′ end, and 531 then filters low quality reads (mean Q<10), reads with >7 ambiguous bases (N's), reads < 50 bp, 532 and repetitive reads with mutual information score > 0.5. The program then trims reads on both 533 ends by removing low quality bases (Q<20), poly-A or poly-T runs of n ≥ 12, and ambiguous 534 bases. The appearance of AGATC at the 3′ end was also removed as we presumed it was an 535 adapter fragment. We removed both reads in a pair if either one of them failed the filters. On 536 average, 2.9% of all reads failed to pass the filter (min 2.5%, max 3.4%). The full command and 537 parameters used for SHEAR can be found in Additional File 2, Method S1. 538 We mapped RNA-seq reads to the Solanum pennellii reference genome using the STAR 539 spliced aligner with default parameters [106]. The reference genome sequence and the genome 540 annotation (Spenn v2.0) were downloaded from solgenomics.net [107]. On average across libraries, 81% of reads mapped uniquely to the reference genome. We counted reads mapped to 542 the reference gene annotation and to unannotated putative S-locus F-box (SLF) genes [28] (a 543 total of 48938 genic regions) using featureCounts v1.4.5-p1 [108]. A total of 663185500 read 544 pairs were counted (72% of raw reads). For genes annotated as 'hypothetical' in Spenn v2.0 545 [107], further sequence alignments were carried out using NCBI BLAST searches 546 (https://blast.ncbi.nlm.nih.gov; [109]). 549 For all tests of differential expression we used linear models implemented by the limma 550 package [110,111] and modules from the edgeR package [112] in R [113]. First, we normalized 551 and transformed the raw reads with voom, a weighted transformation based on the expected 552 relationship between expression mean and variance [111]. We computed t-statistics on the 553 transformed expression values for each gene using an empirical Bayes adjustment of standard 554 errors with the eBayes function [110]. 555 We initially searched for differential expression among stylar tissues by carrying out 556 separate pairwise comparisons. Specifically, we compared the following pairs of style treatments 557 in LA1777: intrapopulation-pollinated (compatible) against unpollinated (UP) styles, self-558 pollinated (incompatible) against UP styles and interpopulation (incompatible) against UP styles.

559
We visualized the genome-wide patterns of expression through a PCA of the normalized 560 mean read counts per gene (cpm reads in the library) with the prcomp function (implemented in 561 R; [113]). The PCA showed that there were negligible differences at the genome-wide scale 562 among all style treatments within a population ( Fig. 3; Additional File 2, Fig. S1). Because of 563 these high consistencies in their gene expression profiles, all style treatments (UP, self-, intra-and interpopulation pollinated) within a population were considered identical in our linear 565 models that focus on the differences between UI-competent (LA1777) and UI-compromised 566 (LA0407) styles. 567 We identified genes that are differentially expressed in UI-competent tissues using a 568 linear model with a single fixed effect (collinear with the population of origin) for pollen (n=6) 569 and styles (n = 24) separately. From these models, we took genes as differentially expressed if 570 they showed large differences in expression (> 3-fold change), with statistical significance at a 571 false discovery rate (FDR) of 5%. Further, we considered only tissue-specific genes: we required 572 genes to have a significant (FDR < 5%) tissue effect in a general linear model Y ~ P + T + e 573 (where P is a population effect, T is a tissue effect, and e is the error term). For stylar-side 574 factors, we included genes as differentially expressed if they were upregulated in styles with 575 respect to pollen, and vice versa for pollen-side factors. This last filter ensured that differences in 576 styles were unlikely to be contributed by pollen in the styles of pollinated samples.

577
A number of a priori pollen and pistil UI candidate genes were selected based on 578 information from previous publications [23,26,28,31,54]. Some of the SLF genes have not yet 579 been annotated in the S. pennellii genome, so we used locus numbers and sequences from Li and 580 Chetelat (2015) [28] to find these genes in our dataset. We identified expression levels of all a 581 priori candidates, and carried out statistical tests similar to the ones described above to determine 582 whether they were differentially expressed.    Gene annotation from Spenn_v2.0 was enhanced, when possible, for genes noted as hypothetical proteins (marked as [HP]).  A linear discriminant function was trained on the expression values of LA1777 and LA0407, and then used to classify other accessions as UI or non-UI.     Figure S1. The genome-wide patterns of expression in styles from two populations of Solanum habrochaites (LA1777 in orange and LA0407 in green), summarized by a principal components analysis. Most variance among samples is due to differences between populations (horizontal axis), and we found no large differences among experimental treatments at this scale (unpollinated styles, stars; self-pollinated, circles; intrapopulation-pollinated, triangles; interpopulation pollinated, squares).