Distinct Regulatory Changes Underlying Differential Expression of TCP Genes Associated with 32 Petal Variations in Zygomorphic Flowers of Petrocosmea (Gesneriaceae)

91 floral zygomorphy mainly through specific expression in dorsal petals of a flower. Their 93 expression changes are usually related to morphological diversity of zygomorphic flowers. However, it is still a challenge to elucidate the molecular mechanism underlying their expression differentiation. It is also unknown whether CINCINNATA ( CIN )-like TCP genes, locally 96 controlling cell growth and proliferation, are involved in the evolution of floral zygomorphy. To 97 address these questions, we selected two closely related species, i.e. Petrocosmea glabristoma and 98 P. sinensis with distinct petal morphology, to conduct expression, hybridization, mutant, and 99 allele-specific expression (ASE) analyses. The results show that the size change of the dorsal 100 petals between the two species is mainly mediated by the expression differentiation of CYC1C and 101 CYC1D , while the shape variation of all petals is related to the expression change of CIN1 . In 102 reciprocal F1 hybrids, the expression of CYC1C , CYC1D , and CIN1 conforms to an additive 103 inheritance mode, consistent with the petal phenotypes of hybrids. Through ASE analyses, we find 104 that the expression differentiation of these TCP genes is underlain by distinctly different types of 105 regulatory changes. We suggest that highly redundant paralogs with identical expression patterns 106 and interspecific expression differentiation may be controlled by remarkably different regulatory 107 pathways because natural selection may favor different regulatory modifications rather than 108 coding sequence changes of key developmental genes in generating morphological diversity. 109 110 111 112 113 114 115 116 117 118 119 whether different regulatory pathways or their interplay account the expression differentiation floral symmetry related genes therefore the diversification dorsal zygomorphic flowers. In this study, two Petrocosmea species, i.e. P. glabristoma and P. sinensis as the representatives of two clades, were selected to investigate possible regulatory pathways driving the morphological divergence of zygomorphic flowers by employing multiple approaches, including comparative expression analyses of CYC - and CIN -like TCP genes, hybridization, mutant, and ASE assays. The results showed that the expression differentiation of CYC1C , CYC1D and CIN1 coordinately promoted the morphological divergence of the dorsal petals that is characteristic of different types of zygomorphic flowers in the two species. In addition, both CYC - and CIN -like genes conformed to an additive inheritance mode in reciprocal F1 hybrids, 193 consistent with the phenotypes of the hybrids. We further found that cis -regulatory changes, 194 trans -regulatory changes, and a complex of both changes contributed to the expression differentiation of CYC1C , CYC1D , and CIN1 , respectively. This report will shed significant light on how regulatory changes underlie expression differentiation of key developmental genes and therefore morphological variations of the dorsal petals in zygomorphic flowers. Real-time PCR was also performed to investigate the expression patterns of all CIN -like 267 genes in all floral organs of the second whorl in both parents and F1 hybrids. The results showed 268 that CIN1 expression levels were similar in the dorsal, lateral, and ventral petals in both P.


INTRODUCTION 121
During the evolution of both animals and plants, a common phenomenon is that 122 spatial-temporal expression modifications of key developmental regulators, as a linkage between 123 genotype and phenotype, usually contribute to the origin of morphological novelties (Doebley et 124 al., 1997;Cong et al., 2002;Gompel et al., 2005;Prud'homme et al., 2006Prud'homme et al., , 2007Werner et al., 125 et al., 2003;Palatnik et al., 2003;Crawford et al., 2004;Koyama et al., 2007;Nag et al., 2009). 151 Therefore, it would be interesting to decipher whether CIN-like gene expression changes bring 152 about morphological transitions in non-model organisms. It is also an important issue whether 153 CYC-and CIN-like TCP genes act together to control specific morphologies and generate 154 successive variation forms during the evolution of zygomorphic groups. 155 Generally, gene expression can be regulated at transcriptional, posttranscriptional, and 156 translational levels. Of these, transcriptional regulation determines when and where a gene is 157 expressed, at what level, and under what environmental conditions (Wray et al., 2003). During the 158 transcriptional regulation process, two types of regulatory units exist: one is cis-regulatory 159 elements or sequences that flank every gene; the other is trans-acting factors that encode 160 elsewhere and physically interplay with cis-elements to either promote or repress the transcription 161 of the gene (Wray et al., 2003). Correspondingly, gene expression modifications can arise from 162 either cis-regulatory changes that affect transcription in an allele-specific manner, or 163 trans-regulatory alterations that modify the expression or activity of related factors (Wray et al., 164 2003;Wittkopp et al., 2004). Recently, a burgeoning amount of data have shown that distinct 165 portions of cis-and trans-regulatory changes contribute to gene expression divergence in different 166 organisms by quantifying allele-specific expression (ASE) levels in a common, heterozygous 167 diploid individual (Wittkopp et al., 2004(Wittkopp et al., , 2008Kiekens et al., 2006;Stupar et al., 2007;Zhuang 168 and Adams, 2007;Main et al., 2009;von Korff et al., 2009;Bell et al., 2013;Cubillos et al., 2014). 169 In spite of advances regarding cis-and trans-actions, the detailed relationship between regulatory 170 variations of key developmental genes and phenotypic divergences, especially some important 171 morphological novelties, has not been extensively explored yet, particularly in plants (Rosas et al., 172 2010). 173 Lamiales is a large clade in angiosperms believed to be ancestrally zygomorphic and further 174 specialized and diversified (Cubas, 2004). As one of the earliest diverging groups in Lamiales, 175 Gesneriaceae has diverse forms of zygomorphic flowers (Cubas, 2004). Recent studies have 176 shown that the spatial-temporal expression alterations of CYC-like genes are highly correlated 177 with the diverse morphologies of zygomorphic flowers in Gesneriaceae (Du and Wang, 2008;Gao 178 et al., 2008;Song et al., 2009;Yang et al., 2012). The genus Petrocosmea, comprising 42 species, 179 is a mid-sized group in Gesneriaceae with very similar vegetative traits but extremely diverse 180 www.plantphysiol.org on August 20, 2017 -Published by Downloaded from Copyright © 2015 American Society of Plant Biologists. All rights reserved. zygomorphic flowers that are mainly reflected in size and shape variations of the dorsal petals in 181 the second floral whorl (Li and Wang, 2004;Qiu et al., 2011Qiu et al., , 2015. Petrocosmea therefore 182 represents an ideal group to address whether different regulatory pathways or their interplay 183 account for the expression differentiation of floral symmetry related genes and therefore the 184 morphological diversification of the dorsal petals in zygomorphic flowers. 185 In this study, two Petrocosmea species, i.e. P. glabristoma and P. sinensis as the 186 representatives of two clades, were selected to investigate possible regulatory pathways driving 187 the morphological divergence of zygomorphic flowers by employing multiple approaches, 188 including comparative expression analyses of CYC-and CIN-like TCP genes, hybridization, 189 mutant, and ASE assays. The results showed that the expression differentiation of CYC1C,190 CYC1D, and CIN1 coordinately promoted the morphological divergence of the dorsal petals that 191 is characteristic of different types of zygomorphic flowers in the two species. In addition, both 192 CYC-and CIN-like genes conformed to an additive inheritance mode in reciprocal F1 hybrids, 193 consistent with the phenotypes of the hybrids. We further found that cis-regulatory changes,

RESULTS 201
Floral Morphology of Parents and F1 Hybrids 202 P. glabristoma and P. sinensis belong to clade D and clade E in the genus Petrocosmea, 203 respectively (Qiu et al., 2015;Supplemental Fig. S1). They have very similar vegetative traits but 204 exhibit distinctly different zygomorphic flowers that are mainly manifested in the second whorl of 205 floral organs (Fig. 1, A and B). First, the dorsal petals of P. glabristoma flowers are significantly 206 short compared to those of P. sinensis. The dorsal petal length of P. glabristoma flowers is 4.83 ± 207 0.21 mm on average, and the length ratio of dorsal to ventral petals is 0.38 ± 0.01. In P. sinensis 208 flowers, the dorsal petal length is 8.18 ± 0.33 mm, and the length ratio of dorsal to ventral petals is 209 0.84 ± 0.03 (Table I). Second, the curvature degree of their dorsal petals is also distinctly different 210 from each other. As shown in Fig. 1, the dorsal petals of P. glabristoma are extended upward and 211 almost perpendicular to the corolla tube, while those of P. sinensis are strongly reflexed backward 212 ( Fig. 1, A and B, the side view). In addition to dorsal petals, their lateral and ventral petals are also 213 clearly different from each other in shape. For example, both lateral and ventral petals of P. 214 sinensis are strongly curved with lateral margins reflexed backward, while those of P. glabristoma 215 are explanate (Fig. 1, A' and B', the dorsal view). 216 Reciprocal F1 hybrids were constructed using P. glabristoma (Pg) as the maternal parent and 217 P. sinensis (Ps) as the paternal parent, or vice versa, to see how the petal phenotypes in the two 218 parents are inherited and to rule out parent-of-origin effects. Cross Pg × Ps and reciprocal cross Ps 219 × Pg have very similar floral phenotypes (Fig. 1,C and D). Their flowers exhibit an intermediate 220 phenotype between P. glabristoma and P. sinensis that is mainly manifested in the length of the 221 dorsal petals (Fig. 1). First, the dorsal petal length of Pg × Ps and Ps × Pg hybrids is 6.22 ± 0.47 222 mm and 5.74 ± 0.35 mm, respectively, both falling somewhere between the two parents (Table I). 223 Second, the length ratio of dorsal to ventral petals is respectively 0.58 ± 0.04 and 0.56 ± 0.03 for 224 Pg × Ps and Ps × Pg, also intermediate between P. glabristoma and P. sinensis (Table I). The 225 curvature degree of dorsal, lateral, and ventral petals in both Pg × Ps and Ps × Pg hybrids is more 226 similar to P. glabristoma, but still located between the two parents ( Fig. 1)

Sequence and Expression Analyses of CYC-and CIN-like Genes in Parents and F1 Hybrids 231
To know whether the morphological differences of the dorsal petals between P. glabristoma 232 and P. sinensis are related to CYC-like gene activity, we isolated four CYC-like genes, i.e. CYC1C, 233 CYC1D, CYC2A, and CYC2B from both species. Three CIN-like genes, CIN1, CIN2A,and CIN2B,234 were also isolated due to their potential roles in locally controlling cell growth and proliferation of 235 petals. Phylogenetic analyses showed that CYC-and CIN-like genes isolated here were clustered 236 with A. majus CYC/DICH and CIN, respectively (Supplemental Fig. S2). Sequence analyses 237 showed that Petrocosmea CYC-like genes are characterized by the presence of both TCP and R 238 domains, while CIN-like genes only contain the TCP domain (Supplemental Fig. S3). In addition, 239 all CIN-like genes isolated herein contain the miR319 recognition sequence (Supplemental Fig.  240 S3), a typical feature shared by CIN-like genes (Martín-Trillo and Cubas, 2010). Furthermore, 241 there is a very high amino acid sequence similarity between each pair of orthologous genes 242 (Supplemental Fig. S4 and S5), hinting that these CYC-and CIN-like genes are functionally 243 conserved between the two species. 244 Real-time PCR was conducted to survey the expression patterns of all CYC-like genes in the 245 second whorl of floral organs in P. glabristoma, P. sinensis, as well as in F1 hybrids. As shown in 246 and P. sinensis with very low expression signal in the lateral and ventral petals ( Fig. 2A). By 248 contrast, transcription signals of CYC2A and CYC2B were almost undetectable in all petals of P. 249 glabristoma and P. sinensis flowers ( Fig. 2A), consistent with a previous report (Gao et al., 2008). 250 Therefore, we only conducted further expression analyses of CYC1C and CYC1D in the dorsal 251 petals. For CYC1C, PgCYC1C was strongly expressed in the dorsal petals of P. glabristoma, while 252 PsCYC1C had a much lower expression level in the corresponding regions of P. sinensis (Fig. 3A). 253 The expression ratio of PgCYC1C and PsCYC1C was 1.97 ± 0.20 (Table II), consistent with the 254 morphological variations of the dorsal petals between P. glabristoma and P. sinensis. In Pg × Ps 255 hybrid, the expression signal of CYC1C was obviously weaker than in P. glabristoma, but stronger 256 than in P. sinensis. Similarly, its expression level in Ps × Pg hybrid also fell between the two 257 parents (Fig. 3A). Similar to CYC1C, CYC1D also showed a significant expression difference with 258 a ratio of 2.11 ± 0.27 between PgCYC1D and PsCYC1D ( Fig. 3A; Table II Real-time PCR was also performed to investigate the expression patterns of all CIN-like 267 genes in all floral organs of the second whorl in both parents and F1 hybrids. The results showed 268 that CIN1 expression levels were similar in the dorsal, lateral, and ventral petals in both P. 269 glabristoma and P. sinensis, but were much different between the two species (Fig. 2B). Different 270 from CIN1, CIN2A and CIN2B had similar expression levels between P. glabristoma and P. 271 sinensis in all petals (Fig. 2B), ruling out the possibility of their involvement in the petal shape 272 variation between the two species. Therefore, we further carried out expression analyses of only 273 CIN1 in both parents and hybrids. The results showed that PgCIN1 and PsCIN1 exhibited an 274 extreme expression divergence in the dorsal petals between P. glabristoma and P. sinensis with a 275 ratio of 0.09 ± 0.01 ( Fig. 3B; Table II), in line with the distinct curvature degree of the dorsal 276 petals that is upward extended in P. glabristoma while extremely backward reflexed in P. sinenis. 277 In addition, CIN1 also showed an extreme expression difference in both lateral and ventral petals 278 between P. glabristoma and P. sinensis (Fig. 3B), consistent with the different curvature status of 279 these petals between the two species. In Pg × Ps hybrid, CIN1 expression levels were higher than 280

Allele-specific Expression Analyses 300
Given that CYC1C, CYC1D, and CIN1 were differentially expressed in petals of P. 301 glabristoma and P. sinensis, we further conducted allele-specific expression analyses in reciprocal 302 www.plantphysiol.org on August 20, 2017 -Published by Downloaded from Copyright © 2015 American Society of Plant Biologists. All rights reserved.

F1 hybrids to determine what regulatory changes account for their expression differentiation 303
between the two parents. First, we identified the SNP sites for each orthologous pair by 304 sequencing the PCR products of CYC1C, CYC1D, and CIN1 amplified from about 30 different 305 genomic DNA samples of P. glabristoma and P. sinensis, respectively, which were further 306 confirmed by genotyping the reciprocal F1 hybrids (Supplemental Fig. S6). For each gene, all 307 sequences from 30 individuals of each species were identical, indicating the homozygosity of the 308 two parents. 309 In ASE assays, allele expression ratios in the reciprocal F1 hybrids for each gene were 310 normalized by genomic DNA in which allelic ratios were assumed to be present in equal amounts. 311 For CYC1C, the expression levels of PgCYC1C and PsCYC1C alleles in both Pg × Ps and Ps × Pg 312 hybrids were significantly different from each other with the ratios of 2.13 ± 0.15 and 1.98 ± 0.14, 313 respectively ( Fig. 5A; Table II). Because CYC1C showed the same expression bias between 314 parents and hybrids, cis-regulatory changes were inferred to account for, at least predominantly, 315 the observed expression changes between the two parents. For CYC1D, PgCYC1D and PsCYC1D 316 alleles showed no significant expression difference in both Pg × Ps and Ps × Pg hybrids with the 317 ratios of 1.17 ± 0.06 and 1.07 ± 0.04, respectively ( Fig. 5B; Table II). Since CYC1D showed 318 biased expression in parents but equal expression level in hybrids, the expression difference of 319 PgCYC1D and PsCYC1D genes between the two parents could be mainly explained by 320 trans-regulatory changes. 321 PgCIN1 and PsCIN1 alleles showed an allele-specific expression imbalance in the dorsal 322 petals of the reciprocal F1 hybrids with expression ratios of 0.64 ± 0.01 and 0.59 ± 0.01 in Pg × Ps 323 and Ps × Pg hybrids, respectively ( Fig. 5C; Table II). They also showed the extreme allele-specific 324 expression imbalance in both lateral and ventral petals of the reciprocal F1 hybrids ( Fig. 5C; Table  325 II). Therefore, cis-regulatory changes should be involved in the expression difference of PgCIN1 326 and PsCIN1 genes in the dorsal, lateral, and ventral petals between the two parents. Nevertheless, 327 the allele expression differences in all petals of the reciprocal F1 hybrids were much lower than 328 the expression bias in the corresponding regions of the two parents ( Fig. 5C), indicating that the 329 expression differentiation of CIN1 between the two parents might also be attributed to 330 trans-regulatory changes. The above results also indicate that there is little parent-of-origin effect 331 on the expression of CYC1C, CYC1D, as well as CIN1 in the reciprocal F1 hybrids. abundance of allele-specific transcripts in F1 hybrids (Pg F1 /Ps F1 ; Fig. 5) according to Wittkopp et 336 al. (2004Wittkopp et 336 al. ( , 2008. Expression ratio for each gene was log 2 transformed at first. Then, the expression 337 ratio between parents (log 2 Pg/Ps) was plotted against the allelic expression ratio in F1 hybrids 338 (log 2 Pg F1 /Ps F1 ). The relative contribution of cis-and trans-regulatory changes was determined 339 based on the position of the point, as a plot of the relative allelic expression bias against the 340 relative parental expression differentiation (Fig. 5). For CYC1C, the allelic expression ratio in the 341 Pg × Ps hybrid and the expression ratio between parents were almost identical (log 2 Pg/Ps = 342 log 2 Pg F1 /Ps F1 ) with the point fallen on the diagonal (y = x) (Fig. 5D). Similarly, the point for the 343 allelic expression ratio in the Ps × Pg hybrid and the parental expression ratio fell on the y = x line 344

DISCUSSION 358
Petrocosmea is divided into five clades according to molecular phylogenetic analyses (Qiu et 359 al., 2015;Supplemental Fig. S1). In clade A, the two dorsal petals of the upper lip are slightly 360 shorter than the three petals of the lower lip. During the successive morphological evolution and 361 diversification process, the two dorsal petals tend to be shortened and specialized in size and shape 362 from clade A through B and C to clade D, with coordinated variations of correlative characters. upstream enhancer plays an important role in the domestication process of modern maize from its 381 wild ancestor teosinte (Doebley et al., 1997;Clark et al., 2006;Studer et al., 2011). 382 Floral zygomorphy as a key novelty contributes significantly to the explosive radiation of 383 angiosperms (Dilcher, 2000). The modified expression of CYC-like TCP genes usually plays a 384 decisive role in generating novel morphologies in floral symmetry (Hileman and Cubas, 2009 or ventralized peloric flowers (Luo et al., 1996(Luo et al., , 1999. Similar phenomena have been widely 388 observed in other angiosperms (Cubas et al., 1999b;Citerne et al., 2006;Busch and Zachgo, 2007;389 Pang et al., 2010;Yang et al., 2012). Furthermore, morphological transitions in zygomorphic 390 flowers and between zygomorphy and actinomorphy in closely related species are usually related 391 to the expression changes of CYC-like genes (Zhang et al., 2010;Howarth et al., 2011;Busch et 392 al., 2012;Zhong and Kellogg, 2015). As representatives of two different clades in Petrocosmea, P. in other Gesneriaceae plants (Yang et al., 2012;Liu et al., 2014). Therefore, we conclude that the 404 expression differentiation of CYC1C and CYC1D, specifically controlling the dorsal petals, brings 405 about, at least partially, the morphological divergence of the dorsal petals between P. glabristoma 406 and P. sinensis. 407 Different from CYC-like genes, CIN-like genes mainly control the local morphology of 408 leaves and petals with curvature produced upon their up-or down-regulation (Nath et al., 2003;409 Palatnik et al., 2003;Crawford et al., 2004;Koyama et al., 2007Koyama et al., , 2010aKoyama et al., , 2011Nag et al., 2009;410 Gupta et al., 2014). In this study, P. glabristoma has upward extended dorsal petals and explanate 411 lateral/ventral petals, while P. sinensis has highly backward reflexed dorsal petals and obviously Hybridization usually results in phenotypic variations that at least in part attribute to 430 modifications of gene expression Chen, 2010). A large number of 431 studies have shown that the majority of genes differentially expressed between two parents are 432 often expressed at additive levels in the F1 hybrid offspring (Guo et al., 2004(Guo et al., , 2006Stupar and 433 Springer, 2006;Swanson-Wagner et al., 2006;Stupar et al., 2007). In this study, expression levels

TCP Genes 449
Differences in gene expression can be caused by either cis-regulatory changes that result 450 from changes in enhancer or promoter sequences, changes in the transcribed region, or changes in 451 chromatin structure, or trans-regulatory changes that are caused by genetic and epigenetic changes 452 affecting the activity or availability of proteins and RNAs that mediate gene expression (Wittkopp 453 et al., 2004(Wittkopp 453 et al., , 2008. While cis-regulatory changes alter gene expression in an allele-specific manner, 454 trans-regulatory variations influence expression of both alleles in a diploid cell. It has been 455 hypothesized that natural selection may favor one type of change over the other (Wittkopp et al., 456 2008). 457 In this study, we examined allele-specific expression patterns of CYC-and CIN-like genes in 458 F1 hybrids to elucidate the regulatory mechanisms underlying their expression modifications 459 between P. glabristoma and P. sinensis. ASE analyses preclude the possibility of parent-of-origin 460 (or genome imprinting) effect accounting for allele expression differences because of identical 461 allele-specific expression patterns of both CYC-and CIN-like genes in the reciprocal F1 hybrids. 462 The results further demonstrate that the expression differentiation of these TCP genes is underlain 463 by distinctly different types of regulatory changes. The expression differentiation of CYC1C is 464 only (or at least predominantly) caused by cis-regulatory changes, while that of CYC1D is 465 predominantly underlain by trans-acting variations. In angiosperms, especially in eudicots, the 466 CYC2 clade genes have undergone frequent duplications that gave rise to multiple copies of 467 CYC2 genes in most members of zygomorphic clades (Howarth and Donoghue, 2006). Of these, a 468 pair of CYC2 genes are usually involved in controlling the dorsal petals redundantly in 469 zygomorphic flowers with others having no or transient expression signals (Luo et al., 1996(Luo et al., , 1999470 Citerne et al., 2006;Feng et al., 2006;Gao et al., 2008;Wang et al., 2008;Song et al., 2009;471 Zhang et al., 2010;Howarth et al., 2011;Yang et al., 2012). A recent study shows that a double 472 positive autoregulatory feedback loop has evolved repeatedly in different zygomorphic clades of 473 angiosperms to maintain the expression of these gene pairs (Yang et al., 2012) mechanisms would be wholly different from each other as mentioned here (Fig. 6). 478 Natural selection driving morphological diversification has received much attention at the 479 level of protein evolution, but very little for gene expression evolution (Fraser, 2011  In addition, genome duplication initially generates redundant gene copies followed by 497 evolution and functional divergence. Many duplicate genes can retain a considerable degree of 498 expressional or functional similarity over very long periods, which is usually assumed due to 499 purifying selection maintaining the ancestral function in both duplicates (Lynch and Conery, 2000;500 Kondrashov and Kondrashov, 2006;Dean et al., 2008;Vavouri et al., 2008). This may be true with 501 respect to the preservation of the CYC-like TCP genes in a genome for a long period of time in 502 zygomorphic clades. However, the high redundancy of CYC1C and CYC1D with identical 503 expression patterns and interspecific expression differentiation in Petrocosmea species is merely 504 shaped by the regulatory context rather than their own coding sequences. This similarity is 505 maintained by remarkably different regulatory mechanisms (Fig. 6). Therefore, our findings herein 506 pose a challenge to the traditional way of evaluating and understanding the role of selection in 507 fixation of gene duplications and their evolutionary fates. 508 Unlike CYC-like genes, CIN1 expression divergence between P. glabristoma and P. sinensis portion of expression changes. It has been reported that miR319 is involved in preventing aberrant 511 activity of CIN-like genes expressed from its native promoter in both leaves and petals of A. 512 thaliana by direct mRNA cleavage (Palatnik et al., 2003;Nag et al., 2009). In TCP gene family, 513 only CIN lineage genes contain sequences matching the miR319 target sequence (Palatnik et al., 514 2003;Martín-Trillo and Cubas, 2010). Therefore, trans-regulatory changes of CIN1 might be 515 related to the modification of miR319 activity and function. This scenario might be sound because 516 both PgCIN1 and PsCIN1 contain potential target sequences of miR319. It is also possible that the 517 interaction of modified miR319 activity and changed cis-elements might have finally given rise to 518 the expression divergence of CIN1 between P. glabristoma and P. sinensis (Fig. 6). It would be 519 important to further decipher how the miR319 activity is co-opted to the CIN1 regulatory pathway 520 and whether there is a complex combinatorial mechanism underlying a possible dynamic 521 regulatory interaction or network in and between CYC-and CIN-like genes, and if, how they are 522 modified coordinately in creating the expression differentiation of these genes responsible for 523 related morphological evolution in zygomorphic flowers. 524

CONCLUSIONS 526
We here report that CYC-and CIN-like TCP genes have undergone expression differentiation, 527 strongly correlated with the morphological divergence between P. glabristoma and P. sinensis. It is 528 the first report that CYC-and CIN-like TCP genes act together in patterning a specific morphology 529 of floral zygomorphy with their coordinate expression differentiation giving rise to phenotypic 530 variations (Fig. 6). In F1 hybrids, their expression conforms to an additive pattern, consistent with 531 the petal phenotypes. We further reveal a new phenomenon that significantly distinct regulatory 532 pathways underlie the identical expression patterns and interspecific expression differentiation of 533 highly redundant paralogous genes, in which cis-regulatory changes contribute to the CYC1C 534 expression differentiation while trans-acting variations account for the CYC1D expression 535 divergence. The remarkablely differential expression of CIN1 between the two species is mediated 536 coordinately by cis-and trans-regulatory changes. Given that the role of selection in fixation of  The F1 hybrid seeds were surface-sterilized in 70% ethanol for 1 min, rinsed with sterile 556 water once, disinfected with 2.5% sodium hypochlorite for 3-5 min, and finally rinsed five times 557 with sterile water. The sterilized seeds were germinated on MS medium in a growth chamber at 26℃ 558 under a photoperiod of 10 h light /14 h dark (100 µmol m -2 s -1 ). The plantlets were transferred to 8 559 cm pots and grown under the same conditions as the parents. 560

DNA Extraction and Gene Isolation 562
Genomic DNA of P. glabristoma and P. sinensis was extracted from fresh leaves following a 563 modified CTAB procedure of Doyle and Doyle (1987). A forward primer F1 and a reverse primer 564 Gao et al., 2008;Pang et al., 2010) were used to amplify partial coding sequences of CYC-like 565 genes. To isolate CIN-like genes, a pair of primers PCIN-F and PCIN-R were designed according 566 to the sequences of A. majus CIN and Arabidopsis thaliana CIN-like genes. The full-length coding 567 sequences of CYC-and CIN-like genes were further isolated using Genome Walking Kit (TaKaRa). 568 All primers were listed in Supplemental Table S1. The PCR products were cloned into the and CIN-like genes. Neighbor-joining analyses were performed using the amino acid sequences of 575 Petrocosmea CYC-and CIN-like genes, CYC/DICH/CIN genes from A. majus (Luo et al., 1996(Luo et al., , 576 1999Nath et al., 2003), one CYC-like genes (AtTCP1) and eight CIN-like genes 577 (AtTCP2/3/4/5/10/13/ 17/24) from A. thaliana (Martín-Trillo and Cubas, 2010), and four CYC-like 578 genes (PhCYC1C/1D/2A/2B) from Primulina heterotricha (Yang et al., 2012). The Oryza sativa 579 PCF1 and PCF2 (Kosugi and Ohashi, 1997) were used as outgroups. The sequences were first 580 aligned using Clustal X (Thompson et al., 1997) and further adjusted manually with BioEdit (Hall, 581 1999). MEGA6 software (SudhirKumar) was used to construct Neighbor-joining trees using the 582 following settings: p-distance, pairwise deletion, with the bootstrap values calculated for 1000 583 replicates. 584

RNA Extraction and Real-time PCR 586
Dorsal, lateral, and ventral petals of P. glabristoma, P. sinensis, Pg × Ps, and Ps × Pg flowers 587 before anthesis were dissected, immediately frozen in liquid nitrogen, and stored at -80℃ until 588 extraction. Total RNA was extracted using an SV Total RNA Isolation System (Promega) 589 according to the manufacturer's instructions. cDNA was synthesized using a RevertAid H Minus 590 First-Strand cDNA Synthesis Kit (Thermo). Gene-specific primers were used to amplify all genes, 591 and ACTIN was amplified as a reference gene (see Supplemental Table S1 online). It is 592 noteworthy that in real-time PCR, different alleles for each TCP gene in F1 hybrids were 593 indiscriminately amplified using a pair of primers that anneal to the conserved regions. The 594 specificity of all primers was confirmed by sequencing PCR products. The SYBR Premix Ex Taq 595 (TaKaRa) was used to perform real-time PCR with ROX as a reference dye using the StepOne 596 Plus Real-Time PCR System (AB Applied Biosystems). The PCR conditions were: initial 597 denaturation at 95℃ 30 sec, and 40 cycles of 95℃ 10 sec and 60℃ 40 sec, after which 598 dissociation curves were recorded using 1 cycle of 95℃ 30 sec, 60℃ 30 sec, and 95℃ 30 sec. The 599 relative expression levels were determined by normalizing the PCR threshold cycle number of 600 each gene with that of ACTIN. The data were the average of three biological replicates with each 601 including three technical replicates. The expression ratio of each gene between paternal and 602 maternal parents (Pg/Ps) was calculated and compared using the Fisher's Least Significant 603 Difference (LSD) test to identify the parental expression imbalance (H 0 : Pg/Ps = 1, P < 0.01). All 604 statistical analyses were performed using SPSS software v. 14.0 (SPSS Inc.). 605

SNP Identification and ASE Assays by Single-base Primer Extension Experiment 607
For SNP identification, at least 30 individuals from either P. glabristoma or P. sinensis were 608 sampled to isolate genomic DNA. DNA fragments for CYC1C, CYC1D, and CIN1 were amplified 609 using real-time PCR primers and the PCR products were directly sequenced. The homozygosity of 610 two parents was identified as the identity of all sequences from different individuals of a species 611 for each gene. SNP sites were identified as single-base differences between two parents. The SNP 612 sites were further confirmed by genotyping reciprocal F1 hybrids in which two peaks would 613 appear in a position corresponding to two different reads respectively from two parents in the 614 same position using the Chromas software (Chromos MFC Application). One fixed SNP site for 615 each gene was selected for ASE analyses. 616 For ASE assays, the gene-specific segment surrounding the SNP site was first amplified from 617 the cDNA pool of reciprocal F1 hybrids using real-time PCR primers. Single-base extension 618 experiments were conducted using the ABI P RISM ® SNaPshot TM Multiplex Kit (AB Applied 619 Biosystems) according to the manufacturer's instructions. Briefly, following PCR thermal cycling, 620 unincorporated primers and dNTPs were removed by adding 5 units of Shrimp Alkaline 621 Phosphatase (TaKaRa) and 1 unit of Exonuclease I (TaKaRa) to each 15 µl PCR product. 622 Reactions were mixed briefly and incubated at 37℃ for 60 min, followed by 75℃ for 15 min to 623 inactivate the enzyme. The PCR products were then subjected to a primer extension assay using 624 single-base extension primers that were designed to anneal to the amplified DNA fragments 625 adjacent to the SNP sites (Supplemental Table S1). Primer extension reactions were carried out in 626 a total volume of 10 µl containing 5 µl of SNaPshot Multiplex Ready Reaction Mix, 0.4 µL of 5 627 µM extension primer, 2 µl of PCR product, and 2.6 µl of double distilled water. Thermal cycling 628 conditions for extension reactions were carried out with the following program: initial 629 denaturation at 94℃ 2 min, and 25 cycles of 96℃ 10 sec, 50℃ 5 sec, and 60℃ 30 sec. After cycling, the unincorporated fluorescent ddNTPs were removed by adding 1 unit of SAP and 631 incubating for 60 min at 37℃, followed by 15 min at 65℃ for enzyme inactivation. The resulting 632 primer extension products were analyzed on an ABI 3730 capillary electrophoresis DNA 633 instrument using the Peak Scanner Software v 1.0 (AB Applied Biosystems) according to the 634 manufacturer's protocol. The expression ratio of two alleles was measured by comparing the 635 height of peaks. Allelic ratios of genomic DNA from F1 hybrids assumed to be present in equal 636 amounts (1:1) were used to normalize allelic ratios of cDNA samples to adjust PCR bias 637 (Wittkopp et al., 2004(Wittkopp et al., , 2008Zhuang and Adams, 2007). Standard errors were calculated from 638 three independent biological replicates (including three technical replicates for each). For each 639 gene, the expression ratio of two alleles in F1 hybrids (Pg F1 /Ps F1 ) was compared with interspecific 640 expression bias of two orthologs in two parents (Pg/Ps) to identify cis-regulatory divergence by 641 using the Fisher's LSD test (H 0 : Pg F1 /Ps F1 = 1, P < 0.01). The statistical analyses were performed 642 using SPSS software v. 14.0. 643 644

Determination of Cis-and Trans-Regulatory Changes 645
Gene expression modifications can result from either cis-or trans-regulatory changes. To 646 discriminate the two changes, we compared the expression ratio of two alleles in heterozygous F1 647 hybrids with the expression ratio of two orthologs between two homozygous parents. In the same 648 genetic background of hybrids, two alleles inherited from different parents are exposed to a 649 common set of trans-acting factors. Therefore, any observed allelic expression imbalance in 650 hybrids is characteristic of cis-regulatory divergence. Trans-acting variation can then be deduced 651 if genes show expression differences in two parents but no allelic expression imbalance in hybrids. 652 To further determine relative contributions of cis-and trans-regulatory changes to gene 653 expression differences between two parents, we compared the expression ratio of two orthologs 654 between parents (Pg/Ps) with the expression ratio of two alleles in F1 hybrids (Pg F1 /Ps F1 ) 655 according to Wittkopp et al. (2004Wittkopp et al. ( , 2008. Expression ratios for each gene were log 2 transformed 656 at first. Then, the expression ratio between parents (log 2 Pg/Ps, the horizontal axis) was plotted 657 against the allelic expression bias in F1 hybrids (log 2 Pg F1 /Ps F1 , the vertical axis). If the hybrid and 658 parental expression ratios are identical (log 2 Pg/Ps = log 2 Pg F1 /Ps F1 ) and points fall on the diagonal 659 (y = x), cis-regulatory changes are believed to completely explain the expression differences 660 The following supplemental materials are available. 673 Supplemental Figure S1. The majority rule consensus Bayesian tree generated from analysis 674 of combined cpDNA and nDNA data.