Abstract Dissecting complex connections between cytogenetic traits (ploidy levels) and plant invasiveness has emerged as a popular research subject in the field of invasion biology. Although recent work suggests that polyploids are more likely to be invasive than their corresponding diploids, the molecular basis underlying the successful invasion of polyploids remains largely unexplored. To this end, we adopted an RNA‐seq and sRNA‐seq approach to describe how polyploids mediate invasiveness differences in two contrasting cytotypes of Solidago canadensis L., a widespread wild hexaploid invader with localized cultivated diploid populations. Our analysis of the leaf transcriptome revealed 116,801 unigenes, of which 12,897 unigenes displayed significant differences in expression levels. A substantial number of these differentially expressed unigenes (DEUs) were significantly associated with the biosynthesis of secondary metabolites, carbohydrate metabolism, lipid metabolism, and environmental adaptation pathways. Gene Ontology term enrichment‐based categorization of DEU‐functions was consistent with this observation, as terms related to single‐organism, cellular, and metabolic processes including catalytic, binding, transporter, and enzyme regulator activity were over‐represented. Concomitantly, 186 miRNAs belonging to 44 miRNA families were identified in the same leaf tissues, with 59 miRNAs being differentially expressed. Furthermore, we discovered 83 miRNA‐target interacting pairs that were oppositely regulated, and a meticulous study of these targets depicted that several unigenes encoding transcription factors, DNA methyltransferase, and leucine‐rich repeat receptor‐like kinases involved in the stress response were greatly influenced. Collectively, these transcriptional and epigenetic data provide new insights into miRNA‐mediated gene expression regulatory mechanisms that may operate in hexaploid cytotypes to favor successful invasion. Keywords: invasiveness, miRNA, ploidy, polyploid, Solidago canadensis L., transcriptome 1. INTRODUCTION Amid growing evidence of biological invasion impacts on global biodiversity, ecosystem functioning, species conservation, and even social and economic activities (Dyer et al., [28]2017; Rejmánek, [29]2015), there is mounting interest in searching the determining elements underlying the successful invasion of alien species. A key element associated with successful invasion of alien species is their capacity for rapid adaptation to environmental challenges following introduction (Huang et al., [30]2017). Identifying critical traits that benefit this rapid environmental adaptation is therefore a talking point for conservation concern as it can create greater opportunities to predict the invasion risk related to diverse alien species. A steady stream of ecological research in recent years has identified a variety of shared ecological traits related to invasiveness among invasive alien species, such as increased growth and fecundity, wide ecological tolerance, high fitness, and strong clonal propagation (Richardson & Pyšek, [31]2008). As ecological traits have not generated an explicit recognition pattern in plant invasiveness, some invasion biologists have switched from ecological traits to genetic or genomic traits in comprehending patterns of plant invasiveness (Hodgins, Lai, Nurkowski, Huang, & Rieseberg, [32]2013; Prentis & Pavasovic, [33]2013; Rius & Darling, [34]2014). For instance, Rius and Darling ([35]2014) showed that genetic admixture may act as a genuine driver to play central roles in the successful invasion of genetically admixed individuals. Likewise, another related study performed on invasive alien plants showed a statistical connection between ploidy level and invasiveness, concluding that polyploid plants were more likely to be invasive than diploids and that ploidy level (and chromosome number) was positively related to plant invasiveness (Pandit, Pocock, & Kunin, [36]2011; Pandit, White, & Pocock, [37]2014). Polyploids, recognized as organisms that possess more than two complete sets of chromosomes in their somatic cells, have often been suggested to represent a powerful driver in the speciation, evolution, and adaptation of plants, with far‐reaching ecological and evolutionary consequences. On the other hand, polyploids often appear to be over‐represented in invasive plants (Thébault, Gillet, Müller‐Schärer, & Buttler, [38]2011) and have cumulatively been acknowledged as a latent advantageous attribute of plant invaders (Pandit et al., [39]2011; te Beest et al., [40]2012). Furthermore, the influences of polyploids can act as a cascade process that directly or indirectly mediates virtually all aspects of plant genetics, morphology, physiology, life history, and ecology (Levin, [41]1983). Therefore, many evolutionary biologists believe that polyploids provide introduced plants with new features that permit them to invade largely varied environments or expand their geographical range. However, this hypothesis has not yet been proved efficiently. Ecological studies have conventionally focused on comparing clearly relevant diploid and polyploid in their native and introduced ranges (Hahn, Buckley, & Müller‐Schärer, [42]2012), but few studies have elucidated the molecular basis underlying the invasiveness difference in alien plants with different ploidy levels (cytotypes). Therefore, the genetic and epigenetic impact imposed by ploidy alteration remains elusive. Studies on natural and synthetic polyploids have repeatedly revealed that rapid and dynamic changes at the genetic, gene expression, and epigenetic levels occur after polyploid formation (Chen, [43]2007; Jackson & Chen, [44]2010; Sun, Wu, et al., [45]2017). Likewise, evidence is also mounting that epigenetic modifications can change gene expression and reconstruct gene expression networks thus resulting in pronounced phenotypic alterations (Hao, Lucero, Sanderson, Zacharias, & Holbrook, [46]2013; Song & Chen, [47]2015) and allowing polyploids to occupy new habitats, grow vigorously and improve adaptation in novel environments (Madlung, [48]2013). As an extensive type of epigenetic modifications in nonmodel organisms, miRNA has attracted considerable concern due to its regulatory mechanisms for gene expression. Additionally, miRNA is highly conserved in evolution but becomes activated in polyploidization (Axtell, [49]2008; Ha et al., [50]2009). More importantly, alterations in miRNA expression can mediate their target‐gene expression at the post‐transcriptional level, and this effect is viewed as one of the main reasons for phenotypic changes of polyploids (Chen, [51]2007; Ha et al., [52]2009). Accordingly, elucidating the divergences in gene and miRNA expression between different ploidy levels (cytotypes) of alien plants and how they influence phenotypic differentiation is crucial to explain how polyploids might have contributed to successful invasion. Herein, we evaluate the effect of polyploids on gene and miRNA expression while also considering the potential roles of miRNA‐mediated gene expression regulation in driving differences in invasiveness between diploid and hexaploid cytotypes of Solidago canadensis L. Specifically, the objectives of our current work are as follows: (a) to characterize the initial expression profiling of genes and miRNAs in two identified ploidy levels, that is, diploid and hexaploid cytotypes of S. canadensis at a genome‐wide scale; (b) to determine key candidate genes and miRNA regulators that may contribute to the successful invasion of hexaploid cytotypes based on gene and miRNA expression divergences, as well as over‐represented functional categories of these candidates; and (c) finally to explore the strong evidence for the potential genetic roles of epigenetic and transcriptional alterations in the successful invasion of hexaploid cytotypes. To this end, we adopted an RNA‐seq and sRNA‐seq approach to investigate the divergences of gene and miRNA expression between diploid and hexaploid cytotypes of S. canadensis. Furthermore, we constructed a co‐expression network of differentially expressed genes and miRNAs to shed light on the regulatory action of miRNAs. Taken together, our work provides new insights into miRNA‐mediated gene expression regulatory mechanisms that may be useful to explain the successful invasion of hexaploid cytotypes. 2. MATERIALS AND METHODS 2.1. Study species Solidago canadensis (Asteraceae), a perennial weed native to North America (Werner, Bradbury, & Gross, [53]1980) where it exists in a diploid, tetraploid, or hexaploid cytotype (Melville & Morton, [54]1982), has invaded a wide geographical range globally, including New Zealand, Australia, Europe, and Asia (Abhilasha, Quintana, Vivanco, & Joshi, [55]2008; Szymura, Szymura, Wolski, & Swierszcz, [56]2018). Introduced into eastern China in the 1930s as an ornamental plant, S. canadensis began to escape cultivation and spread in the 1980s. Currently, it has become highly abundant and has noticeably affected the diversity and richness of native plant species (Wang, Jiang, Zhou, & Wu, [57]2018). However, it is worthwhile noting that, only hexaploid cytotypes of S. canadensis have long been convincingly reported to occur widely in the introduced range in China and become invasive thus far (Wang, [58]2016). Their corresponding diploid cytotypes (also called “Huang Ying” in China) were cultivated mainly in Yunnan Province in southwestern China as an important cut‐flower plant. An earlier experiment carried out with common gardens showed that the growth of hexaploid cytotypes of S. canadensis was more vigorous than their related diploids, offering clear advantages for the successful invasion of hexaploid cytotypes (Li, [59]2011). Additionally, the roots, stems, and leaves of hexaploid cytotypes were morphologically and anatomically distinct from their diploids (Wang, [60]2007). Overall, the contrasting invasive propensities and geographical and phenotypic differentiation between hexaploid cytotypes and their related diploids make S. canadensis an excellent study system to answer such questions as how polyploids both affect gene and miRNA expression and alter molecular pathways that may be responsible for the successful invasion of hexaploid cytotypes and whether miRNA plays key roles in reprogramming the transcriptional expression. 2.2. Population sampling and chromosome counting Invasive populations of S. canadensis largely cluster around the Yangtze River Delta, which occupies its main distribution range in China (Figure [61]1, Figure [62]A1 in Appendix), and cultivated populations were narrowly cultivated in Yunnan Province in southwestern China. Therefore, we sampled invasive populations (separated from each other by at least 3 km) in 42 locations throughout Jiangsu, Zhejiang, Anhui, Hubei, and Shanghai and one cultivated population from Yunnan Province (Table [63]1; Figure [64]1). Finally, 449 sampled individuals representing 43 populations were collected and subsequently transplanted into pots with commercial soil and grown for 3 months in the greenhouse of Wuhan University under natural photoperiod conditions. Chromosome counting was performed according to the modified carbol fuchsin squash method (The detailed methods are presented in Supporting Information Appendix [65]S1). Figure 1. Figure 1 [66]Open in a new tab Map of China sampling sites for 43 populations of S. canadensis described in Table [67]1. The circles indicate sampling locations. The blue and yellow circles show the diploid and hexaploid populations used in gene and miRNA expression analyses, respectively Table 1. Geographical coordinates and sample size of 43 populations of S. canadensis in China No. Population code Location Geographical coordinates Status No.of samples Chromosome number Latitude (N) Longitude (E) 1 MH 1 Minhang District, Shanghai City N31°08′41.40″ E121°23′19.51″ Invasive 10 54 2 MH 2 Minhang District, Shanghai City N31°09′54.98″ E121°20′53.35″ Invasive 10 54 3 MH 3 Minhang District, Shanghai City N31°13′02.09″ E121°18′33.51″ Invasive 10 54 4 SJ 1 Songjiang District, Shanghai City N31°05′37.02″ E121°11′39.99″ Invasive 10 54 5 SJ 2 Songjiang District, Shanghai City N31°06′01.76″ E121°12′17.63″ Invasive 10 54 6 PD Pudong District, Shanghai City N31°15′14.57″ E121°38′25.04″ Invasive 10 54 7 GY Guanyun County, Lianyungang City, Jiangsu Province N34°23′34.49″ E119°14′15.33″ Invasive 8 54 8 XP Xinpu District, Lianyungang City, Jiangsu Province N34°38′24.10″ E119°11′06.86″ Invasive 8 54 9 YD 1 Yandu District, Yancheng City, Jiangsu Province N33°18′14.48″ E120°06′26.04″ Invasive 8 54 10 YD 2 Yandu District, Yancheng City, Jiangsu Province N33°18′15.30″ E120°06′45.48″ Invasive 8 54 11 DF Dafeng District, Yancheng City, Jiangsu Province N33°10′44.02″ E120°22′23.67″ Invasive 8 54 12 TC 1 Taicang City, Jiangsu Province N31°32′16.26″ E121°07′48.35″ Invasive 10 54 13 TC 2 Taicang City, Jiangsu Province N31°26′51.43″ E121°06′33.75″ Invasive 10 54 14 KS Kunshan City, Jiangsu Province N31°13′49.32″ E121°01′35.66″ Invasive 10 54 15 BH 1 Binhu District, Wuxi City, Jiangsu Province N31°32′01.06″ E120°09′22.39″ Invasive 10 54 16 BH 2 Binhu District, Wuxi City, Jiangsu Province N31°29′51.43″ E120°24′50.49″ Invasive 10 54 17 BH 3 Binhu District, Wuxi City, Jiangsu Province N31°29′27.02″ E120°27′01.11″ Invasive 10 54 18 WJ Wujin District, Changzhou City, Jiangsu Province N31°31′26.74″ E120°03′33.15″ Invasive 10 54 19 QX Qixia District, Nanjing City, Jiangsu Province N32°06′50.22″ E118°52′17.00″ Invasive 10 54 20 XW 1 Xuanwu District, Nanjing City, Jiangsu Province N32°05′46.33″ E118°53′07.09″ Invasive 10 54 21 XW 2 Xuanwu District, Nanjing City, Jiangsu Province N32°05′21.80″ E118°50′13.38″ Invasive 10 54 22 PK Pukou District, Nanjing City, Jiangsu Province N32°04′20.65″ E118°36′43.51″ Invasive 10 54 23 JG Jianggan District, Hangzhou City, Zhejiang Province N30°17′29.55″ E120°14′22.75″ Invasive 12 54 24 XS 1 Xiaoshan District, Hangzhou City, Zhejiang Province N30°11′33.37″ E120°16′23.98″ Invasive 10 54 25 XS 2 Xiaoshan District, Hangzhou City, Zhejiang Province N30°16′30.02″ E120°17′05.09″ Invasive 11 54 26 XS 3 Xiaoshan District, Hangzhou City, Zhejiang Province N30°07′19.58″ E120°15′41.75″ Invasive 9 54 27 BJ 1 Binjiang District, Hangzhou City, Zhejiang Province N30°10′18.65″ E120°08′13.93″ Invasive 9 54 28 BJ 2 Binjiang District, Hangzhou City, Zhejiang Province N30°09′26.83″ E120°08′07.79″ Invasive 12 54 29 XH Xihu District, Hangzhou City, Zhejiang Province N30°08′28.35″ E120°04′20.74″ Invasive 10 54 30 NH 1 Nanhu District, Jiaxing City, Zhejiang Province N30°45′04.50″ E120°44′21.24″ Invasive 13 54 31 NH 2 Nanhu District, Jiaxing City, Zhejiang Province N30°43′32.96″ E120°44′53.83″ Invasive 10 54 32 HN Haining City, Zhejiang Province N30°26′47.84″ E120°23′35.61″ Invasive 9 54 33 JH 1 Jinghu District, Wuhu City, Anhui Province N31°21′35.87″ E118°22′52.54″ Invasive 12 54 34 JH 2 Jinghu District, Wuhu City, Anhui Province N31°22′10.33″ E118°22′13.17″ Invasive 12 54 35 SS 1 Sanshan District, Wuhu City, Anhui Province N31°13′15.10″ E118°13′19.19″ Invasive 12 54 36 SS 2 Sanshan District, Wuhu City, Anhui Province N31°12′54.60″ E118°16′31.88″ Invasive 12 54 37 LM Lion Mountain District, Tongling City, Anhui Province N30°55′23.64″ E117°51′06.76″ Invasive 12 54 38 GC 1 Guichi District, Chizhou City, Anhui Province N30°36′19.83″ E117°30′10.12″ Invasive 12 54 39 GC 2 Guichi District, Chizhou City, Anhui Province N30°37′35.69″ E117°29′16.50″ Invasive 12 54 40 HS 1 Hongshan District, Wuhan City, Hubei Province N30°32′53.42″ E114°31′17.84″ Invasive 7 54 41 HS 2 Hongshan District, Wuhan City, Hubei Province N30°32′36.65″ E114°24′53.26″ Invasive 9 54 42 HS 3[68]a Hongshan District, Wuhan City, Hubei Province N30°32′22.40″ E114°25′01.12″ Invasive 14 54 43 CG[69]a Chenggong District, Kunming City, Yunnan Province N24°55′05.42″ E102°47′51.01″ Cultivated 20 18 [70]Open in a new tab ^a Population used in the analyses of gene and miRNA expression. 2.3. Sample preparation, cDNA and small RNA library construction and sequencing Based on chromosome survey on the above 43 populations, diploid (population code: CG) and hexaploid cytotypes (HS 3) were selected as the experimental materials for comparison to investigate gene and miRNA expression profiling in this work (Figure [71]2a,b). Leaves were collected from three independent comparable potted‐seedlings creating three biological replicates for diploid (D) and hexaploid cytotypes (H). The top three to four fully expanded leaves were gently removed in the morning, covered by aluminum foil, frozen in liquid nitrogen immediately, and subsequently stored at −80°C until total RNA extraction. The cDNA and small RNA library were constructed following the methods provided by Beijing Genomics Institute (BGI, Shenzhen, China) (Supporting Information Appendix [72]S1). Figure 2. Figure 2 [73]Open in a new tab Morphological and cytological divergences between diploid (sampled from CG population, 2n = 2x = 18) and hexaploid cytotypes (HS 3 population, 2n = 6x = 54) of S. canadensis. Plant morphology of diploid and hexaploid cytotypes in vegetative stage (a) and reproductive stage (b). Chromosome numbers of diploid (c) and hexaploid cytotypes (d). Population names follow Table [74]1 2.4. De novo assembly and unigene annotation Raw reads were filtered by removing those reads that contained adaptors, unknown nucleotides (more than 5%), and low‐quality bases (more than 20% of the bases with a quality score less than 15). De novo assembly of all processed reads was performed by Trinity (version: v2.0.6, Grabherr et al., [75]2011), with parameters set as follows: ‐min_contig_length 200; ‐CPU 8; ‐min_kmer_cov_4; ‐min_glue_4; ‐bfly_opts'‐ V5; ‐edge‐thr=0.1; and ‐stderr'. Then, the constructed transcripts from the Inchworm, Chrysalis, and Butterfly modules of Trinity were further clustered into nonredundant unigenes by using TGICL (version: v2.0.6, Pertea et al., [76]2003) to eliminate the redundant Trinity‐generated transcripts, with parameters set as follows: ‐I40‐c10‐v25‐O'‐repeat_stringency 0.95‐minmatch 35‐minscore35'. To construct a uniform transcriptome reference, all assembled unigenes from six samples of two cytotypes of S. canadensis were pooled together and further clustered to generate “All‐Unigene” for subsequent assembly evaluation, unigene annotation, and expression analysis. The “All‐Unigene” sequences were aligned by BLASTx to a series of protein databases to gain unigene annotation. See Supporting Information Appendix [77]S1 for more unigene annotation details. 2.5. Unigene quantification and differentially expressed unigene (DEU) analysis The Bowtie2 program (version: v2.2.5, Langmead & Salzberg, [78]2012) was used to map clean reads from each sample to assemble “All‐Unigene” with the following parameters: ‐q; ‐phred 64; ‐sensitive; ‐dpad 0; ‐gbar 99999999; ‐mp 1,1; ‐np 1; ‐score‐minL,0, ‐0.1‐I1‐X 1000; ‐no‐mixed; ‐no‐discordant; ‐p 1‐k 200, and RSEM (version: v1.2.12, Li & Dewey, [79]2011) was applied to calculate the read counts mapped to each unigene with the default parameter. Then, fragments per kilobase of transcript per million fragments mapped (FPKM) was applied to normalize the expression value. Differential gene expression analysis was performed using the DESeq2 R package as described by Love, Huber, and Anders ([80]2014) for comparisons between diploid and hexaploid cytotypes with three biological replicates. An absolute value of log[2]fold‐change ≥2 and an adjusted p‐value <0.001 was set as the threshold to identify DEUs. Following this, identified DEUs were subjected to GO and KEGG analyses. The regulated unigenes were assigned GO terms by the Blast2‐GO program (version: v2.5.0, Conesa et al., [81]2005), and their enrichment was performed for testing over‐represented GO categories using the GOseq R package with a corrected p‐value (FDR analog) setting of ≤0.05. DEUs were further assigned KO (KEGG Orthology) numbers using the KEGG database, and their enrichment was performed as mentioned for GO. 2.6. Transcription factor (TF)‐encoding gene prediction To identify putative TF candidate genes, getorf (version: EMBOSS: 6.5.7.0, Rice, Longden, & Bleasby, [82]2000) was used to find and extract open reading frames (ORFs) from all assembled unigene sequences with the minimum size parameter set as 150, and then the sequences of ORFs were searched against the plant transcription factor database (PlnTFDB; version: 3.0) using hmmsearch (version: 3.0, Mistry, Finn, Eddy, Bateman, & Punta, [83]2013) with the default parameters. 2.7. miRNA identification and differentially expressed miRNA (DEM) analysis Raw reads were filtered by removing low‐quality contaminated reads as well as adaptor sequences, and then generated clean reads in the range of 18–30 nt were chosen for mapping to the S. canadensis mRNA transcriptome by SOAP with default settings. Subsequently, sequences with a perfect match were compared to Rfam 11.0 and NCBI GenBank databases to eliminate noncoding RNAs, including rRNA, scRNA, snRNA, snoRNA, tRNA, and repeats. Given that sequences from S. canadensis were not included in miRBase, the remaining unique reads were searched against currently annotated plant miRNAs (Viridiplantae) available in the miRBase 22.0 database using the BLASTn program to identify the known miRNAs. Transcripts per million was used to normalize the read count of each identified miRNA based on the following formula: Normalized expression = Actual miRNA count × 10^6/Total count of clean reads. After normalization, differential expression analysis of miRNA was performed using DEGseq as described by Wang, Feng, Wang, Wang, and Zhang ([84]2010). An absolute value of log[2]fold‐change ≥1 and a q‐value <0.001 was set as the threshold to identify DEMs. 2.8. Prediction of miRNA targets To predict the potential genes targeted by miRNAs, the Targetfinder (version: 1.5, Fahlgren & Carrington, [85]2010) in combination with psRobot (version: 1.2, Wu, Ma, Chen, Wang, & Wang, [86]2012) software was applied to predict as many miRNA targets as possible from the assembled S. canadensis unigene set (116,801 “All‐Unigene”) with default parameters. Additionally, the expression level of predicted miRNA targets was taken from the inventory of assembled “All‐Unigene.” GO terms were also evaluated using a similar method. 2.9. Visualization of miRNA‐target interaction network To unravel complex links between candidate miRNAs and unigenes, we proposed a strategy that integrated expression data of DEMs and DEUs to visualize the miRNA‐target interaction network and further discover key miRNAs. Here, we defined coherent miRNA targets as those presenting opposite expression patterns compared with those of the miRNAs, showing that the expression of unigenes was negatively correlated with that of miRNAs (Ye, Wang, & Wang, [87]2016). To construct the miRNA‐target interaction network, three separate steps were performed. First, DEMs and DEUs were screened following the method mentioned above. Second, predicted targets of up‐regulated miRNAs (down‐regulated miRNAs) overlapped with identified down‐regulated unigenes (up‐regulated unigenes) to obtain coherent miRNA targets. Finally, acquired coherent miRNA targets and DEMs were subjected to visualization of the miRNA‐target interaction network by Cytoscape. 2.10. Candidate unigene and miRNA validation via qRT‐PCR Eighteen promising candidate unigenes and six miRNAs observed to be differentially expressed were chosen for qRT‐PCR to validate the reliability of RNA‐seq and sRNA‐seq results with the following selection criteria: (a) up‐ or down‐regulated unigenes discussed in this paper (i.e., Expansin, ARGOS); and (b) miRNA‐target interaction pairs that were negatively correlated in expression levels. qRT‐PCR was implemented in triplicate on an ABI Step One Plus Real‐Time PCR System (Applied Biosystems) with unigene‐ and miRNA‐specific sense and anti‐sense primer (Table [88]A1 in Appendix, Supporting Information Appendix [89]S1). A homolog of GAPDH (Unigene25510_All) was co‐amplified to normalize the expression values of unigenes and miRNAs in each sample using the double‐standard curve method. 3. RESULTS 3.1. Gene expression profiling in diploid and hexaploid cytotypes of S. canadensis The inspection of chromosome numbers revealed that two cytotypes were ascertained among the 449 individuals of S. canadensis examined. For the cultivated population, all individuals were observed to be diploid cytotypes with a chromosome number of 2n = 2x = 18 (Figure [90]2c). For the invasive populations, all individuals were observed to be hexaploid cytotypes with a chromosome number of 2n = 6x = 54 (Table [91]1; Figure [92]2d). However, tetraploid cytotypes with a chromosome number of 2n = 4x = 36 or mixed‐cytotypes reported by Li ([93]2011) were not found in the current work. To explore key candidate genes behind the invasiveness differences in diploid and hexaploid cytotypes, we generated the first transcriptomic profile of S. canadensis. A total of 334.79 million (M) raw reads were produced and subjected to Seq‐QC collating, which resulted in 289.45 M (86%) clean reads with Q20 values ranging from 98.88% to 98.94%. Then, clean reads from six libraries were de novo assembled separately into unigenes by Trinity. These assembled unigenes were pooled together and further clustered into a reference transcriptome (116,801 “All‐Unigene”) with an average length of 1,056 bp, a N50 value of 1,610 bp, and a GC content of 39.20% (Table [94]A2 in Appendix). These numbers are comparable to those generated in other polyploid studies (e.g., Vigna et al., [95]2016; Zhou et al., [96]2015) and imply a high‐quality assembly. Furthermore, we also found that the length distribution of the assembled “All‐Unigene” ranged from 224 to 23,608 bp with a total length of 123,376,557 bp, of which 32,942 (28.20%) unigenes ranged from 300 to 500 bp, 32,696 (27.99%) unigenes ranged from 500 to 1,000 bp, 33,468 (28.65%) unigenes ranged from 1,000 to 2,000 bp, and 17,695 (15.15%) unigenes had lengths longer than 2,000 bp (Figure [97]A2 in Appendix). Out of the 116,801 “All‐Unigene” acquired above, expression of 12,428 unigenes was found only in diploid cytotypes, and expression of 19,520 unigenes was observed only in hexaploid cytotypes. These seem to represent a suite of ploidy‐dependent unigenes, which means a specific role of these ploidy‐dependent unigenes in contrasting invasiveness differences. Subsequently, to identify notably changed unigenes, we applied the aforementioned filter criterion and noticed that 12,897 unigenes displayed at least a four‐fold change in expression levels, with the majority of them (6,768 out of 12,897) down‐regulated in hexaploid cytotypes (Supporting Information Table [98]S1). After that, these DEUs were further subjected to investigation of the specific regulated pathways in which they were involved. However, it must be underlined here that our work has revealed novel unigenes whose functions are unknown, which will be the long‐running theme of future research. Furthermore, qRT‐PCR analysis performed for eighteen DEUs confirmed the mRNA changes detected by RNA‐seq (Figure [99]A3 in Appendix). Further, the identified 2,644 putative TF‐encoding genes in this work were assigned to 58 TF families, of which MYB members (337) were over‐represented, followed by MYB‐related (270), AP2‐EREBP (211), bHLH (134), and WRKY family members (120). In addition, we discovered 381 TF‐encoding genes that were differentially expressed, with the majority being up‐regulated in hexaploid cytotypes. Notably, among these differentially expressed TF‐encoding genes, almost all the members of the bHLH group were found to be up‐regulated (12/15 genes) in hexaploid cytotypes. In addition, MYB (32/53), MYB‐related (31/47), ARF (10/16) and Trihelix (9/12) members exhibited a similar trend. However, the majority of TF‐encoding genes belonging to the WRKY (18/26) and NAC (13/17) families were down‐regulated in hexaploid cytotypes (Figure [100]3). These TF genes had differential expression patterns, implying a variety of regulatory modes. Figure 3. Figure 3 [101]Open in a new tab A bar graph representing the differential expression of TF‐encoding genes in diploid and hexaploid cytotypes of S. canadensis. Yellow indicates the up‐regulated genes and blue down‐regulated genes in hexaploid cytotypes 3.2. Functional and pathway analysis of ploidy‐responsive unigenes in S. canadensis To better understand the functionality of unigenes differentially expressed in response to ploidy, we mapped the above‐mentioned DEUs to the GO and KEGG databases to perform functional analyses and found that a total of 4,545 (35.24%) unigenes from the 12,897 DEUs were successfully classified into three major functional categories: biological process (3,171), molecular function (3,536), and cellular component (2,764). Then, the three major categories were further assigned to 50 terms (Figure [102]A4a in Appendix), including 21 terms in the biological process, 14 terms in the molecular function, and 15 terms in the cellular component categories. The most abundant GO term related to biological process was “metabolic process” represented by 2,374 DEUs, followed by “cellular process,” “single‐organism process,” “localization,” “biological regulation,” and “response to stimulus” represented by 2,243, 1,741, 568, 501, and 462 DEUs, respectively. In the molecular function category, the two main representative distributions were “catalytic activity” (2,501) and “binding” (1,904). Other GO terms, such as “transporter activity,” “enzyme regulator activity,” “structural molecule activity,” “antioxidant activity,” “receptor activity” and “channel regulator activity,” associated with the biosynthesis of secondary metabolites were also enriched. With respect to the cellular component category, a large proportion of DEUs were clustered in “cell” (1,824), “cell part” (1,808), “membrane” (1,516), and “organelle” (1,264). In addition, we also found that 8,666 DEUs were assigned to 133 unique KEGG pathways, with 5,705 representing metabolism pathways, 2,006 pathways in genetic information processing, 469 pathways in cellular process, 418 pathways in environmental information processing, and 415 pathways in organismal systems (Figure [103]A5 in Appendix). Notably, there were only two pathways that were significantly over‐represented under “organismal systems,” that is, “circadian rhythm‐plant” and “plant‐pathogen interaction.” The most represented pathway in DEUs was “metabolic pathways,” followed by “biosynthesis of secondary metabolites,” “plant‐pathogen interaction,” “RNA transport,” and “spliceosome.” Subsequently, the hypergeometric distribution was calculated to identify significantly enriched pathways in which DEUs were involved. A total of eight pathways associated with metabolism were significantly enriched, with a Q value ≤0.05 (Table [104]A3 in Appendix). It was conspicuous that unigenes related to “metabolic pathways (Pathway ID: ko01100)” were significantly enriched among the DEUs, implying that they may operate in the metabolic adaptation mechanism of hexaploid cytotypes. Additionally, unigenes for carbohydrate metabolism of “pentose and glucuronate interconversions (Pathway ID: ko00040)” were enriched. Moreover, unigenes for lipid metabolism of “fatty acid degradation (Pathway ID: ko00071)” were enriched. Additionally, DEUs involved in the metabolism of terpenoids and polyketides, particularly “carotenoid biosynthesis (Pathway ID: ko00906),” and “sesquiterpenoid and triterpenoid biosynthesis (Pathway ID: ko00909)” were enriched. Finally, “biosynthesis of secondary metabolites (Pathway ID: ko01110),” “isoflavonoid biosynthesis (Pathway ID: ko00943),” and “flavone and flavonol biosynthesis (Pathway ID: ko00944)” were enriched, signifying considerable modulation of unigenes responsible for the regulation of plant secondary metabolites. 3.3. miRNA expression profiling in diploid and hexaploid cytotypes of S. canadensis A total of 179.9 M 50‐base pair (bp) single‐end raw reads were produced and subjected to Seq‐QC collating, which resulted in 166.6 M (92.6%) clean reads with lengths ranging from 18 to 30 nt (Table [105]A4 in Appendix). The sRNA length distribution in six libraries showed that the majority of reads were distributed between 20 and 24 nt in length, which corresponds to the size from Dicer‐like digestion products. In addition, the most abundant sequence in all six libraries was 24 nt sRNA (average 37.24% vs. 42.31% in D vs. H), followed by 21 nt sRNA (average 20.06% vs. 23.11% in D vs. H) (Figure [106]A6 in Appendix), which was in agreement with the typical size distribution of sRNAs reported in other plant species, such as Arabidopsis (Rajagopalan, Vaucheret, Trejo, & Bartel, [107]2006), Oryza sativa (Morin et al., [108]2008), and Citrus trifoliata (Song et al., [109]2010). We identified 186 miRNAs belonging to 44 miRNA families in two cytotypes of S. canadensis and found that the identified families included a changing count of miRNA members (Supporting Information Table [110]S2). Among the detected miRNAs, the miR166 family possessed the largest number of members, with 26 members that were discriminated by the divergences in nucleotide sequences, followed by miR171, miR167, miR168, miR396, miR156, miR169, miR159, miR319, miR164, miR393, and miR160 families, with 14, 12, 11, 11, 10, 10, 8, 8, 6, 6, and 5 members, respectively. miR398, miR399, and miR858 included four members, and miR390, miR395, and miR403 included three members. Of the remaining 26 miRNA families, 12 families, such as miR157, miR161, and miR162 families, comprised two members, and 14 miRNA families were represented only by a single member each. A further analysis showed that 59 miRNAs were differentially expressed, of which 38 miRNAs were up‐regulated and 21 miRNAs were down‐regulated in hexaploid cytotypes relative to their diploids. Among the DEMs, sca‐miR395c, sca‐miR8155, and sca‐miR6173 were markedly down‐regulated with log[2]fold‐change values of −3.23 (q = 1.53e−42), −3.15 (q = 1.40e−04) and −2.45 (q = 6.35e−11), respectively, and sca‐miR166p, sca‐miR528, and sca‐miR396a were markedly up‐regulated with log[2]fold‐change values of 5.16 (q = 7.39e−12), 5.13 (q = 9.36e−07), and 5.04 (q = 5.30e−11), respectively. Notably, for the miR160 and miR169 family, sca‐miR160e, sca‐miR169b, sca‐miR169e, sca‐miR169f, sca‐miR169g, and sca‐miR169h were up‐regulated specifically in hexaploid cytotypes, while sca‐miR160b and sca‐miR169d were up‐regulated specifically in diploid cytotypes. These observations suggested that different members from the same miRNA family had different regulatory modes, probably associated with the cooperative and redundant regulation activity of miRNAs. qRT‐PCR analysis performed for six DEMs confirmed the miRNA changes detected by sRNA‐seq (Figure [111]A7 in Appendix). In addition, correlation between qRT‐PCR results and sequencing results were also calculated. We acquired a significant Pearson “r” close to 0.85 (p < 0.001) (Figure [112]A8 in Appendix), which strongly suggested that our transcriptome and sRNA sequencing data were credible. 3.4. Unigenes involved in growth‐related pathways are targeted by DEMs Our analysis revealed 1,801 unigenes from 116,801 assembled S. canadensis “All‐Unigene” were predicted as targets of 184 miRNAs, of which 884 putative targets were predicted to be cleaved by 58 DEMs. Moreover, a meticulous inspection of the DEMs and their corresponding targets indicated that (a) miR5139a had the highest target abundance (179), and (b) the genes such as CL10163.Contig1_All, CL13112.Contig1_All, and Unigene2861_All had the highest miRNA abundance (4). To understand in depth the group of unigenes targeted by DEMs, GO functional analysis of the predicted targets was carried out. Under the biological process category of GO classification, unigenes involved in terms such as “cellular process,” “metabolic process,” “single‐organism process,” “response to stimulus,” etc. were abundantly enriched as the targets of DEMs. Under the molecular function category, unigenes displaying “catalytic activity,” “binding,” “transporter activity,” etc. were targeted by miRNAs. Moreover, “cell,” “membrane,” “organelle,” etc. related unigenes were discovered to be clustered into the cellular component category as targets (Figure [113]A4b in Appendix). In further pathway analysis of 884 putative targets, “cutin, suberin and wax biosynthesis (Pathway ID: ko00073),” “protein processing in endoplasmic reticulum (Pathway ID: ko04141),” “plant hormone signal transduction (Pathway ID: ko04075),” “selenocompound metabolism (Pathway ID: ko00450)” and “cysteine and methionine metabolism (Pathway ID: ko00270)” pathways were significantly enriched with a Q value ≤0.05. These results suggest that miRNAs were more likely to activate plant primary metabolism and make contributions to the improved vigor shown by hexaploid cytotypes, as it has been noted earlier that hexaploid cytotypes typically exhibited enhanced growth in comparison with diploids. 3.5. Integrative analysis of gene and miRNA expression confirms that environmental adaptation‐related unigenes are centrally targeted To detect which biological processes or pathways within a cell were most likely regulated by miRNAs, we integrated overall gene and miRNA expression data to identify miRNA‐target interacting pairs that were negatively correlated in log[2]fold‐change between DEMs and target mRNA expression. As a result, 83 miRNA‐target interacting pairs with the involvement of 24 DEMs and 69 targets were visualized by Cytoscape. For each such pair, we then classified 83 miRNA‐target interacting pairs into two categories depending on the expression patterns of DEMs as either up‐regulated or down‐regulated, respectively, for 47 miRNA‐target pairs involved in 10 down‐regulated miRNAs and 34 up‐regulated targets; or 36 miRNA‐target pairs involved in 14 up‐regulated miRNA and 35 down‐regulated targets (Figure [114]4). Furthermore, we have also taken note that the coherent miRNA targets included (a) several TFs that were predicted to be targeted by miRNA regulators, for example, sca‐miR164d targets FAR1, sca‐miR530 targets MYB, sca‐miR396a targets Trihelix, and sca‐miR5139a targets VOZ1‐like, suggesting that these miRNAs may operate to enhance the adaptation of hexaploid cytotypes through an integrative miRNA‐TF‐mRNA regulatory network; (b) receptor‐like protein kinases (RLKs) that were predicted to be targets of multiple miRNAs such as sca‐miR161a, sca‐miR5139a, sca‐miR5139b, and sca‐miR8155. This target is an important enzyme gene and functions in regulating plant growth, development, signal transduction, immunity, and stress responses (Sun, Li, Wang, Zhang, & Wu, [115]2017). Notably, these RLK genes were remarkably up‐regulated in hexaploid cytotypes, suggesting that their regulator miRNAs may play key roles in the environmental adaptation of hexaploid cytotypes; (c) unigenes associated with methylation and ubiquitination processes, such as histone‐lysine N‐methyltransferase (CL7649.Contig3_All), ubiquitin‐protein ligase (CL2235.Contig13_All), U‐box domain‐containing protein (CL2207.Contig4_All), and F‐box protein (Unigene1223_All), that were predicted to be targets of sca‐miR396d, sca‐miR444a, sca‐miR393d, and sca‐miR5139a, suggesting that these unigenes may be subjected to miRNA‐mediated DNA methylation and ubiquitination; and (d) two dirigent protein genes (CL6884.Contig1_All and CL6884.Contig3_All) that were predicted to be targets of sca‐miR169d. This target was an unspecific oxidizing enzyme gene for radical formation that functions in lignan biosynthesis, which was previously reported to be an integral regulator of plant secondary metabolism (Effenberger et al., [116]2015). Notably, these dirigent protein genes were remarkably up‐regulated in hexaploid cytotypes, suggesting that its regulator sca‐miR169d may play key roles in plant secondary metabolism (Table [117]A5 in Appendix). Figure 4. Figure 4 [118]Open in a new tab miRNA‐gene interaction network of S. canadensis. In this network, oval nodes represented unigenes and triangle nodes represented miRNAs. The negative correlation was denoted by a line. The yellow and blue color mean up‐regulation and down‐regulation and the highest to lowest fold changes are marked from yellow to blue 4. DISCUSSION A large number of works have investigated ecological and evolutionary elements responsible for successful invasion (Hahn et al., [119]2012; Thébault et al., [120]2011). However, research into the molecular basis for invasiveness in invasive plants is just getting started. Here, we found 12,897 unigenes and 59 miRNA regulators with divergences in expression between diploid and hexaploid cytotypes. Intriguingly, among them were an over‐representation of unigenes and coherent miRNA targets relevant to metabolism, plant growth and development, and stress responses, implying that these modified genetic and epigenetic attributes may harbor both biochemical and ecological advantages that were beneficial to the successful invasion of hexaploid cytotypes. 4.1. Unique gene and miRNA expression characteristics might have contributed to the invasiveness of hexaploid cytotypes In Arabidopsis thaliana, only 0.1% differences in gene expression between diploid and autotetraploid were detected (Yu et al., [121]2010). In newly synthesized autotetraploid Paspalum notatum, 0.6% of genes were differentially expressed compared to its diploid (Martelotto et al., [122]2005). Similarly, the analysis of 21,081 genes in Citrus limonia autotetraploids revealed less than 1.1% differences in comparison with diploids (Allario et al., [123]2011). In contrast, many researchers have observed a more noticeable transcriptomic divergence between allopolyploids and their parents in several plants (Li et al., [124]2014; Ye et al., [125]2016). Remarkably, here we detected >10% transcriptomic differences as a consequence of hexaploid cytotype formation. Two factors may account for this dramatic change. First, S. canadensis is a polyploid, and assembling its transcriptome has been exceedingly difficult because it principally comprises highly similar repeats, thereby causing several contigs that often represent nonoverlapping fragments of the same unigene. Second, given that polyploid effects on gene expression might be induced by genome doubling and/or hybridization, we speculate that the expression pattern of hexaploid cytotypes of S. canadensis should be that of an allohexaploid. Furthermore, a significant caveat in the interpretation of these results is that we have only sequenced one population per ploidy level, and genetic differentiation among different geographic populations of the same ploidy could also be contributing to gene and miRNA expression differences. Thus, further work is essential to explore the genetic relationship between cytotypes. Furthermore, to test whether the invasiveness difference between plants was reflected in the changes of their gene expression patterns, we examined evidence for successful invasion in introduced populations across multiple invasive plants. Here, common‐garden studies comparing native and introduced populations of Cirsium arvense, Centaurea diffusa, and Mikania micrantha, as well as comparisons between S. canadensis and invasive taxa of the Asteraceae, have been performed, and obvious similarities have emerged (Guggisberg, Lai, Huang, & Rieseberg, [126]2013; Guo et al., [127]2018; Hodgins et al., [128]2015; Turner, Nurkowski, & Rieseberg, [129]2017). In these studies, introduced populations notably differ from their native populations with regard to stress response. In line with this observation, the significantly different regulation of stress response genes, such as receptor‐like proteins, is of particular interest between introduced and native populations because these genes mediate plant cellular defense pathways. Similar, several stress response genes associated with secondary metabolism, such as the cytochrome P450 gene family, were also found to be significantly expressed in the present work. Therefore, it is reasonable to speculate that these stress response genes might have crucial functions in invasive characteristics. However, we also noticed that genes involved in photosynthesis were exclusively enriched in M. micrantha (Guo et al., [130]2018). Hence, it seems that the pattern of gene expression across different invasive plants is dependent on a specific plant, and thus, it is difficult to generalize a rule of gene expression during invasion. Likewise, the relative amount of miRNAs was higher in a derived hexaploid wheat (BBAADD) than in the parental tetraploid Triticum turgidum ssp. durum (BBAA) and diploid Aegilops tauschii (DD) (Kenan‐Eichler et al., [131]2011). Analogously, the number of miRNA or miRNA families in cultivated allotetraploid cotton G. hirsutum (AADD) was markedly greater than those in its two diploid ancestors, G. raimondii (DD), and G. arboreum (AA) (Xie & Zhang, [132]2015). Ghani et al. ([133]2014) also reported that the percentages and expression levels of miRNAs increased in allodiploid (AB) and allotetraploid (AABB) relative to the parents Brassica rapa (AA) and Brassica nigra (BB). In the present work, the number and expression levels of miRNAs in hexaploid cytotypes were greater than those in their diploids, which was consistent with the findings of the above‐mentioned studies. These results suggest that an increase in ploidy was generally coupled with an obvious increase in the percentages and expression levels of miRNAs. 4.2. Several regulatory mechanisms seem to operate gene expression properly in hexaploid cytotypes How does hexaploid cytotypes regulate the differential expression of unigenes? Several mechanisms could be associated with this regulation. miRNAs work as regulators for controlling target‐gene expression, thereby affecting a variety of aspects of phenotype, growth, development, and stress response (Ha et al., [134]2009). Here, we showed a subset of key candidate miRNA regulators within diploid and hexaploid cytotypes and used these DEMs to predict putative targets using two different target‐prediction software. To shed light on the regulatory action of these DEMs, we compared these predicted targets with DEUs based on GO functional classification (Figure [135]A4a,b in Appendix) and found that biological processes were highly likely to be regulated by miRNAs, such as (a) a considerable proportion of the enriched unigenes were clustered in “biological process”; (b) processes such as “metabolic” and “cellular” were abundantly enriched; and (c) “single‐organism,” “localization,” “biological regulation,” and “response to stimulus” were also adequately reflected. Such observations suggested that unigenes described in the above‐mentioned terms were most likely targeted by miRNAs. However, unigenes associated with “cell killing,” “locomotion,” and “rhythmic process” were enriched only in DEUs, implying that although unigenes associated with the foregoing processes were differentially expressed, this regulation of gene expression cannot be attributed to the miRNA‐induced cleavage of targets. In contrast, no term was only enriched under the same category for the targets of DEMs. These observations suggested that few specific biological processes were regulated by miRNAs. Similarly, “channel regulator activity” and “electron carrier activity” enriched in the “molecular function” category were only amid DEUs. In addition, under the category of “cellular component,” unigenes related to terms “cell,” “membrane” and “organelle” were overwhelming in this comparison, whereas no unigene associated with terms “nucleoid,” “virion,” and “virion part” was enriched among the targets of DEMs, indicating that these unigenes are closely regulated at the transcriptional level and may not be prominently influenced by miRNA‐induced gene silencing. In addition to miRNAs, other accessional regulation manners, such as DNA methylation, may also function to regulate gene expression. There is impressive evidence that an allopolyploid's intergenomic interactions between two divergent genomes were projected to incur DNA methylation changes, eventually causing the differential expression of genes, which can potentially lead to profound phenotypic consequences (Chen, [136]2007; Salmon & Ainouche, [137]2010). DNA methylation changes between an allopolyploid and its parents have been very well reported. For instance, in Spartina allopolyploids, a high level of epigenetic regulation might explain the morphological plasticity and its larger ecological amplitude (Salmon, Ainouche, & Wendel, [138]2005). Additionally, Madlung et al. ([139]2002) reported that changes in DNA methylation would result in the development of altered morphologies in synthetic allotetraploids. Although DNA methylation alterations are principally observed in allopolyploids, activation, or repression of gene expression has also been shown to correlate with DNA methylation variation in autopolyploid Arabidopsis (Yu et al., [140]2010) and Cymbopogon (Lavania et al., [141]2012). In the present work, a large number of DEUs related to epigenetic regulation were investigated in two cytotypes of S. canadensis (Supporting Information Table [142]S3). In particular, transcriptome analysis defined eleven unigenes (annotated as DNA (cytosine‐5)‐methyltransferase1 gene, for example, CL12526.Contig3_All, CL5231.Contig1_All, and CL5231.Contig2_All) that displayed striking changes in gene expression. Interestingly, almost all the DNA (cytosine‐5)‐methyltransferase1 genes were found to be significantly down‐regulated (10/11 genes) in hexaploid cytotypes, and these observations should be further investigated because such genes could potentially participate in the maintenance of CG methylation. Additionally, to answer developmental and environmental alterations, chromatin composed of DNA and histones in eukaryotic cell nuclei is modulated by several histone modifications. Among these modifications, histone demethylation regulates gene expression mainly by demethylating histone lysine residues (Shi & Tsukada, [143]2013). Recent studies have identified Jumonji (JmJ) proteins to be involved in histone demethylation and closely related to the reproductive process. The loss‐of‐function mutations of the rice gene JmJ706 resulted in spikelet development defects (Sun & Zhou, [144]2008). Here, a total of eight DEUs (e.g., CL1566.Contig6_All, CL1566.Contig9_All, and CL1566.Contig12_All) were annotated as JmJ genes, and all of them were up‐regulated in hexaploid cytotypes, which might partly suggest that the JmJ histone demethylase unigenes may alter the expression of a large number of target genes and contribute to the variation in physiology, biochemistry and phenotype between diploid and hexaploid cytotypes. However, further study is needed. Taken together, these data clearly state that complicated and overlapping gene expression regulatory mechanisms may have evolved in hexaploid cytotypes to guarantee suitable transcriptional control in response to environmental stimuli. 4.3. Potential roles of transcriptional alterations in the successful invasion of hexaploid cytotypes Polyploids play recognized roles in driving organ size and growth of plants. The leaf is the main photosynthetic organ, and its size strongly affects the energy capture, photosynthetic capacity, and physiological activities of plants (Baute et al., [145]2017; Niinemets, Portsmuth, & Tobias, [146]2006). The coordination of cell proliferation and expansion is a crucial determinant that serves a critical function in precisely controlling leaf size and growth caused by cell ploidy (Baute et al., [147]2017; Marshall et al., [148]2012; Sugiyama, [149]2005), which have been previously suggested to be regulated by a number of genes encoding transcription factors, modification proteins, plant hormones, and cell wall protein. The Growth‐Regulating‐Factor (GRF) protein, a plant‐specific transcription factor, has been confirmed to affect leaf growth by positively regulating cell proliferation, cell expansion, and adaxial‐abaxial patterning (Omidbakhshfard, Proost, Fujikura, & Mueller‐Roeber, [150]2015). In addition, the GRF protein has also been shown to perform transcription regulation functions by interacting with GRF‐Interacting Factor (GIF) protein (Debernardi et al., [151]2014). In this work, five unigenes (unigene52239_All, CL15674.Contig2_All, CL8349.Contig2_All, CL1635.Contig4_All, and CL1635.Contig1_All) annotated as GRF and two unigenes (Unigene15938_All and CL6058.Contig2_All) annotated as GIF were found to be differentially expressed and may form functional complexes potentially implicated in leaf size and growth. Furthermore, other transcription factors, such as the TCP transcription factor (e.g., CL7816.Contig2_All, CL9341.Contig3_All), were also identified as regulators of leaf size. Except for TFs, regulatory proteins act as important regulators of leaf size and growth by influencing cell proliferation. EBP1, an ortholog of ErbB3‐binding protein from humans, regulates leaf size and growth by cell proliferation. Some studies highlight that the expression of EBP1 correlates with plant organ size, growth, and stress tolerance (Cao et al., [152]2009; Horváth et al., [153]2006). In the present work, one ortholog of EBP1, CL16506.Contig4_All, was found to be significantly up‐regulated in hexaploid cytotypes, which may be responsible for the larger leaves, faster growth, and better stress resistance of hexaploid cytotypes. Similarly, F‐box proteins, which are members of regulatory protein families that affects leaf size (Baute et al., [154]2017), are abundantly expressed. Moreover, earlier studies showed that auxin mediated the expression of multiple genes (e.g., ARGOS and ARF) to affect plant organ size and growth (Schruff et al., [155]2006; Wang, Zhou, Xu, & Gao, [156]2010). Auxin‐Regulated Gene involved in Organ Size (ARGOS), a gene deeply induced by auxin, participate in organ size regulation. Wang, Zhou, et al. ([157]2010) also pointed out that overexpression of a Chinese cabbage (Brassica rapa) BrARGOS gene in Arabidopsis elevates the size of plant organs. In this work, an ortholog of BrARGOS, CL15040.Contig2_All was detected, and the up‐regulated expression may have similar functions in the organ giantism observed in hexaploid cytotypes. Auxin Response Factor (ARF), a transcription factor, functions in plant size, growth, and stress adaptation by transcriptionally activating and repressing the expression of auxin response genes (Zhao, Zhang, Ma, & Wang, [158]2016). Here, 16 ARF encoding genes were found to be differentially expressed, which might contribute to invasiveness differences between diploid and hexaploid cytotypes. Lastly, abundant studies have shown correlations between expansin gene expression and cell wall remodeling, growth and stress response, and phenotype changes in plants (Goh, Sloan, Malinowski, & Fleming, [159]2014; Lee & Choi, [160]2005; Li et al., [161]2013), which supports the roles for expansin as an important cell wall protein in plant cell wall modification, growth promotion, and stress tolerance. As expected, the overexpression of expansin genes has remodeled leaf structure, which confers them enhanced tolerance to abiotic stresses (Cho & Cosgrove, [162]2000; Kwon et al., [163]2008). In the present work, thirteen unigenes encoding expansins were differentially expressed, and eleven of them were more highly expressed in hexaploid cytotypes of S. canadensis than in diploid cytotypes. These results suggested that the activation of expansins may be a rapid growth and adaptation mechanism of hexaploid cytotypes in novel heterogeneous environments. Furthermore, polyploids can also profoundly affect plant metabolism qualitatively and quantitatively, furnishing the chance for increased metabolic activity through transcriptional divergence, which eventually results in alterations in the levels of secondary metabolites (Fasano et al., [164]2016). There are multiple studies on the induction of polyploids to promote the production of specific secondary metabolites. For instance, autotetraploids of Catharanthus roseus produced more vindoline, catharanthine, and vinblastine than their diploids (Xing et al., [165]2011). Echinacea purpurea autotetraploids showed that the induction of polyploids resulted in higher caffeic acid derivatives and alkamides (Xu et al., [166]2014). Evidence of the influence of polyploids on chemical profiles has also been recorded in allopolyploids. Banyai et al. ([167]2010) reported that allotetraploid Artemisia annua produced more terpenoids or triterpene‐type compounds than diploids. Supporting the role of plant secondary metabolism in polyploid‐mediated invasiveness differences “biosynthesis of secondary metabolites (Pathway ID: ko01110)” was found to be the most significantly enriched pathway with a Q value far below 0.05 in the pathway enrichment analysis of DEUs in the present work. Taking the above into account, we propose that polyploids are more likely to remodel the transcriptome and metabolome in hexaploid cytotypes, resulting in ploidy‐specific metabolic adaptation. Moreover, a marked number of DEUs encoding enzymes related to plant metabolism were observed, which further supports this plausible explanation. The synthesis of secondary metabolites primarily contains the oxidation, reduction, and cyclization steps, in which unigenes encoding enzymes of cytochrome P450 (CYPs) and uridine diphosphate glucuronosyl transferases (UGTs) play crucial roles in catalyzing these reactions (Zhang et al., [168]2016). Based on the functional annotation of DEUs, a total of 120 core enzyme unigenes encoding CYPs were differentially expressed (Table [169]A6 in Appendix). In addition, CYPs are one of the largest superfamilies of enzyme proteins (Darabi, Seddigh, & Abarshahr, [170]2017). A large number of CYPs are involved in a wide range of biosynthetic reactions and biochemical pathways, leading to the synthesis of UV protectants (flavonoids and anthocyanins), defensive compounds (isoflavonoids, phytoalexins, hydroxamic acids, and terpenes), fatty acids, hormones (gibberellins and brassinosteroids), signaling molecules (oxylipins, salicylic acid, and jasmonic acid), accessory pigments (carotenoids), and structural polymers such as lignins (Darabi et al., [171]2017; Schuler & Werck‐Reichhart, [172]2003). In the present work, many CYP‐related unigenes were identified, such as CYP93A (e.g., CL361.Contig6_All, CL1330.Contig7_All, CL16738.Contig1_All), CYP76B (e.g., CL3689.Contig2_All, CL1852.Contig2_All), and CYP71 (e.g., CL6714.Contig2_All, Unigene23159_All, CL15279.Contig1_All), which respectively participated in the biosynthesis of isoflavonoids, flavonoids, and sesquiterpenoids and triterpenoids, which may act as defensive compounds that protect against oxidative damage under abiotic stress. Furthermore, a great deal of UGT genes that participated in flavonoid biosynthesis, such as UGT73, UGT74, UGT76, UGT83, UGT85, and UGT89, were also identified (Table [173]A7 in Appendix). Given these findings, it is attractive to investigate the potential model whereby polyploids impact the metabolome in hexaploid cytotypes of S. canadensis. In conclusion, important candidate unigenes and miRNA regulators that contributed to the successful invasion of hexaploid cytotypes of S. canadensis have been investigated in the current work, and we have also further inferred ploidy‐related regulation of DNA methylation as an additional modulatory event that occurs to modulate transcriptome reprogramming to drive invasion success. Furthermore, a model for depicting the events involved in ploidy alteration in S. canadensis is summarized in Figure [174]5. Collectively, this work not only describes which molecular processes and functional pathways are likely vital in the successful invasion of polyploids but also offers a valuable dataset for future functional experiments aiming to determine which of these candidate unigenes and miRNA regulators truly underlie the differences in invasiveness between diploid and hexaploid cytotypes. Figure 5. Figure 5 [175]Open in a new tab Graphical summary of molecular responses to ploidy alteration in S. canadensis. Dotted lines and dashed boxes represent the putative regulations CONFLICT OF INTEREST None declared. AUTHOR CONTRIBUTIONS C.C.X. collected plants, performed greenhouse experiment, analyzed the collected data, and wrote this manuscript; Y.M.G. participated in data analysis; J.B.W. conceived the research and contributed to data interpretation and revisions of this manuscript. Supporting information [176]Click here for additional data file.^ (17.5KB, docx) [177]Click here for additional data file.^ (3.7MB, xlsx) [178]Click here for additional data file.^ (29.5KB, xlsx) [179]Click here for additional data file.^ (33.1KB, xlsx) ACKNOWLEDGMENTS