Molecular basis underlying the successful invasion of hexaploid cytotypes of Solidago canadensis L.: Insights from integrated gene and miRNA expression profiling

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.


| INTRODUC TI ON
Amid growing evidence of biological invasion impacts on global biodiversity, ecosystem functioning, species conservation, and even social and economic activities (Dyer et al., 2017;Rejmánek, 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., 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, 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, 2013;Prentis & Pavasovic, 2013;Rius & Darling, 2014). For instance, Rius and Darling (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, 2011;Pandit, White, & Pocock, 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, 2011) and have cumulatively been acknowledged as a latent advantageous attribute of plant invaders (Pandit et al., 2011;te Beest et al., 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, 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, 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, 2007;Jackson & Chen, 2010;. 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, 2013;Song & Chen, 2015) and allowing polyploids to occupy new habitats, grow vigorously and improve adaptation in novel environments (Madlung, 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, 2008;Ha et al., 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, 2007;Ha et al., 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.

| Study species
Solidago canadensis (Asteraceae), a perennial weed native to North America (Werner, Bradbury, & Gross, 1980) where it exists in a diploid, tetraploid, or hexaploid cytotype (Melville & Morton, 1982), has invaded a wide geographical range globally, including New Zealand, Australia, Europe, and Asia (Abhilasha, Quintana, Vivanco, & Joshi, 2008;Szymura, Szymura, Wolski, & Swierszcz, 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, 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, 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, 2011). Additionally, the roots, stems, and leaves of hexaploid cytotypes were morphologically and anatomically distinct from their diploids (Wang, 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.

| 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 1, Figure 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 1; Figure 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 S1).

| 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 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 S1).

| 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., 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., 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 S1 for more unigene annotation details.

| Unigene quantification and differentially expressed unigene (DEU) analysis
The Bowtie2 program (version: v2.2.5, Langmead & Salzberg, 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, 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 (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., 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.

| 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 (2010). An absolute value of log 2 fold-change ≥1 and a q-value <0.001 was set as the threshold to identify DEMs.

| Prediction of miRNA targets
To predict the potential genes targeted by miRNAs, the Targetfinder (version: 1.5, Fahlgren & Carrington, 2010) in combination with psRobot (version: 1.2, Wu, Ma, Chen, Wang, & Wang, 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.

| 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, 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.

| 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 A1 in Appendix, Supporting Information Appendix 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.

| 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 2c). For the invasive populations, all individuals were observed to be hexaploid cytotypes with a chromosome number of 2n = 6x = 54 (Table 1; Figure 2d). However, tetraploid cytotypes with a chromosome number of 2n = 4x = 36 or mixed-cytotypes reported by Li (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 A2 in Appendix). These numbers are comparable to those generated in other polyploid studies (e.g., Vigna et al., 2016;Zhou et al., 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 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 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 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) F I G U R E 1 Map of China sampling sites for 43 populations of S. canadensis described in 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), MYBrelated (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 3). These TF genes had differential expression patterns, implying a variety of regulatory modes.

| Functional and pathway analysis of ploidyresponsive 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 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 F I G U R E 2 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 1 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   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 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.

| 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 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 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, 2006), Oryza sativa (Morin et al., 2008), and Citrus trifoliata (Song et al., 2010).

| 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 miR- 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   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.

| Integrative analysis of gene and miRNA expression confirms that environmental adaptationrelated 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 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, 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 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., 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 A5 in Appendix).

| D ISCUSS I ON
A large number of works have investigated ecological and evolutionary elements responsible for successful invasion (Hahn et al., 2012;Thébault et al., 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.
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., 2018). Hence, it seems that the pattern of  (Xie & Zhang, 2015). Ghani et al. (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.

| 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., 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 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 F I G U R E 3 A bar graph representing the differential expression of TF-encoding genes in diploid and hexaploid cytotypes of S. canadensis. Yellow indicates the upregulated genes and blue down-regulated genes in hexaploid cytotypes 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, 2007;Salmon & Ainouche, 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, 2005). Additionally, Madlung et al. (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  and Cymbopogon (Lavania et al., 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   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 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. 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., 2006;Wang, Zhou, Xu, & Gao, 2010). Auxin-Regulated Gene involved in Organ Size (ARGOS), a gene deeply induced by auxin, participate in organ size regulation. Wang, Zhou, et al. (2010)  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., 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., 2011). Echinacea purpurea autotetraploids showed that the induction of polyploids resulted in higher caffeic acid derivatives and alkamides (Xu et al., 2014).

| Potential roles of transcriptional alterations in the successful invasion of hexaploid cytotypes
Evidence of the influence of polyploids on chemical profiles has also been recorded in allopolyploids. Banyai et al. (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 . Based on the functional annotation of DEUs, a total of 120 core enzyme unigenes encoding CYPs were differentially expressed (Table A6 in Appendix). In addition, CYPs are one of the largest superfamilies of enzyme proteins (Darabi, Seddigh, & Abarshahr, 2017 Furthermore, a model for depicting the events involved in ploidy alteration in S. canadensis is summarized in Figure 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.

ACK N OWLED G M ENTS
This work was supported by the National Natural Science Foundation of China (31570539, 31370258).

CO N FLI C T O F I NTE R E S T
None declared.