Transcriptional analyses reveal the molecular mechanism governing shade tolerance in the invasive plant Solidago canadensis

Abstract Solidago canadensis is an invasive plant that is capable of adapting to variable light conditions. To elucidate the shade tolerance mechanism in S. canadensis at the molecular level, transcriptome analyses were performed for leaves growing under natural light and three shade level conditions. Many differentially expressed genes (DEGs) were found in the comparative analysis, including those involved in photosynthesis, antioxidant, and secondary metabolism of phenol‐ and flavonoid‐related pathways. Most genes encoding proteins involved in photosynthesis, such as photosystem I reaction center subunit (Psa), photosystem II core complex protein (Psb), and light‐harvesting chlorophyll protein (Lhca and Lhcb), and reactive oxygen species (ROS) scavenging‐related enzymes, such as superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT), were upregulated with the shade levels. Furthermore, most of the DEGs related to secondary metabolite synthesis were also upregulated in the shade conditions. Our study indicates that S. canadensis can respond to shade stress by modulating the expression of several photosynthesis‐related, free radical scavenging‐related, and secondary metabolism‐related genes; thus, this species has the ability to adapt to different light conditions.

to respond to altered environmental conditions by changing its phenotype, a phenomenon that has been recorded for many species that assume different response traits in variable environments (Ortega-Mayagoitia, Hernandez-Martinez, & Ciros-Perez, 2018). This characteristic plays a central role in biological invasions by permitting individuals to colonize different environments and establish viable populations (Du, Liu, Yan, Li, & Li, 2017).
Light is one of the most prominent environmental factors for plant growth (Wang, Chen, Botella, & Guo, 2019). Light not only affects the rate of photosynthesis, antioxidant enzyme activity, and the production of biomass, but also it also regulates all aspects of development in plants (Li, Xin, et al., 2017;Wu, Chen, Sun, Deng, & Chen, 2016). Shade can cause certain metabolic changes as well as imbalances in photosynthesis and carbohydrate production, potentially limiting plant growth and stability (Huang, Zhang, & Liu, 2018;Sajad et al., 2018). Plants are able to respond to changes in light decreased conditions by adjusting a range of physiological and morphological characteristics, such as the specific leaf area, leaf size, biomass allocation pattern, and chlorophyll content (Feng & van Kleunen, 2014;Legner, Fleck, & Leuschner, 2014;Liu et al., 2016).
Plants can respond to light decreased environments by altering gene expression levels, for example, photosystem-related genes are differentially expressed, and most transcription factor (TF) genes are induced under shade conditions in Camellia sinensis (Wu et al., 2016).
Gibberellin biosynthesis genes are found to play a role in shade tolerance in Lolium perenne (Li, Katin-Grazzini, et al., 2017). The expression of genes related to photosystem I and photosystem II could be affected by various light conditions (Rogowski et al., 2019). Shade avoidance and tolerance responses in plants could be reflected in the alteration of expressed genes, such as those involved in hormone signaling and pigment biosynthesis (Ranade, Delhomme, & García-Gil, 2019).
Invasive plants can spread across diverse habitats via individual plasticity (Williams, Mack, & Black, 1995). When dispersed in a shaded environment with low light intensities, invasive plants frequently develop a greater leaf area and biomass to enhance their light capture and utilization efficiency to adapt to the reduced light conditions (Du et al., 2017). This phenomenon reveals that plants have developed some regulatory processes to adapt to shade conditions and may enable their invasion into wider areas.
Solidago canadensis L. (Asteraceae), native to North America, has successfully invaded many areas worldwide. The success of invasion can be related to the plant's high phenotypic plasticity (Du et al., 2017). S. canadensis has a relatively wide niche with different soil properties in Europe (Szymura & Szymura, 2013) and can rapidly spread in non-native regions with different climatic conditions (Li, Du, Guan, Yu, & van Kleunen, 2016). Currently, S. canadensis has been naturalized in many areas and is considered one of the most destructive and widespread invasive species (Wang, Jiang, Zhou, & Wu, 2018), distributed in multiple areas, including farmland, the sides of roads, and upland forests.
To date, the biological processes of S. canadensis in response to shading stress have been unclear. The present study provides a comprehensive description of the differential gene expression patterns in this plant under different degrees of shading. We examine the consequences of three different light intensity levels of shade treatment on the morphological, physiological, and transcriptional characteristics of S. canadensis. Based on the comprehensive analysis, candidate shade stress-responsive genes are identified in comparison with control natural light conditions. These findings will contribute to understanding how S. canadensis adapts to the shade environment and may help to elucidate how it has become an invasive species.

| Plant materials and treatment
In this study, we transplanted 36 rhizomes, respectively, origi- growth status and randomly divided them into four groups. One of these groups was used as control, and the other three groups were covered with one layer, two layers, and three layers of shade net, respectively. In the control treatment (L), plants were cultivated under natural light. In the three different levels of shade treatment with one-, two-and three-layer shade nets, the natural irradiance availability was reduced to 50% (L 1 ), 25% (L 2 ), and 10% (L 3 ), respectively, via black shade nets that were clamped to a metal frame at a height of 150 cm (Semchenko, Lepik, Gotzenberger, & Zobel, 2012). Light intensities were measured by a handheld luxmeter (SMART Digital Lux Meter AS 813). After cultivation for 21 days, we collected leaves from the four treatment groups, with each containing three replicates. The leaves in each plant were collected from the top third to fifth leaves, and these leaves were prepared for RNA extraction, quantitative real-time PCR analysis, and chlorophyll content measurement.

| Measurement of chlorophyll content
One hundred milligrams of each fresh leaf sample from the four treatment was ground with a sterile pestle and morter. The samples were transferred into 15-ml tubes followed by the addition of 10 ml 96% ethanol. After mixing well, the tubes were placed at 4°C for two days. The samples were centrifuged at 17925 g at 4°C for 1 min, and then, the OD value of the supernatant was determined at wavelengths of 470, 665, and 649 nm with an MAPADA UV-1200 spectrophotometer (Shanghai Meipuda Instrument Co. Ltd.). The contents of chlorophyll a, b, and total carotenoid for all samples were calculated according to previous studies as follows: C a = 13.95 × A665-6.88 × A645; C b = 24.96 × A649-7.32 × A665; (Lichtenthaler & Wellburn, 1985).

| mRNA isolation, library construction, and RNA-seq
The leaves of 12 samples from the four treatments were collected for mRNA isolation. After being fully ground, total RNA from frozen leaves was isolated using TRIzol reagent. RNA integrity and purity were checked using an Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit) with a minimum integrity value of 7.5. The cDNA library of each sample was constructed as follows. First, mRNA was isolated from total RNA with magnetic beads containing Oligo (dT) primer.
Subsequently, the enriched mRNAs were fragmented in fragmentation buffer, and then, first-strand and second-strand cDNA were synthesized and amplified using an ABI Step One Plus Real-Time PCR System. Finally, these cDNA library preparations were sequenced using an Illumina HiSeq X Ten from both the 5′ and 3′ ends, and 150bp paired-end reads were generated.

| De novo assembly and functional annotation
Raw reads containing low-quality, adaptor-polluted, and a high content of unknown bases were removed and satisfied with the number of Q20 above 90% before downstream analyses. After filtering out low-quality and ambiguous nucleotides, de novo assemblies were prepared from the clean reads using Trinity (v2.0.6; Grabherr et al., 2011). A transcriptome assembly, generated from all mixed samples, was used as a reference for further gene expression analysis.

| Differentially expressed gene analysis
All clean reads from each sample were mapped to the assembled unigene dataset with Bowtie2 software (Langmead & Salzberg, 2012), and the gene expression levels were calculated with RSEM software (Li & Dewey, 2011). To provide a relatively accurate assessment of transcript abundance, the fragments per kilobase of exon model per million mapped reads (FPKM) value were used as a measure to normalize the unigene expression levels (Mortazavi, Williams, Mccue, Schaeffer, & Wold, 2008). The differentially expressed genes (DEGs) between each group were identified with DEseq2 software with a |log 2 (fold change)| ≥ 1 and adjusted p-value < .05 as the threshold for significantly different expression (Love, Huber, & Anders, 2014).
GO enrichment analysis of all DEGs was performed using hyper, a function of R based on the hypergeometric distribution (Young, Wakefield, Smyth, & Oshlack, 2010). Furthermore, REVIGO was used to visualize the enriched GO term in multiple ways to facilitate interpretation (Supek, Bosnjak, Skunca, & Smuc, 2011). For KEGG pathway enrichment analysis, all DEGs were mapped to pathways in the KEGG database and identified significantly enriched pathways based on a corrected p-value < .05.

| Gene expression network construction
Weighted gene correlation network analysis (WGCNA) is used to construct the expressed gene network and identify stress-related candidate key genes. This method has been widely applied in model organisms and nonmodel organisms for gene expression studies (Fu et al., 2014;Lou et al., 2017). The gene network was constructed using the R package WGCNA (Langfelder & Horvath, 2008), and all the DEGs between the shade treatment and control were applied to WGCNA. These collected genes were clustered to form different modules, which were assembled by unique color. A significant correlation for each module was satisfied with coefficients >0.80 and p-values < .01. Modules that were significantly correlated with shade stress were extracted and exported into Cytoscape (ver. 2.8, http://www.cytos cape.org/) for network visualization (Smoot, Ono, Ruscheinski, Wang, & Ideker, 2011). Hub genes were identified based on the connection degree of each node in the network.

| Statistical analysis
In order to investigate whether the chlorophyll and carotenoid content of plants were changed under different light conditions, the experiment of chlorophyll content measurement for each treatment was repeated three times in this study. The data of chlorophyll and carotenoid content were analyzed by one-way analysis of variance, and the mean (±SD) was the average calculated from three biological replicates. Letters on the bars represent a significant difference at p-value ≤ .05 (LSD test) among different treatments compared with the control group.

| Validation of gene expression by qRT-PCR
The quantitative real-time polymerase chain reaction (qRT-PCR) technique was used to validate gene expression from RNA-Seq. RNA was isolated with TRIzol reagent from each independent biological sample. In order to prevent the presence of DNA in extracted RNA samples, all the RNA samples were treated with DNase I to digest DNA. cDNAs were obtained by reverse transcription of each RNA sample. Specific primers for 12 randomly selected genes were designed using Primer 5 software (Lalitha, 2000). The SYBR Green Master Mix (Applied Biosystems) was used by quantitative PCR. The relative expression of selected genes was normalized by the 2 −∆∆Ct method based on the 18S ribosomal gene as a reference.

| Changes in morphological and chlorophyll content under shade stress
To shading net (L 3 ) were smaller than those of the L 2 but higher than those of the one-layer shading net (L 1 ) treatment groups. In the three different shading treatment groups, the leaf angle and total content of chlorophyll increased with the reduction of light intensity and were significantly higher than the control group. In addition, there was a significant difference in total chlorophyll content between the L with L 1 , L 2 , and L 3 groups ( Figure 2c). The content of chlorophyll a and chlorophyll b also increased with the reduction of light intensity, and most groups were also significantly different from each other ( Figure 2a,b). The carotenoid content in the three treatment groups increased compared with the L group, while there were no significant differences among these three groups ( Figure 2d). In contrast, the ratio of Chl a/Chl b and carotenoid/total Chl in each group decreased with reduced light intensity (Table S1). These results suggest that S. canadensis can exhibit different phenotypes and pigment levels under shade stress compared with natural growth condition.

| Overview of transcriptome profiles
To understand the molecular basis of the phenotypic differences in shade conditions, we generated transcriptome profiles based on RNA-Seq analyses. The leaf samples of three different shade conditions and one natural light condition were used for RNA-Seq analysis with three biological replicates. In total, 12 libraries were constructed and sequenced on an Illumina HiSeq X Ten platform, with each library satisfied with a Q20 percentage greater than 98% (ranging from 98.41% to 98.82%) and single read length of 150 bp ( FPKMs in the L, L 1 , L 2 , and L 3 group, respectively. By redundancy analysis of annotated genes from these four groups, a total of 83,575 unigenes were identified, with 67,237 unigenes detected in all four groups and 1,204 unigenes that were unique to the L group, 592 to the L 1 group, 651 to the L 2 group, and 2,488 to the L 3 group ( Figure   S1). The expressed genes were divided into four ranges according to the FPKM, and approximately 21.3% of the genes were expressed in the range of 0-1, 41.3% in the range of 1-5, 30.7% in the range of 5-30, and 6.7% were above 30.

| Comparison of DEGs under shade conditions
When we examined the number of DEGs for each of these three

| GO functional enrichment analysis
To further explore the function of DEGs, GO functional enrichment was conducted. A total of 3,466 DEGs from the three compared groups were evaluated with GO enrichment and enriched into 49 GO terms classified into three main GO categories ( Figure S2). Many genes were significantly enriched in GO terms related to photosynthesis, reactive oxygen species (ROS), and secondary metabolism synthesis (Table S3). For example, in the category of cellular component, GO terms, such as "photosystem," "photosynthetic membrane," "thylakoid part," "photosystem I," "photosystem II" and "chloroplast," which related to the photosystem were significantly enriched. In the  category of molecular function, GO terms, such as "germacrene-A synthase activity," "pigment binding," "peroxidase activity" and "sesquiterpene synthase activity," which related to ROS response and secondary metabolite synthesis were significantly enriched. In the category of biological processes, GO terms, such as "photosynthesis," "photosynthesis, light harvesting," "protein-chromophore linkage photosynthesis," "light reaction oxidation-reduction process," "light harvesting in photosystem I" and "hydrogen peroxide metabolic process," which related to photosynthesis and ROS response were enriched. All the significantly enriched GO terms were analyzed by REVIGO to further identify representative nonredundant GO terms ( Figure 5). The expression of most of these significant enriched genes was induced under the shade conditions, for example, many photosynthesis-related genes were upregulated in GO terms such as "photosystem I," "photosynthetic membrane," "chloroplast thylakoid photosystem II" and "photosynthesis, light harvesting" (Figure 6). These results suggested that genes involved in these GO terms might play important roles in the response to shade stress.
During this process, S. canadensis might respond to shading stress by regulating the expression of these functional genes.

| Pathway enrichment analysis
According to the KEGG pathway analysis, 81, 117, and 125 path- propanoid biosynthesis," "carbon metabolism," "carbon fixation in photosynthetic organisms," "starch and sucrose metabolism," "stilbenoid, diarylheptanoid and gingerol biosynthesis," "isoflavonoid biosynthesis," "flavonoid biosynthesis" and "plant hormone signal transduction" were significantly enriched (Table S4). The KEGG enrichment analysis indicated that the number of significant enrichment pathways increased as the light intensity was reduced, and the DEGs that were significantly enriched in these pathways are likely to play essential roles in coping with shading stress in S. canadensis.

| DEGs related to photosynthesis
Many DEGs related to the photosynthesis process were collected.
For example, photosystem I reaction center subunit-related genes were also differentially expressed in the shade condition groups when compared with the control (Figure 9a). The number of these photosynthesis-related genes increased with reduced light intensity.
These results suggested that S. canadensis may cope with shading stress by regulating the expression of photosynthesis-related genes.

| Antioxidant enzyme-related DEGs after shading
In this study, some antioxidant enzyme-related genes, such as SOD, CAT, and POD, were differentially expressed in the shade conditions compared with the control, and most of them were upregulated

| DEGs enriched in secondary metabolismrelated pathways
According to the KEGG pathway enrichment analysis, many DEGs from the three compared groups were enriched in secondary metabolism-related pathways, such as "phenylpropanoid biosynthesis" and "flavonoid biosynthesis." The number of DEGs in these two pathways increased with the reduction of light intensity and was significantly enriched in L 2 versus L and L 3 versus L. In contrast, the proportion of downregulated genes among these three comparisons increased in the "phenylpropanoid biosynthesis" pathway, in which 34 DEGs were collected and most were upregulated in the three shade treatment groups compared with the L group, such as the PAL, C4H, 4Cl, and C3′H genes (Figure 9c). In the "flavonoid biosynthesis" pathway, 23 DEGs were collected, and most of these DEGs were also upregulated in the three shade treatment groups compared with the L group, such as the CHS, DFR, F3′H, and FLS genes (Figure 9d).
In addition, some of these secondary metabolite synthesis-related DEGs in this pathway in the L 3 versus L group were downregulated compared with L 2 versus L. These DEGs will respond to changes via induction or suppression in response to different levels of shading.
These results indicated that shading stress might affect the secondary metabolism level in S. canadensis.

| Gene expression network analysis of shade stress-related DEGs
In this study, 3,466 DEGs from three comparison groups were analyzed by WGCNA. These DEGs were parsed into 23 gene expression modules, with the number of genes in each module ranging from 43 to 871 (Figure 10a). The module-trait correlation analysis showed five significant coregulatory modules (correlation coefficient > 0.8, p < .01), with three in the L group shown in turquoise, dark red, and red and two in the L 3 group shown in dark green and blue (Figure 10b).  (Table 2). These hub genes may play a key role in the response of S. canadensis to shade stress.

| Gene expression validation by qRT-PCR
The gene expression results obtained from the high-throughput sequencing were validated by quantitative reverse transcription-PCR (qRT-PCR) analysis. A total of 12 randomly selected key candidate genes from each of four different treatment groups were used for validation. As anticipated, the qRT-PCR results for the different groups ( Figure 11) were consistent with the sequencing data. The primer sequences used in this study are listed in Table S5.

| D ISCUSS I ON
Previous studies have elaborated that invasive species consistently show adaptive potential and rapidly adapt to novel environments (Bertelsmeier & Keller, 2018;Liu, Li, & Xiong, 2015). Furthermore, these species often have the ability to express different phenotypes in response to changing or novel environmental conditions and display higher plasticity, even in low-resource systems (Funk, 2008).
S. canadensis, an invasive plant, could adapt to low-light conditions in a shade stress environment by adjusting the morphology, growth, photosynthetic physiology, and other characteristics, and it exhibits a certain plasticity in response to low light intensity (Du et al., 2017).

| Alteration of phenotype may facilitate the growth of S. canadensis under shade condition
Leaf area and leaf angle insertion are considered to be responses to the shading environment and increase with shade application (Ajmi et al., 2018). Under limiting irradiance, leaves developed large areas to capture more light sources and help the plants maintain better performance (Liu et al., 2016;Yang, Sun, Zhang, Cochard, & Cao, 2014). Furthermore, leaf insertion angles are important in determining the transmission of radiation for vegetation canopies, and the direction of leaves depends greatly on the light exposure (Raabe, Pisek, Sonnentag, & Annuk, 2015). Under shade stress, to absorb more light, the plant leaf angle insertion can be extended to 100°, while under normalized light, to avoid excess sunlight, the leaf angle insertion always ranges from 30° to 60° (Larbi et al., 2015). In our study, all the plants under shading treatment had a larger leaf area and leaf angle insertion than those cultivated under natural light.
Such a strategy will assist the capture of light energy by S. canadensis and improved acclimation to the shading environment.

| Alteration of physiology may contribute to the shade tolerance of S. canadensis
Chlorophyll is the most important photosynthetic pigment in photosynthetic processes, and the chlorophyll content is always significantly affected by shade stress (Alridiwirsah, Harahap, Akoeb, & Hanum, 2018). Under shade stress, plants increase their light use efficiency by increasing leaf chlorophyll a and chlorophyll b, total chlorophyll content, and decreasing chlorophyll a/b ratios (Dai et al., 2009;Gregoriou, Pontikis, & Vemmos, 2007), and if the available light source is reduced in a shade environment for photosynthesis, the chlorophyll content will be increased (Muhidin et al., 2018). Our results revealed a significantly increased chlorophyll content (chl a, chl b, and chl a + b) in the different treatment groups. These results were consistent with the linear increase in shade stress, concurrent with a significant decrease in the chl a/b ratio. The carotenoid content increased in the three treatment groups with no significant differences among them. The ratio of carotenoid content with total chlorophyll was decreased with the reduction of light radiation. This phenomenon indicates that the photosynthesis system in S. canadensis has a wide range of strategies to adapt to shading stress and increase its light use efficiency by regulating the content of chlorophyll.

| Photosynthesis-related genes in S. canadensis play a vital role in the shade environment
To maximize light use for photosynthesis, plants possess fully developed strategies based on short-and long-term responses to adapt to low light intensities (Mathur, Jain, & Jajoo, 2018). Photosynthesis is the key function for plant growth and is highly sensitive to light conditions. In photosynthetic systems, the pigment protein complexes photosystem I and photosystem II are interconnected through the cytochrome (Cyt) b 6 f complex and electron transfer chain (Mamedov, Govindjee, Nadtochenko, & Semenov, 2015).  (Silveira & Carvalho, 2016). The results for the three comparison groups revealed significant differences in photosynthesis-related GO terms and KEGG pathways. In the present study, most of the photosynthesis-related DEGs were enriched in the "photosynthesis," "photosynthesis-antenna proteins," and "carbon fixation in photosynthetic organisms" pathways. Moreover, the "photosynthesis"  Lhca2, and Lhca3, contains LHC trimers that are also essential for coping with different light conditions (Puthiyaveetil, Ibrahim, & Allen, 2012

| Expression of ROS response enzyme-related genes maintains growth health under shade stress
Reactive oxygen species always act as important molecules and play pivotal roles in plant for abiotic stress (Chan, Yokawa, Kim, & Song, 2016;Li et al., 2018;Ma et al., 2011;Qi et al., 2018). Furthermore, ROS were also characterized as important signaling molecules that mediate plant hormone signaling transduction process, such as mitogen-activated protein kinase (Wang, Du, Li, Ren, & Song, 2010;Wang, Du, Zhao, Miao, & Song, 2013) and hormone signaling (Chen, Li, Li, Song, & Miao, 2017;Dong, Wei, & Chen, 2013;Song, Miao, & Song, 2014). The activities of antioxidant enzymes, such as superoxide dismutase (SOD), catalase (CAT), and peroxidase (POD), were significantly changed under different light conditions in Camellia sinensis (Wu et al., 2016). When plants are exposed to high irradiance conditions, the levels of antioxidant enzymes such as POD, SOD, and CAT are significantly altered (Shao et al., 2014). In this study, with increasing shade stress, a number of antioxidant enzyme-related DEGs also increased, such as POD, SOD, CAT, APX, and GPx.
Therefore, under shade stress, antioxidant enzyme activity may be altered in S. canadensis by regulating the expression of related genes to scavenge ROS, thereby maintaining normal plant growth.

| TF gene regulation may promote acclimation to shade stress in S. canadensis
To adapt to adverse environmental conditions, many plants have developed different mechanisms in response to abiotic stresses.

Most transcription factors (TFs) can act as regulatory proteins that
play an important role in these mechanisms (Wang, Yang, et al., 2018). For example, the TIFY family is a novel plant-specific gene family that participates in the regulation of various plant-specific biologic processes (Sun et al., 2017). The AP2/EREBP superfamily is one of the largest and most conserved gene families in plants, and expansion of the AP2/EREBP superfamily may have enhanced the wide adaptability of cotton (Liu & Zhang, 2017). Furthermore, the WRKY and MYB TF families play a vital role in diverse regulatory and multiple stress response processes in plants (Li et al., 2019;Wang, Ren, et al., 2019;Wang, Yue, et al., 2019). In this study, AP2/EREBP and TIFY genes predominated in the three comparison groups of differentially expressed TF-coding genes followed by others, such as MYB-, NAC-, WRKY-, and bHLH TF-coding genes, which were also differentially expressed and mostly upregulated in each group. Therefore, the AP2/EREBP and TIFY TF genes may play an important role in the process of shading stress in S. canadensis, and they may promote an increase in tolerance to stressful environments.

| CON CLUS IONS
Invasive species consistently display the trait of phenotypic plasticity to adapt to various environments. The results from this study confirm that shade levels strongly affect phenotypic changes and physiological responses in S. canadensis. The expression levels of photosynthesis and antioxidant enzyme-related genes were altered in response to shade stress. This phenomenon suggests that

F I G U R E 11
The result of qRT-PCR for randomly selected genes. The blue, cyan, rose, and yellow bar represent the expression of selected genes in L, L 1 , L 2, and L 3 group respectively S. canadensis has enhanced photosynthesis and ROS scavenging abilities to absorb large amounts of light and ensure normal growth.
While the expression of secondary metabolism-related genes was also affected under shade stress, this phenomenon was indicative of an adaptive strategy to retain fitness under harsh conditions. The expression network analyses revealed the important roles of candidate genes for shading tolerance in S. canadensis. This shading tolerance might facilitate S. canadensis to adapt to various environments and successfully develop into an invasive species.

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

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The mRNA-seq data as fastq files from the experiment of S. canadensis were deposited in NCBI SRA database under the accession number PRJNA577291.