Transcriptome analysis reveals plant response to colchicine treatment during on chromosome doubling

Colchicine was commonly used to artificially double chromosomes while the transcriptome changes in colchicine treated plants has rarely been characterized. To understand the molecular mechanism of colchicine on chromosome doubling, we characterized transcriptome data of diploid orchardgrass root after colchicine treatment. Our results showed that 3381 of differentially expressed genes (DEGs) were mainly affected by water stress, 1258 DEGs that were expressed significantly in sample DacR5tr but not in DacR5ck were considered to be mainly affected by colchicine and combination of water and colchicine. These DEGs mainly regulated by colchicine were enriched to gene ontology (GO) accessions of cation binding, catalytic activity, membrane and transporter activity, and enriched to Kyoto Encyclopedia of Genes and Genome (KEGG) pathways of phenylpropanoid biosynthesis, phenylalanine metabolism, plant hormone signal transduction and starch and sucrose metabolism. Genes related to microtubule, spindle, chromosomal kinetochore, vesicle, cellulose and processes of cytoplasm movement, chromatid segregation, membrane and cell wall development were inhibited by colchicine. Our results revealed that colchicine restrained the microtubules and inhibited gene expression of cytokinesis, which might slow down the cell activity, delay the cell into anaerobic respiration, resulting in apoptosis at late stage, and relieving of waterlogging.

of alpha-and beta-tubulin dimers, and then pull chromosomes toward two polar ends of the cell. Colchicine can associate with tubulin dimers to reduce the amount of microtubules in the synthesis process 21,23 , which may prevent spindle formation, shorten the length of spindle fibers, even result in non-division of cells. Although it has widely been accepted that colchicine blocks the cell division by inhibiting the synthesis of microtubules, the exact mechanism on chromosome doubling remains unclear. Furthermore, the mechanism of colchicine was mostly studied on human health, rarely on plants 24,25 .
Dactylis L., commonly known as orchardgrass, belongs to Festucaea in the Poaceae family, is an important perennial, cool-season and tall-growing grass, and consists of diploids, tetraploids and a hexaploid. Due to being widely used as forage and hay, Dactylis has been spread almost over every continent. Orchardgrass is one of the top four perennial forage grass, with about 14 million pounds of seed harvested each year 26 . To date, Dactylis species commonly used as forage were tetraploids, such as D. glomerata subsp. glomerate, subsp. hispanica and subsp. marina 27 . Although there are many diploid species in orchardgrass, many of them are threatened by climate warming and habitat degradation. To supply and expand the gene pool for Dactylis breeding programs, diploids have been used for breeding through artificial chromosome doubling of diploids [28][29][30] .
To double the chromosome of plants using colchicine, the common process is to soak plant tissue in colchicine solution 16,31 . Water, colchicine, and potential synergism of water and colchicine could cause transcriptome changes in treated individuals, which has rarely been characterized. To explore the mechanism of colchicine on chromosome doubling, and reveal the effect of colchicine on the plant as one stress, diploid orchardgrass, D. glomerata subsp. smithii L., was treated with colchicine. Transcriptome analysis was performed to reveal the change of gene expression in response to the treatment.

Results
Transcriptome assembly and gene annotation. Total of 48999576, 50310774 and 51255434 clean reads were obtained from 51481126, 52867138 and 53823216 raw reads in three Dactylis samples of DacR0 (0 h, nontreated), DacR5ck (5 h, water treated) and DacR5tr (5 h, colchicine treated). Full-length of transcriptomes was reconstructed using the clean reads via Trinity because of the lack of a reference genome. Total of 251227 transcripts were finally assembled, of which, 172524 unigenes were selected since those were the longest transcripts of each gene. The length of unigenes ranged from 201 to 15543 base pairs (bp), with mean length of 583 bp, median length of 331 bp. The length of N50 and N90 were 843 bp and 249 bp, respectively.
BLAST search against the NR and SwissProt databases found that, out of the 172524 unigenes, 67653 genes were predicted as CDSs with an average length of 603 bp, ranging from 30 bp to 15102 bp. In addition, ESTScan showed that 92127 unigenes that were not hit in BLAST were predicted as CDSs with an average length of 305 bp, ranging from 51 bp to 6609 bp.
After genes were functionally annotated, the genes were described and estimated according to the biological function via GO project. Biological processes, cellular components and molecular functions were applied to screen the genes' attributes. Among the unigenes, 50502 genes were attached to at least one GO term. In the biological process term, genes were mainly attached to cellular process (52.82%), metabolic process (51.49%) and single-organism process (38.82%). The cell, cell part and organelle, were found as the first three elements in cellular component with a proportion of 28.67%, 28.65% and 19.34%, respectively. In the molecular function term, the first three major elements were binding (52.66%), catalytic activity (42.40%) and transporter activity (5.86%).
Total of 21567 unigenes were annotated in KEGG database, and were classified into 131 secondary pathways of 5 primaries. Scanning in the secondary pathways, translation of genetic information processing term held most of the unigenes (2989 genes, 13.86%), followed by carbohydrate metabolism of metabolism term (2323 genes, 10.77%).

DEGs analysis in DacR5tr vs DacR0 and DacR5ck vs DacR0.
The normalized read count of each unigene was applied to differential expression analysis. After filtering, 3381 and 3582 differentially expressed genes (DEGs) were screened out in comparison of DacR5ck vs DacR0, and DacR5tr vs DacR0, respectively (Fig. 1a,b). Compared with untreated DacR0, 1352 and 1357 DEGs were up-regulated in DacR5ck and DacR5tr, whereas 2029 and 2225 DEGs were down-regulated in DacR5ck and DacR5tr, respectively. Total of 1978 DEGs were common between DacR5ck vs DacR0 and DacR5tr vs DacR0 (Fig. 1c), in which 1361 down-regulated DEGs and 561 up-regulated DEGs were detected in the both comparisons, while 56 DEGs were inconsistent between the two comparisons.
KEGG pathway enrichment analysis showed the highest number of DEGs was involved in phenylpropanoid biosynthesis, followed by plant hormone signal transduction and phenylalanine metabolism in DacR5ck vs DacR0. However, the largest part of KEGG pathway in DacR5tr vs DacR0 was ribosome, followed by the phenylpropanoid biosynthesis, starch and sucrose metabolism. In addition, these four pathways were enriched as down-regulation.
The fold change of DEGs showed that two alcohol dihydrogen genes were most upregulated in DacR5ck vs DacR0, followed by the genes related to ATP metabolic process and oxidation-reduction process. The genes with the highest up-regulation in DacR5tr vs DacR0 were related to proteolysis and oxidation-reduction process.
GO and KEGG enrichment analysis were processed with these DEGs. Total of 2225 DEGs were involved in 2990 GO accessions, among which, DEGs were mainly enriched in biological process and molecular function terms. Among the DEGs, up-regulated DEGs were enriched mainly in cation binding and metal ion binding terms, while down-regulated DEGs were enriched mainly in catalytic activity, membrane and oxidation-reduction process (Table 1). KEGG enrichment analysis showed that DEGs were involved in pathways of plant hormone signal transduction, phenylpropanoid biosynthesis, starch and sucrose metabolism, apoptosis, chemical carcinogenesis, etc. Among the 20 KEGG pathways, 10 were mainly down-regulated and 10 were up-regulated DEGs ( Table 2). DEGs related to cell division. DEGs were deeply scanned to reveal the expression of genes related to cell division by comparing the transcriptome of DacR5tr vs DacR5ck. The results of GO enrichment showed that more than 170 GO accessions identified involved in cell division processes, which were descripted as organization or processes of microtubule, molecular motor, spindle, kinetochore, cell cycle checkpoint, Golgi, membrane, vesicle, chromosome segregation, cell wall, cell cycle, etc (Part of the GO terms was showed in Table 3).
Total of 21 GO accessions related to microtubule were revealed, including processes of microtubule binding, motor activity, microtubule anchoring, microtubule polymerization and depolymerization, regulation of polymerization and depolymerization, and organization of microtubule cytoskeleton and organizing center. Total of 16 DEGs were descripted as microtubule-related genes, in which 7 DEGs were up-regulated, 9 DEGs were down-regulated. For almost every GO accession of microtubule, the number of up-regulated DEGs was less than that of down-regulated DEGs. Four GO accessions related to motor activity were scanned out, one up-regulated and three down-regulated DEGs. Total of 6 unigenes were spindle-related DEGs, playing roles in 7 GO accessions including spindle assembly, organization, attachment of spindle microtubules to kinetochore, etc. Three up-regulated and 5 down-regulated kinetochore-related DEGs were uncovered, and involved in 10 GO accessions. Two down-regulated cell cycle checkpoint DEGs were involved in one GO accession. Total of 14 Golgi-related GO accessions were revealed, 4 up-regulated and 10 down-regulated DEGs, which were described as Golgi membrane, Golgi transport, etc. The transport and organization of vesicle contained 19 up-regulated and 33 down-regulated DEGs, which were involved in 20 GO accessions. Eight up-regulated and 10 down-regulated segregation-related DEGs distributed in 17 GO accessions, which were involved in nuclear division, sister For DEGs involved in the GO accessions mentioned, the fold changes of the 16 microtubule-related genes ranged from −5.2 (down-regulated) to 8.9 (up-regulated), most of them were up-or down-expressed at a level from 1 to 2. However, these 16 DEGs were expressed 3.4 fold low to 3.4 fold high between DacR5tr and DacR0. In the microtubule-related genes, 4 of 16 DEGs were found to regulate the motor activity according to GO accessions, which were related to movement in cell, and displayed 5. Calcium signaling pathway ko04020 6 3 9 Flavonoid biosynthesis ko00941 8 0 8 Stilbenoid, diarylheptanoid and gingerol biosynthesis ko00945 5 1 6 Flavone and flavonol biosynthesis ko00944 4 0 4 Degradation of aromatic compounds ko01220 4 0 4
In 5hTR of the accession PI 237607, the expression level of c100325_g1, c101064_g1, c102023_g1, c96405_g1, c101317_g1, c93529_g2 and c49608_g1, were higher, while c71496_g1 and c101667_g3 were lower than that in the 5hck. Nine of the 11 differentially expressed genes in 5htr vs 5hck of PI 237607 were detected with the same expression trend as those in PI 441032. However, the c89024_g1 and c88478_g2 genes showed different expression trends in these two accessions, which were down-regulated in PI 441032 but up-regulated in PI 237607 after treatment (Fig. 2). To reveal the gene expression changes under a long-time treatment, the 11 DEGs were also tested in samples of PI 441032 and PI 237607 under 24 h treatment of water (24 hck) and colchicine (24htr). After 24h treatment, all of the genes were down-regulated in the PI 441032, while most of the genes were down-regulated in the PI 237607 (Fig. 2).

Discussion
Colchicine is one antimitotic agent, and used frequently for chromosome doubling. However, it is also a toxic chemical, possessing great threat to normal plant life. In treatment of Platanus acerifolia using colchicine, higher frequency of chromosome doubling was correlated with higher concentration of colchicine, but accompanied with a lower survival rate. Under colchicine treatment, survival rate of Miscanthus plants significantly depended on each genotype, colchicine concentration and exposure time 35 . High polyploidization rate and high lethality prompted us to select a suitable treatment method 16,[35][36][37] . In our study, 2.5 mM colchicine (0.1%, w/v) and 2% DMSO was used to characterize the changes of genes expression in response to colchicine treatment for 5 h and 24 h.
To double the chromosome of plants using colchicine, scientists often soaked plant roots in colchicine solution 16,31 . Effects on treated plants were from water, colchicine and potential synergism of water and colchicine, which could cause transcriptome changes. In our study, we analyzed the transcriptome of D. glomerata subsp. smithii under 0 h treatment, 5 h water treatment and 5 h colchicine treatment. 3381 DEGs were affected by water. Compared with DacR0, 3582 DEGs were screened out from DacR5tr which was treated by colchicine solution, thus these DEGs were mainly affected by water, colchicine, and potential synergism of water and colchicine. In  comparison of DacR5tr with DacR5ck, 2247 DEGs were identified to be mainly affected by both colchicine and potential synergism of water and colchicine. Comparing these 2247 DEGs with DacR0, 976 were significant differentially expressed in DacR5ck under water stress, and were expressed with high or low level in DacR5tr after colchicine was added (Fig. 1e). These DEGs were affected by water or colchicine. The remaining 1258 DEGs that were expressed significantly in DacR5tr but not in DacR5ck, were considered to be mainly affected by colchicine and potential synergism of water and colchicine. Colchicine is a bioactive alkaloid, one of the antimitotic agents. Its effect on human health and diseases was frequently reported. With suitable dosage, colchicine has been used to inhibit cancer, induce apoptotic cell death, and as an anti-inflammatory for disrupting tubulin. Its effect on the immune system included inhibition of neutrophil chemotaxis, superoxide production, etc 24 . Its well-known effect on plant cells was only its inhibition of microtubules in mitosis 38,39 . In our study, transcriptomes were analyzed to explore its effects on gene expression. According to GO analysis, DEGs were mainly enriched in 16 GO terms with data of DacR5tr vs DacR5ck. DEGs included in these 16 GO terms had lower or higher expression in DacR5tr than in DacR5ck, whose changes were mainly affected by colchicine. Specially, DEGs related to cation binding especially metal ion were induced to a high level of expression, while more DEGs related to the other GO terms were restrained at low level expression under colchicine stress, including the major terms of catalytic activity, membrane, transporter activity, etc. The genes in KEGG pathways of phenylpropanoid biosynthesis, phenylalanine metabolism, and plant hormone signal transduction were mainly down-regulated, while genes in metabolism processes of starch, sucrose, cysteine, methionine, galactose, and chemical carcinogenesis were up-regulated. Phenylpropanoids biosynthesize the amino acids of phenylalanine and tyrosine, and are found throughout the plant kingdom. Phenylpropanoids are not only indispensable for structural polymers, but also involved in responses to biotic or abiotic stresses 40 . Plant hormones play an important role in regulating plant growth and development by the transmission of molecular signals in cells 41,42 . Decrease in expression of genes related to phenylpropanoid, phenylalanine, plant hormone and increase in expression of more chemical carcinogenesis genes caused by colchicine might partly explain why high mortality rate was accompanied with chromosome doubling when colchicine was used.

Gene ID Gene function annotation
Mitosis involves replication of chromosomes and division of the mother cell into two daughter cells leading the individual to grow 22,43 . In the mitotic process of plant cells, organization of microtubule, motor protein, spindle, chromosomal kinetochore and vesicle play important roles, while several processes of movement of cytoplasm, segregation of chromatids, check of metaphase, development of membrane and cell wall around the cell plate in middle of cells, decide the final cell division 44,45 . In our research under colchicine treatment, more than 170 GO terms were found to be related to microtubules, motor proteins, kinetochores, cytokinesis, etc. According to previous studies, assembly of microtubules was inhibited by binding of colchicine with alpha-and beta-tubulin dimers, without reducing disassembly process 46 . Assembly and depolymerization of microtubules were less frequent than its disassembly, resulting in the reduction of microtubules. However, in our research, more genes in GO terms related to microtubules, microtubule binding, tubulin complexes and microtubule associated complexes were down-regulated, while GO terms with more up-regulated genes were mainly related to microtubule anchoring, cytoskeleton and organizing center. Total of 6 spindle-related genes were found, including 3 up-regulated and 3 down-regulated genes. Gene expression for assembly of spindle was increased, while genes for spindle microtubule and its attachment to kinetochore were restrained. Thus, colchicine influenced the microtubule and spindle by not only reducing density of available tubulin dimers but also regulating the related genes' expression. The expression of genes related to segregation of sister chromatids was down-regulated, resulting in inhibition of sister chromatid segregation under colchicine treatment. The metaphase checkpoint guarantees the exact segregation of chromosomes and promotes cell cycle to anaphase 47 . However, the expression of checkpoint genes was down-regulated by colchicine, which possibly leads the cell to proceed to anaphase unsuccessfully. Furthermore, inhibitions caused by colchicine also occurred in movement of cell components and cytokinesis. Relatively, the cell plate and cell wall in the middle of the cell was not easy to be formed with less vesicles and less biogenesis of cellulose under the inhibition of colchicine. As a result, the cell division cannot be completed under colchicine.
In our study, water stress was actually a waterlogging stress. Under waterlogging, the plant encounters the condition of hypoxia (deficiency of O 2 ) or anoxia (absence of O 2 ). Growing in oxygen deficient conditions, plant cell activity was energized by root anaerobic respiration. As a result, ethanol was accumulated. Furthermore, it was reported that anaerobic polypeptides (APNs) including lactate dehydrogenase and alcohol dehydrogenase (ADH) were increased under low oxygen environment as ADH was the major APNs 48 . Total of 26 ADH DEGs from water stressed Dactylis were expressed, 8 down-regulated and 18 up-regulated. However, 16 of the 26 ADH DEGs were down-regulated under colchicine treatment in comparison with DacR5ck. The highest expression level of ADH was 10-fold up-regulated in DacR5ck vs DacR0, while it was only 3-fold up-regulated in DacR5tr vs DacR0, suggesting that colchicine might relieve the stress from waterlogging. Moreover, KEGG pathways showed that more apoptosis related genes were down-regulated in DacR5tr than those in DacR5ck. In this study, treatment was performed with 2.5 mM colchicine (0.1%, w/v) mixed with 2% dimethyl sulfoxide for 5 h. Compared to previous studies 49,50 , concentration of colchicine used in our treatment was low, and treatment time was short, which may be moderate to orchardgrass young plants. In addition, our results suggested that colchicine restrained the microtubules and inhibited the gene expression of cytokinesis, which might slowdown the cell activity, delay the cell into anaerobic respiration, resulting in apoptosis at a late stage and relieving the stress from waterlogging.
In summary, the inhibition of colchicine on microtubules was suggested by gene expression. Numbers of genes related to spindle, chromosomal kinetochore, vesicle, cellulose and processes of cytoplasm movement, chromatid segregation, membrane and cell wall development were also affected by colchicine. Multiple organizations and processes worked together to prevent cell division. Our results suggested that colchicine might have a buffer effect on relieving waterlogging stress. For further research, we suggested that solid medium mixed with colchicine be used for chromosome doubling to eliminate the effect caused by water, such as agarose culture medium, in this case, how genes respond to colchicine treatment during chromosome doubling can be precisely figured out.

Methods
Plant materials and treatment. Seeds of diploid plant of D. glomerata subsp. smithii Link (2n = 2x = 14) were kindly provided by the GRIN (USDA-ARS Germplasm Resources Information Network). Two accessions, PI 237607 from Spain and PI 441032 from United Kingdom, were used. Seeds were germinated under room temperature, and transplantation of seedlings was followed after plant expanding with 3-4 leaves. Seedling plants grew in greenhouse with regular irrigation. Plants were identified via observation of morphology and cytology.
After one month, young plants were used for treatment. To keep the plants away from physical damage or at least minimize the physical damage during collection, whole young plants with soil around the roots were moved softly from the container. Soil was washed off softly with running water. Dozens of clean plants were separated randomly into 2 groups after their roots were cut off about 3 cm from the bottom. One of the groups grew in water, whereas another one was treated by 2.5 mM colchicine (0.1%, w/v) mixed with 2% dimethyl sulfoxide (DMSO). Roots of samples in each group were separated randomly into more than 3 sub-groups, each with a few dozen of plants. Roots tissues were harvested from each sub-group at the time of 0 h, 5 h and 24 h, and preserved immediately in liquid nitrogen for RNA isolation. After treatment, all of the samples were transplanted and grew in greenhouse for seeds setting, identification and generating of chromosome doubled plants as new germplasm resource for Dactylis breeding.
RNA isolation and cDNA preparation. Total RNA from root tissues was isolated using TRI Reagent ® (Sigma-Aldrich, USA) following manufacturer's protocol. To eliminate the potential contamination from genomic DNA, RNAs were treated by RNase-free DNase set (QIAGEN, USA). RNA quality was evaluated using electrophoresis on 1.0% agarose gels, while quantity was measured using NanoDrop 2000 Spectrophotometer (Thermo Scientific, USA). cDNA was synthesized using QuantiTect Reverse Transcription Kit (QIAGEN, USA). Total RNA and cDNA of each accession under treatments were prepared and used for Solexa sequencing and real-time PCR, respectively. Three RNAs samples of accession PI 441032 were sequenced using Solexa sequencing, and were named as DacR0 (0 h, non-treated), DacR5ck (5 h, water treated) and DacR5tr (5 h, colchicine treated). The cDNAs of 5 h water-treated and 5 h colchicine-treated samples of PI 441032 were used to validate differentially expressed genes identified in transcriptomic analysis using qPCR. The 24 h water-treated and 24 h colchicine-treated samples of PI 441032 were used to observe trends in gene expression in a longer colchicine treatment time. The cDNAs of 4 samples of PI 237607 were used to examine if the gene expression trend in this accession was similar to those in the PI 441032 using qPCR.
Sequencing and gene annotation. For each sample, at least 10 μg of total RNA with A260/A280 > 2.0, A260/A230 > 2.0 was sent to the Novogene (Beijing, China) for Solexa sequencing commercially. The mRNA was purified from total RNA using oligo (dT) magnetic beads after quality control. Synthesis of cDNA was processed via steps of fragmentation of mRNA, the first strand was synthesized using random hexamer primer and reverse transcriptase, and the second strand was synthesized using RNase H and DNA polymerase I. Final cDNA library was harvested after terminal repair, size selection, adapter ligation, PCR amplification and purification. Library quality was evaluated on Agilent Bioanalyzer 2100 system (Agilent Technologies, USA). The index-coded samples were clustered using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) on cBot Cluster Generation System. Prepared cDNA library was fed into HiSeq platform and paired-end reads were generated.
Raw reads were obtained from the sequencing-received image data, and then clean reads were filtered by discarding adapter-containing, poly-N-containing and low-quality reads through in-house perl scripts. Q20, Q30, GC-content and sequence duplication were calculated for clean data. Because of the absence of a reference genome, the transcriptome was assembled from clean reads using Trinity 51 with min_kmer_cov set to 2 and the other parameters at default.
Identification and enrichment analysis of DEGS. To quantify the gene expression levels for each sample, clean data was mapped back to the assembled transcriptome, read count for each gene was obtained from the mapping results. All these were performed using RSEM 52 . Subsequently, the read counts were adjusted by edgeR program package through one scaling normalized factor. DEGs were respectively screened out from comparisons of DacR5ck vs DacR0, and DacR5tr vs DacR5ck by DEGseq (2010) R package with parameter setting as q-value < 0.005 and|log2 (fold change)| > 1.
Differentially expressed genes were submitted to GO term for enrichment analysis using the GOseq R packages based on Wallenius non-central hyper-geometric distribution 53 . The DEGs were also applied to KEGG pathway enrichment analysis implemented by KOBAS 54 to test the statistical enrichment of DEGs in KEGG pathways.
Detection for gene expression with qPCR. Total of 19 unigenes were randomly chosen for validating the expression of DEGs presented in the transcriptomes and detecting the variation of their expressions in samples using qPCR. Primers for unigenes were designed by Primer3 version 0.4.0 55 and Beacon designer version 8.14 (PREMIER Biosoft, Palo Alto, CA, USA), suitable primers were screened out for the final qPCR experiment by estimating its specificity, amplification efficiency, etc. The cDNA of each sample prepared was diluted 1:10 with nuclease-free water before qPCR analysis. Each reaction was performed in total of 15 μl mixture containing 7.5 μl of 2 × SYBR ® Green ROX qPCR mastermix (QIAGEN, USA), 1 μl of cDNA template, 0.5 μM of each primer, with nuclease-free water supplied to final volume. The qPCR was carried out in 96-well blocks with ABI PRISM ® 7000 Sequence Detection System (Applied Biosystems). qPCR amplification was performed under a program: 10 min at 95 °C, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Dissociation curve was obtained by heating the amplicon from 60 to 95 °C. Technical and biological triplicates were applied for qPCR reaction. Mean relative expression ratio was calculated using the 2 −ΔC(t) method 56 with one stably expressed gene as reference gene.