Low-temperature tolerance and transcriptome analyses during seed germination of Anabasis aphylla

ABSTRACT This study aimed to explore differences in the low-temperature tolerance of a desert plant, Anabasis aphylla, at different seed germination stages during the fall using transcriptome sequencing to identify related genes. The survival rate of the seeds decreased with lower temperatures at different germination stages, and the tolerance of the stages in decreasing order was as follows: imbibition (I) ≈ testa rupture (II) > testa removal (III) > elongation (IV). Superoxide dismutase, peroxidase, and catalase activities and the malondialdehyde content showed a downward trend after stage II. Transcriptome sequencing of seeds, including ripe and dry dominant seeds (W1) and seeds during different germination stages [imbibition (W2), testa removal (W3), and elongation (W4)] was performed. Some genes were downregulated during seed germination, which mainly included cytochrome P450, oleosin, ethylene-responsive transcription factor, and low temperature-induced protein. This study suggests that downregulated of some genes may result in a decrease in low-temperature tolerance.


Introduction
Anabasis aphylla belongs to the Chenopodiaceae family and is the predominant and constructive species in the Junggar Basin of northwestern China. It is highly saline alkalineresistant, cold-resistant, and drought-resistant, and stabilizes the sand and improves the soil (Feng et al. 2011). A. aphylla also plays an important role in the construction of proper plant communities and stabilization of the communities (Zhao et al. 2003). Low temperature is one of the main abiotic stresses affecting the growth, development, and spatial distribution of plants (Sun and Huang 2011) and tends to reduce and delay the germination rate of seeds and even causes germination failure (Wang et al. 2008). However, most desert plants such as Salsola affinis (Wei et al. 2007;Luo et al. 2014), S. korshinskyi (Li et al. 2012), Lepidium perfoliatum , Haloxylon ammodendron (Tian 2014), A. elatior (Han et al. 2011), and A. aphylla (Peng et al. 2018) are highly adaptable to low-temperature environments and can thus germinate under such conditions. Our previous studies showed that the seeds of A. aphylla can germinate in early spring when the ground is moist due to melting snow (Wang et al. 2017;Peng et al. 2018) or precipitation occurs before freezing. In the Junggar Basin, precipitation pulses in the fall and melting snow before freezing increase the surface water retention time in deserts with alluvial fan soil, facilitating seed germination of desert plants. In addition, increased precipitation caused by climate changes and the delay in freezing time also promote seed germination in the fall. However, the minimum temperature of the Junggar Basin in the winter is below −30°C, and the ground is covered with snow for an extended period of time. The seedlings that had germinated in the fall are thus covered by snow during winter, and differences in germination time contribute to variations in the germination stages of the seeds. From the perspective of internal factors, the cold tolerance of different stages of seeds remains unknown, and are there differences in cold tolerance at different stages and what is the effect on gene expression? Hence, evaluation of the overwintering safety of A. aphylla propagules at different germination stages and under low ground surface temperatures is necessary. Moreover, the anti-freezing mechanism of A. aphylla seedlings under climate change conditions will provide a reference for the study of stress resistance in desert plants.
With the rapid development of high-throughput sequencing technology, transcriptome sequencing has been extensively used in studying the differential expression of genes in living beings at different developmental stages. In previous studies, transcriptome sequencing was used to investigate gene expression at different developmental stages of Medicago sativa (Liu et al. 2018) and Tussilago farfara (Jia et al. 2017) to screen and identify related key genes. In this study, low-temperature tests were conducted on A. aphylla seeds at different germination stages to analyze the difference in their lowtemperature tolerance. In addition, transcriptome sequencing of A. aphylla seeds from four germination stages after lowtemperature stress was performed to study differences in gene expression, which provides a reference for research on genes related to low-temperature tolerance.

Experimental materials
Seeds of A. aphylla were collected from a natural population on the southern margin of the Junggar Basin in Xinjiang Uygur Autonomous Region, China, in October 2018 (45°2 2 ′ 43.4 ′′ N, 84°50 ′ 32.5 ′′ E; altitude, 843 m) by selecting plump and fully ripened seeds and storing these under dry conditions. One portion of the A. aphylla seeds at different germination stages was used for low-temperature tolerance and physiological experiments, and the other portion was used for transcriptome sequencing, with three biological replicates for each germination stage.
2.2. Detection of subzero-temperature tolerance at different germination stages Sterilized A. aphylla seeds were germinated at a constant temperature of 25°C to achieve four germination stages: imbibition (I), testa rupture (II), testa removal (III), and elongation (IV). Upon reaching each of the germination stages, the A. aphylla seeds (n = 40 seeds per replicate and three biological replicates for the experiment) were slowly cooled to the target temperatures (i.e. −2°C, −6°C, −10°C, −14°C, −18°C, and −22°C) for a seven-day low-temperature treatment, followed by slowly warming to 25°C. The A. aphylla seeds (at each germination stage) that continued to grow subsequently were counted and considered to have low-temperature tolerance or were undergoing recovery from low-temperature stress. The variation in temperature during slow-cooling and slow-warming was 4°C per 12 h.

Determination of physiological indices at different germination stages
All indicators were measured under fresh weight conditions. Samples germinated to stages I, II, III and IV at a low temperature of 4°C were collected, surface water was aspirated, and 0.2 g of sample at each stage was weighed and then used for the determination of SOD, POD, CAT, and MDA indicators. Three replicates were used. POD activity was determined following the guaiacol method as described elsewhere (Zou 2000). CAT activity was determined according to Gao (2006). SOD activity was determined using a nitroblue tetrazolium assay (Zou 2000). The MDA content was determined by a thiobarbituric acid colorimetric assay (Zhang et al. 2004).
2.4. Total RNA extraction, cDNA library construction, and transcriptome sequencing Transcriptome sequencing of A. aphylla seeds germinated at a low temperature of 4°C, including ripe and dry dominant seeds (W1) and seeds during different germination stages [imbibition (W2), testa removal (W3), and elongation (W4)], was performed. A total amount of 1.5 µg RNA per sample was used as input material for the RNA sample preparations, and three biological duplicates were tested. RNA degradation and contamination were monitored on 1% agarose gels. RNA purity was assessed using a NanoPhotometer® spectrophotometer (IMPLEN, Westlake Village, CA, USA). RNA concentrations were measured using a Qubit® RNA Assay Kit in Qubit® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). RNA integrity was assessed using an RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, San Diego, CA, USA). The qualified total RNA samples were sent to Novogene Co., Ltd. (Beijing, China) for cDNA library construction. Sequencing was performed in an Illumina HiSeq 2500 high-throughput sequencing system.

Functional annotation of genes
Raw data (raw reads) in fastq format were first processed through in-house Perl scripts. In this step, clean data (clean reads) were obtained by removing reads containing adapters or poly-N or were of low quality. At the same time, the Q20, Q30, GC content, and sequence duplication level of the clean data were calculated. All downstream analyses were thus based on high-quality clean reads. Transcriptome assembly was accomplished using Trinity (Grabherr et al. 2011) with min_kmer_cov set to 2 by default and all other parameters set to default values. Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (Protein family); KOG/COG (Clusters of Orthologous Groups of proteins), Swiss-Prot (a manually annotated and reviewed protein sequence database), KO (KEGG Ortholog database), and GO (Gene Ontology).

Quantification of gene expression and differential gene screening
Gene expression levels of each sample were estimated by RSEM (Li and Dewey 2011). Differential expression analysis was performed using the DESeq R package (1.10.1). DESeq provides statistical routines for determining differential expression of digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using the Benjamini and Hochberg approach for controlling the false discovery rate. Genes with an adjusted P-value <.05 found by DESeq were designated as differentially expressed. GO enrichment analysis of the differentially expressed genes (DEGs) was implemented by the GOseq R package based on Wallenius non-central hyper-geometric distribution (Young et al. 2010). We used KOBAS (Mao et al. 2005) software to test the statistical enrichment of the differential expression of genes in KEGG pathways.

Real-time fluorescence-based quantitative PCR (RT-qPCR)
Six genes that were possibly associated with the adaptation of low-temperature stress were selected for RT-qPCR analysis. Total RNA was quantified using a NanoDrop ND-2000 system (Thermo Scientific), and RNA integrity was assessed by agarose gel electrophoresis. RNAs of sufficient quality were treated with DNase I, followed by reverse transcription, cDNA dilution, and RT-qPCR in a 7900HT PCR system (Applied Biosystems). After completion of the reaction, data analysis was performed based on the obtained sample Ct value.

Tolerance to subzero-temperature at different germination stages
The A. aphylla seeds at different germination stages that were subjected to gradual cooling to the target temperatures exhibited variations in tolerance to low temperatures (Table 1). The A. aphylla seeds at stages I and II were tolerant to low temperatures, whereas those at stage III were damaged by low temperature, and compared with stages I and II, the damage was obvious and the survival rate decreased. As germination progressed, the tolerance of the A. aphylla seeds to low temperature further decreased, and the ability of A. aphylla seedlings (from stage IV) to survive after low-temperature stress application drastically decreased. These results indicated that the A. aphylla seeds at the early germination stages had a higher tolerance to low temperature than those at later germination stages (tolerance to low temperature: I ≈ II > III > IV). Figure 1 shows the SOD, POD, and CAT activities and the MDA content of the A. aphylla seeds at different germination stages. As germination progressed, the SOD activity of the A. aphylla seeds gradually decreased. POD and CAT activities and the MDA content initially increased, peaked at stage I, and then decreased during the later stages. In addition, no significant differences in SOD and POD activities and the MDA content in the A. aphylla seeds were observed between germination stages I and II.

Sequencing results and de novo assembly
The transcriptomes were sequenced using the Illumina sequencing system. A total of 639,729,002 raw reads were obtained from 12 samples (n = 3 for each germination stage). After eliminating the low-quality reads, 629,456,368 high-quality clean reads were obtained. The base percentage of Q > 20 was within the range of 96.39-97.33%, and the ratio of the GC content to the total number of bases was within the range of 41.56-43.19%. Through gene sequence alignment, 31,496,252-45,670,086 sequences were successfully aligned, with alignment rates ranging from 71.31% to 76.40% (Supplementary Table 1). After splicing the reads, a total of 582,580 transcripts were assembled, and the assembled transcripts were within the length range of 201-15,661 bp (average: 768 bp); the length of N50 was 1127 bp. After further splicing of the transcripts, 484,941 UniGene  Figure 1).

Statistical results of DEGS
The DEGs of W2 vs. W1, W3 vs. W2, W4 vs. W3, W3 vs. W1, and W4 vs. W1 were counted (Figure 2(A)). The number of downregulated genes was higher than the upregulated genes in W3 vs. W2. A Venn diagram was constructed and used to compare the possible relationship of downregulated genes in W1 vs. W2, W3, and W4, and 2787 downregulated genes were co-expressed in three combinations (W2 vs. W1, W3 vs. W1, and W4 vs. W1) (Figure 2(B)). In addition, compared to the gene expression in W1, the number of downregulated genes was highest in W4, followed by W3 and W2. Comparison of the downregulated genes among the four germination stages identified 104 continuously downregulated genes during germination (Figure 2(C)).

Annotation analysis of downregulated genes in different pair-wise combinations
GO enrichment analysis was performed to determine the association of the downregulated genes with specific biological processes using different pair-wise comparisons (Supplementary Table 2). The results showed that the biological processes of mRNA methylation and mRNA modification were mainly enriched in downregulated DEGs. W2 vs. W1 had relatively more downregulated genes, which significantly enriched the functional categories of metabolic process, single-organism metabolic process, and organonitrogen compound metabolic process. W3 vs. W2 had some downregulated genes that significantly enriched the functional categories of response to water, response to auxin, and response to abiotic stimulus. W4 vs. W3 had downregulated genes that significantly enriched the functional categories of lipid transport, lipid localization, and response to hormone processes. Downregulated genes in W3 vs. W1 mainly enriched the functional categories of metabolic process and single-organism metabolic process, and downregulated genes in W4 vs. W1 mainly enriched the functional categories of organic substance biosynthetic process and cellular biosynthetic process. Supplementary Table 3 shows the pathways that were enriched with downregulated genes using different pair-wise comparisons. A total of 1,261 downregulated genes were identified in W2 vs. W1, which enriched 113 pathways, mainly photosynthesis and oxidative phosphorylation pathways. A total of 1292 downregulated genes were identified in W3 vs. W2, which enriched 109 pathways, mainly spliceosome and plant hormone signaling pathways; 296 downregulated genes were detected in W4 vs. W3 that enriched 88 pathways, all of which were plant hormone signaling pathways; 3387 downregulated genes were identified in W3 vs. W1, which enriched 123 pathways, mainly photosynthesis and photosynthesis-antenna protein pathways; and 3658 downregulated genes were detected in W4 vs. W1, which enriched 122 pathways, mainly spliceosome and oxidative phosphorylation pathways.

Cluster analysis of DEGs
K-means cluster and SOM (Self-Organizing Map) cluster analyses were performed on 22,105 differentially expressed genes in this study. K-means cluster analysis divided 22,105 DEGs into 4 subclusters (Figure 3). The genes of subcluster_1 (2536 genes) and subcluster_2 (7509 genes) had the same expression pattern, and these genes were continuously upregulated in the four stages. The genes of subcluster_3 (8880 genes) and subcluster_4 (3180 genes) had the same expression pattern, and these genes were continuously downregulated in the four stages.
3.7. Gene analysis of continuously downregulated genes in W2 vs. W1, W3 vs. W2, and W4 vs. W3 Supplementary Table 4 shows that the analysis of 104 continuously downregulated genes identified only one gene in W2 vs. W1, and 12 genes in W3 vs. W2 exhibited greater than three-fold downregulated expression (log2FC: <−3). Of those genes, the BURP domain-containing protein was downregulated to the greatest extent (Log2FC: −4.427); 27 genes in W4 vs. W3 had greater than three-fold downregulated expression, and serine carboxypeptidase was downregulated to the greatest extent (Log2FC: −9.0395). By annotating the continuously downregulated genes, we found that the number of genes encoding the late embryogenesis-abundant protein was the most abundant, followed by those encoding oleosin, probable protein phosphatase 2C, beta-galactosidase, and carrot ABA-induced in somatic embryos 3 (Table 2).
3.8. Gene analysis of the downregulated genes in the W2 vs. W1, W3 vs. W1, and W4 vs. W1 combinations of comparison By comparing W1 with W2, W3, and W4, common downregulated genes in different combinations of comparison were statistically analyzed ( Table 3). The results showed many genes encoding cytochrome P450 and the BURP domain-containing protein. Similar to the results of continuously downregulated genes, this experiment also detected genes encoding oleosin, late embryogenesis abundant (LEA) protein, and probable protein phosphatase 2C Figure 4. SOM cluster analyses of differentially expressed genes. The diagram was generated using log2-transformed ratio values. The grey lines in each subcluster represent the relative expression of genes in a cluster under different experimental conditions. The blue line indicates the average of the relative expression of all genes in this cluster under different experimental conditions. The red line is the reference, the expression above the line is upregulated, and the expression below the line is downregulated. and included some genes related to low temperature responses such as dehydration-responsive protein, ethylene-responsive transcription factor, low temperatureinduced protein, and cold-regulated 413 plasma membrane proteins.

RT-qPCR results
Correlation analysis was performed between the RT-qPCR and RNA-seq data (Figure 5), and the results were consistent with a correlation coefficient range of 0.561-0.992, confirming that the data generated by transcriptome sequencing were highly reproducible and accurate.

Adaptation features of A. aphylla seeds to low temperature
Desert areas rarely receive any precipitation and are subjected to high rates of evaporation. Plants under such extreme climate conditions have adapted mechanisms for their survival . Seeds of many desert plants have the characteristics of rapid imbibition and germination (Gutterman 1993;Zhang and He 2008;Li et al. 2014), which are the bases for successful plant colonization under adverse conditions. Zhang and He (2008) showed that seeds of eight Artemisia species found in the desert steppe in China germinate rapidly, with low water consumption. In addition, seeds of Salsola kali germinate 29 min after imbibition (Gutterman 1993). Studies have shown that precipitation is an important factor for seed germination in arid regions (Garnett and Williamson 2010), which are predicted to experience increasing precipitation in the next 50-100 years, mainly in the fall and winter (Piao et al. 2010), thereby showing a trend of warming and humidification. In addition, climate change has been associated with an increase in temperature during winter, and the freezing period has become shorter, which  provides favorable conditions for seed germination in the fall. Similarly, seeds of A. aphylla undergo rapid germination when they come into contact with water. Some of the previously ripened and shed seeds in the wild sprout after autumn rainfall. Due to differences in germination time, the seeds will overwinter in different stages. In this study, A. aphylla seeds undergoing imbibition were rarely affected by low temperature and could survive the winter. As seed germination progressed, the tolerance to low temperature decreased during the later germination stages, and the seeds were vulnerable to damage caused by low temperature. A previous report evaluated the effect of low-temperature stress in three stages (testa rupture, radicle grown to 1.0 mm, and cotyledon emergence) of Lepidium apetalum seeds (Zhao et al. 2010). The result indicates that their tolerance to low temperature gradually decreases as germination progresses, which is consistent with our current findings. Li et al. (2014) showed that during seed germination, the MDA content in A. elatior at different seedling stages drastically decreased to very low levels; SOD activity initially decreased and subsequently increased; and CAT activity gradually increased. In this study, POD and CAT activities and the MDA content of A. aphylla seeds were elevated during imbibition (seed germination stage I). As seed germination progressed, POD and CAT activities gradually decreased, which is similar to that observed in A. elatior. SOD activity gradually decreased as seed germination progressed. Changes in these enzymes may have certain effects on the response to low temperature at different germination stages. MDA is toxic to plant cells. To avoid MDA cytotoxicity in plants, the activities of the protective enzymes, i.e. SOD, CAT, and POD, of A. aphylla seeds were maintained at high levels, suggesting that A. aphylla seeds are highly resistant to free radical oxidation during imbibition.
In the present study, seedlings that rapidly germinate in the fall were less likely to survive overwinter, which is also the reason why some of the A. aphylla seeds germinate in early spring. The tolerance to low temperature of A. aphylla seeds gradually decreased in the later stages of germination, with similar climatic conditions of northern Xinjiang in early spring, i.e. moist soil and gradually elevated temperature. These climatic conditions facilitate the survival of A. aphylla seedlings. In addition, numerous studies have shown that surface salinity and plant inhibitors in the arid zone slow the seed germination rate (Zhang et al. 2005;Wei et al. 2007;Ren et al. 2011). These inhibitors are eliminated when snow melts during early spring. Thus, blockage of seed germination processes of A. aphylla before freezing is beneficial to more A. aphylla overwintering propagules.  This blocking effect has important ecological significance for some desert plants in the context of climate change.

Transcriptome analysis at different germination stages under low-temperature stress
In this study, the low-temperature tolerance of A. aphylla seeds decreased as seed germination progressed. Transcriptome sequencing showed that the number of downregulated genes of A. aphylla seeds increased from imbibition to testa removal stages, and this number was higher than the number of upregulated genes, suggesting an important turning point in seed development. Compared to dry dormant seeds, the number of downregulated genes increased as germination progressed, which is possibly due to the generation of more DEGs in the late stage of seed development and an increase in DEGs involved in the negative regulation process in the late germination stages. With the germination of A. aphylla seeds, some downregulated genes related to stress responses such as response to water, response to abiotic stimulus, response to hormones, response to auxin, and response to endogenous stimulus were generated. A study has shown that when coldsensitive plants encounter low-temperature stress, disturbances in physiological processes such as water status, mineral nutrition, photosynthesis, respiration, and metabolism occur (Lukatkin 2002). The degree of chilling injury in plants is mitigated by preventing the disturbance of water conditions (Bloom et al. 2004). In addition, nitric oxide, cytokinins, auxins, and gibberellins are also reportedly involved in the responses to low temperature (Peleg and Blumwald 2011;Shi et al. 2015). This study found that the products encoded by some downregulated genes are related to material synthesis and stress responses. Among these, genes encoding LEA proteins occupy a large proportion of the downregulated genes. LEA proteins are widely existing hydrophilic proteins in plants that cope with the stress of dehydration caused by drought, salinity, and freezing. These proteins have very strong thermal stability and impart protective effects on plants against chilling, thereby preventing the loss of lactic dehydrogenase (LDH) activity (Nakayama et al. 2008). Sutton et al. (1992) was the first to show that the LEA gene is rapidly expressed in Hordeum vulgare seedlings at low temperature, and therefore, the seedlings rapidly generated stress proteins. In addition, many studies have shown that oleosin (Shimada et al. 2008), plant cytochrome P450 (Divi and Krishna 2010;Yang et al. 2014); BURP domain-containing protein (Wang et al. 2016), and cysteine proteinase inhibitors (Valdés-Rodríguez et al. 2007;Hwang et al. 2010;Zhang and Fu 2018) have a mitigation effect on low-temperature stress. Among these, plant cytochrome P450 is involved in the biosynthesis of brassinosteroids (BRs), which play an important role in increasing seedling tolerance and resistance to environmental stress. Moreover, downregulation of low temperature-related genes such as dehydration-responsive protein, heat shock protein 70, ethylene-responsive transcription factor, cold-regulated 413 plasma membrane protein, and low temperature-induced protein in plants is essential to adapting to low temperatures. Therefore, we believe that downregulation of these genes may be one of the factors that lead to a decline in low-temperature tolerance during the later stages of seed germination.
RT-qPCR analysis identified six genes that were highly expressed in A. aphylla seeds during imbibition, and the expression was significantly decreased during testa removal and elongation stages. Among the six genes, Cluster-4157.197391, Cluster-4157.284657, and Cluster-4157.287619 were extremely downregulated during the elongation stage. Cluster-4157.197391, Cluster-4157.284657, and Cluster-4157.287619 are low temperature-associated transcription factors, namely, basic helix-loop-helix (bHLH), EREBP, and C3H, respectively (Pang et al. 2013;Wang et al. 2013). With the growth of seedlings, the expression of Cluster-4157.191255 and Cluster-4157.207065, which represent the LEA protein and oleosin, respectively, significantly decreased during testa removal and was extremely low during the elongation stage. These genes may have a certain resistance to low temperature in the early stage of seed germination. As germination progressed, the expression of some genes in A. aphylla seeds decreased, which weakened the resistance of A. aphylla to low temperatures during the later stage of seed germination.

Conclusion
Seeds of A. aphylla, a desert plant, are highly adaptable to low temperature and can rapidly germinate after exposure to water during the fall. A. aphylla seeds at different germination stages have different tolerance levels to low temperature in the winter, and their freezing resistance decreases as germination progresses. Some enzymatic activities of A. aphylla seeds also decreased as germination progressed after the testa removal stages. According to the analysis of the transcriptome sequencing results of A. aphylla seeds, the downregulated expression of some genes related to stress relief and cell regulation and protection such as the LEA protein, oleosin, and plant cytochrome P450 reduces the resistance to freezing. In addition, some low temperature-related genes such as low temperature-induced protein, cold-regulated 413 plasma membrane protein, and ethylene-responsive transcription factor were also downregulated, and these genes may be very important for anti-freezing of A. aphylla seeds over winter. Therefore, to understand the survival mechanism of desert plant seeds in harsh climatic environments, this study examined the adaptability of the seeds to low temperature and provided materials for explaining the molecular mechanism of safe overwintering of non-germinating and early stage germinating seeds.

Notes on contributors
Mengwen Peng is a master's degree student, studied at the Agricultural College of Shihezi University, Xinjiang, China. Her main research direction is forest ecology.
Yaling Chang is a master's degree student, studied at the Agricultural College of Shihezi University, Xinjiang, China. Her main research direction is forest ecology.
Guangming Chu is an associate professor of ecology, focuses on forest ecology and conservation and utilization of woody plant resources. His main works include the research on the construction of artificial oasis protection ecological security system (Xinjiang, China, 2012) and the physiological and ecological responses and adaptation strategies of vegetation in arid and semi-arid regions (Xinjiang, China, 2010).
Mei Wang is an associate professor of ecology, focuses on forest ecology and conservation biology. Her main works include the research on the construction of artificial oasis protection ecological security system (Xinjiang, China, 2012) and the ecological and conservation techniques of endangered plant populations in Qinling (Xinjiang, China, 2015).