Physiological Attributes and Transcriptomics Analysis Reveal the Response Mechanism of Helictotrichon Virescens to Low-Temperature Stress


 Background: Helictotrichon Virescens is a perennial herb mainly distributed in high altitude areas of 2000 ~ 4500 m. It is widely cultivated in the Qinghai-Tibet Plateau of China, strongly resistant to cold, and an essential part of wild herbage in this region. However, the molecular mechanism of Helictotrichon virescens response to low-temperature stress is unclear, and the key regulating genes to specific biological processes are poorly known. Results: It was used physiological and transcriptome analysis to study the cold stress response mechanism in Helictotrichon Virescens. Under low-temperature treatment at 0℃, the leaves reduced chlorophyll content and light energy absorption through light antenna complex (LHCs). As the content of chlorophyll, a and b decreased, as evident were the photoprotection effects. Besides, circadian pathway-plant is crucial for the response to cold stress in Helictotrichon Virescens. Among them, DEG encoding LHY and HY5 showed strongly up-regulated expression during cold stress in Helictotrichon Virescens. Conclusions: This study provides an opportunity to fully understand the response process of Helictotrichon Virescens to low-temperature stress. It answers the pertinent questions in perennial herbaceous cold stress research, how leaves integrate light and low-temperature signals to regulate circadian rhythms in crops, and how leaf tissue reduces light absorption and enhances light protection through optical antenna complex (LHCs).

can produce higher seed and forage yield with an average annual seed yield of 0.948t/hm 2 per mu 48 and hay yield of 7.74t/hm 2 . It might be a first choice for improving natural grassland, the 49 construction of artificial grassland, and the greening of the Tibetan Plateau's ecological environment. 50 However, the molecular mechanism of the response to low-temperature stress is still unclear, and 51 the key genes regulating related biological processes are poorly understood. 52 Low-temperature stress is a multi-genic process accompanied by physiological adjustments in 53 the plant. Mechanisms of cold stress response in woody species are complicated than those in 54 of the control group. The results indicated the stability and tolerance of the leaf cell membrane 114 affected by low-temperature stress. Under low-temperature stress, the accumulation of proline (Pro) 115 could reduce the freezing point of cell solute and improve the osmotic potential, thus stabilizing the 116 cell membrane system to prevent freezing dehydration decreased solute exudation. After 12h of 117 low-temperature stress, the Pro content of leaves in the treatment and control groups showed a 118 highly significant difference (Fig.2-d). The Pro contents in the treatment group were 2.16 times 119 higher than the control group. By extending stress time, the Pro content did not change significantly. 120 After 12h of low-temperature stress, the increase of Pro content in the plant saturated. ROS content 121 gradually increased by extending stress time (Fig.2-e) with significant or highly significant 122 differences between treatment and the control groups. Simultaneously, the POD, SOD, and CAT 123 activities significantly enhanced ( Fig.2-f, g, and h), which were conducive to eliminate ROS and 124 reduced ROS damage in plants. 125 In conclusion, the low-temperature stress significantly changed the morphology, physiological 126 and cell solutes of Helictotrichon Virescens, damaged leaf cells tolerance and stability and 127 suppressed its growth. Meanwhile, low-temperature produced some of the osmotic regulation 128 hormones to alleviate low-temperature damage, but specific responses and regulation mechanisms 129 are still unclear. Therefore, it is necessary to conduct further research. 130

Sequence assembly and splicing 131
High Illumina HiSeq TM flux sequencing platform was used for transcriptome sequencing of 132 Helictotrichon Virescens leaves after low-temperature stress of 18 samples and obtained a total of 133 24000000bp raw data. After assembling and removing the redundancy, the high-quality sequence 134 data with a length of 22730, 000bp was obtained. Each sample Clean reads exceeded 6.5 G, 135 indicating a high sequencing depth. The base error rate of the sequenced column for each sample 136 was 0.02%-0.03%, Q20 and Q30 were above 98.00% and 94.00%, respectively, and GC contents 137 were above 50.21%-55.21% (Table. S1). 138 Clean reads have 396,649 transcriptional sequences obtained through Trinity splintering used 139 as reference sequences for subsequent analysis. Table.S2 shows the statistical results of reference 140 alignment. The comparison efficiency between Clean reads of 18 samples and the reference genome 141 was higher, all above 70%, meeting subsequent analysis needs. The transcriptome sequences were 142 assembled into 112,775 Unigene through Corset hierarchical clustering for subsequent correlation 143 analysis. The Transcripts of N50 and N90 of Unigene were 1579, 549, 1426, and 452, respectively; 144 larger N50 and N90 showed the overlapping of transcripts (Table. S3 and S4). 145

Gene functional annotation 146
We annotated the obtained unigenes for gene functions using seven major databases (Nr,Nt,147 Pfam, KOG, Swiss-Prot, KO, GO) ( and other species (28.1%) ( Fig. 3-a). In conclusion, Unigene's highest distribution rate in the NR 158 database, indicating that genomes of the two species were the most similar. However, a few of them 159 were in Brachypodium distachyon, Hordeum vulgare, and Triticum aestivum, indicating the 160 evolutionary overlap between the Helictotrichon Virescens and the above species. In the sequence 161 of Helictotrichon Virescens, 10.6% of genes have a sequence similarity of 95-100%, 36.1% genes 162 have a sequence similarity of 80-95%, 39.1% genes have a sequence similarity of 60-80%, and 14% 163 genes have a sequence similarity of 40-60% ( Fig. 3-b). 164

Analysis of differential gene expression 165
Based on expression profile, the gene expression level was divided into 6 grades ( Fig. 4-a). 166 0~0.1, 0.1~0.3, 0.3~3.57, 3.57~15, 15~60, and ＞60 represented very low, low, moderate, high, and 167 ultra-high and an overexpressed gene expression levels, respectively. By increasing the expression 168 level, its proportion gradually decreased. Compared with the control group, the proportion of 169 extremely low, low, and ultra-high expressed genes was higher in the treatment group. In contrast, 170 the balance of moderate, high, and overexpressed genes were more insufficient. Based on correlation 171 analysis on gene expression values (FPKM) of all samples ( Fig. 4-b), the expression patterns within 172 the group were highly similar, but differences between were greater. DEG-seq was used to screen 173 the differentially expressed genes between the low-temperature treatment and control groups (fold 174 change＞2, FDR＜0.01). Low-temperature treatment of 12 h differentially expressed 18184 genes 175

GO Enrichment of DEGs 182
To better explain the biological function of DEGs, it was performed GO enrichment analysis 183 on the differentially expressed genes after 12h, 36h, and 60h treatment, which showed that 18184 184 It was analyzed the first 20 biological processes with significantly enriched genes after 12h, 194 36h, and 60h treatment. The results showed that 11 Biological processes with differentially 195 expressed genes were co-enriched and related to carbohydrate catabolism and metabolism (Fig. S1). DEGs' expression levels in this study suggest that low temperatures themselves disrupt the circadian 237 rhythm of Helictotrichon virescens. 238

Virescens 240
In order to investigate the effect of low temperature stress on transcription factor expression of 241 Helictotrichon virescens, low temperature and normal temperature treatments were carried out, and 242 the differential expression of transcription factor genes (TFs) was analyzed after transcriptome 243 sequencing. The results showed that 2323 differentially expressed transcription factor genes were 244 obtained under low temperature stress. All transcription factor genes obtained were family classified, 245 and the results showed ( Fig. 6) that all transcription factor genes of Helictotrichon virescens were 246 classified into 84 transcription factor families, the top 10 transcription factor families were AP2/ERF, 247 FAR1, TRAF, C2H2, bHLH, B3, bZIP/GRAS, C3H, AUX/IAA and GARP-G2-like. 248 The expression multiples of transcription factor genes were analyzed, and the transcription 249 factor genes differentially expressed under low temperature stress were classified by Ration (the 250 ratio of FPKM value under low temperature stress to FPKM value under normal temperature stress) 251 as the classification index, the results showed ( Fig. 7) that the distribution of multiple changes of 252 transcription factor genes under different time treatments was basically the same, with the most 253 transcription factor genes with expression multiple less than 10 times, the number of transcription 254 factor genes contained after 12h, 36h and 60h treatment was 2242, 2246 and 2264, respectively, 255 Among them, the number of transcription factor genes whose expression multiple was less than 10 256 times was the most (2264) after 60h of low temperature treatment; After 12h of low temperature 257 treatment, the number of transcription factor genes with expression multiple of more than 100 times 258 was the highest (19); After 36h of low temperature treatment, the number of transcription factors 259 with expression multiples between 50-100 times was the highest (12). These results indicate that 260 under low temperature stress, most transcription factor genes have little changes and have little 261 influence on gene expression regulation. Only a few transcription factor genes have undergone 262 drastic changes, which have a great influence on gene expression regulation. 263 The transcription factors with large differential multiples were analyzed, and the transcription 264 factor genes with the top 15 differential expression multiples were selected for analysis under 265 different treatment times, the results showed (Table 2, 3 and 4) that all transcription factors were 266 up-regulated after 12h of low temperature treatment, Among which, 6 transcription factors were 267 AP2/ERF family, 2 were C2H2 family, and bZIP family had the highest differential expression rate; 268 After 36h of low temperature treatment, all transcription factors were up-regulated, among which 6 269 transcription factors were AP2/ERF family, 3 were C2H2 family, and bZIP family had the highest 270 differential expression rate; After 60h of low temperature treatment, all transcription factors were 271 up-regulated, Among them, 5 transcription factors were AP2/ERF family, 2 were C2H2 family and 272 2 were MYB family. As can be seen from the above transcription factor families with large multiple 273 differences, there are different types of transcription factor families with large multiple differences 274 between treatments at different times, and most of these transcription factors are related to plant 275 resistance, mainly concentrated in the AP2/ERF, MYB, bZIP and C2H2 families. 276

RT-PCR verified the results of transcriptome sequencing 285
To verify RNA-seq data's accuracy, we randomly selected 6 DEGs involved in circadian 286 rhythm, Plant hormone signal transduction, Photosynthesis and Photosynthetic -antenna protein 287 pathway for RT-PCR. The up-regulated or down-regulated DEGs in qPCR were consistent with 288 RNA-seq analysis results (Fig. 8). Simultaneously, plants have ROS scavenging mechanisms, including enzymatic and non-enzymatic 307 systems in which antioxidant enzymes are important ROS scavengers. In this study, ROS content 308 was gradually increased by extending treatment time (Fig. 2-e). Specifically, after 12h and 36h of 309 low-temperature stress, ROS contents were significantly different between treatment and control 310 groups. After 60h of low-temperature stress, ROS content was significantly different between the 311 treatment and control groups. 312 Meanwhile, low-temperature stress up-regulated the genes regulating oxidative stress response 313 (Cluster37118.30992) (Fig. 5-a). Besides, by extending treatment time, antioxidant enzymes such 314 as POD, SOD and CAT increased gradually. Specifically, after 12h of low-temperature stress, the 315 contents of POD, SOD and CAT were significantly different between the treatment and control 316 groups. After 36h and 60h of low-temperature stress, POD, SOD, and CAT contents were 317 significantly different between the treatment and control groups. Simultaneously, under low-318 temperature stress, the gene Cluster37118.16911 regulated Peroxidase activity (POD), and the gene 319 Cluster37118.62042 up-regulated the response to oxygen-containing compound (SOD, etc.), which 320 was consistent with the results of physiological experiments (Fig. 2-f, g and h). These results 321 indicated that POD, SOD and CAT were beneficial to ROS clearance and reduced ROS damage in 322 plants during low temperature and cold stress. with low-temperature stress for 12h (Fig. 2-a and b). After 36h and 60h of low-temperature 343 treatment, the content of chlorophyll a and b was significantly different between the treatment and 344 control groups. Besides, the genes encoding light-capture protein complexes, such as Lhca1、 Lhca2、 345 Lhca3、Lhca5、Lhcb1、Lhcb3、Lhcb4、Lhcb5、Lhcb6 and Lhcb7, were up-regulated after cold 346 treatment. Cold treatment up-regulated the genes encoding PSI subunit (psaA、psaD and psaE) and 347 the genes encoding PSII subunit (psbA、psbC and psbO) (Fig. 5-b). These results indicated that 348 leaves tried to reduce light absorption by reducing chlorophyll content during low-temperature 349 treatment. As much as chlorophyll a and b content decreased, the light protection effect became 350 noticeable. Therefore, it is proposed that photoinhibition is a temporary or persistent form of 351 photoprotection, which plays an important role in the low-temperature stress tolerance in 352 Helictotrichon Virescens. 353

Circadian rhythm plays an important role during cold stress of Helictotrichon Virescens 354
Studies have shown that the circadian clock is the gateway to induce genes related to cold stress. Helictotrichon Virescens (Fig. 5-d). Therefore, it is speculated that the circadian rhythm pathway is 363 important to the low-temperature response of Helictotrichon Virescens. 364 Virescens were classified into 84 transcription factor families (Fig. 6), the top 10 transcription factor 374 families were AP2/ERF, FAR1, TRAF, C2H2, bHLH, B3, bZIP/GRAS, C3H, AUX/IAA and 375 GARP-G2-like, most of the transcription factor genes belong to the transcription factor family 376 related to resistance. According to the analysis results of multiple differences of transcription factor 377 genes in this study (Fig. 7), transcription factor families with greater differential expression also  This study uses Helictotrichon Virescens as experimental material. Seeds were sown in 432 flowerpots (21cm×16cm) containing matrix (peat, pine needles, and yellow clay, volume 3:1:1). 433

Transcriptional regulation of transcription factors is a key part of the response of
When the first leaves expanded fully, we selected eight strains of each material with three 434 replications and transferred them into 1/2 Hogland nutrient solution, then collected the samples from 435 5 weeks old seedlings for the relevant experiment. The collected leaf tissues were immediately 436 frozen in liquid nitrogen on the sampling day and stored at -80℃. 437

Low temperature and cold stress treatment on Helictotrichon. Virescens plants 438
After separating five weeks old seedlings into control and treatment groups, it was treated the 439 seedlings with low temperatures stress in a constant temperature and light in an incubator (MODEL: 440

MLR-352H-PC) of the Cell Genetics Laboratory of Maize Research Institute of Sichuan 441
Agricultural University by following complete block design (CBD). The low-temperature and 442 control groups constantly remained at 0℃ and 25℃, respectively, in an illumination set to 3000Lx. 443 seq and related indexes 445

Relative conductivity (REC) 446
The relative conductivity was measured by referring to the method of Xu [28]. The specific 447 procedure was to take a 2×4 cm leaf from the middle of the first fully expanded leaf of each seedling. 448 Then, mixed the leaves, cut the leaves into 1 cm segments, divided them into 3 portions, transferred 449 them into 10 mL EP tubes, and finally filled the tubes with distilled water. After soaking the blades 450 in water for 3h, the EC1 was determined by the conductivity meter; then, boiled the samples in a 451 water bath for 10min, followed by cooling them at room temperature and measured EC2 by a 452 conductivity meter. 453 Calculation formula of relative conductivity: REC=EC1/EC2×100%. 454

Determination of ROS content and Antioxidant enzyme activity 465
The test kit purchased from Suzhou Keming Biotechnology Co., LTD, was used to determine 466 Proline (Pro), Superoxide dismutase (SOD), and Peroxidase (POD), Catalase (CAT) and Reactive 467 oxide species (ROS). Among them, the determination instrument is UV-visible visible 468 spectrophotometer. 469

RNA extraction and detection 470
Beijing Novogene Technology Co., LTD completed RNA extraction and detection. 471

cDNA library construction and sequencing 472
The Beijing Novogene Technology Co., LTD constructed the cDNA library and sequenced it 473 by adopting the Oligo(dT) magnetic beads enriched mRNA with polyA tail procedures. Briefly, 474 mRNA was randomly fragmented by divalent cation in NEB Fragmentation Buffer. The first strand 475 of cDNA was prepared using random oligonucleotide as a primer, and the second strand was 476 synthesized using DNA polymerase I. The ends of purified double-stranded cDNA were repaired by 477 adding a tail and sequencing connector and screened 250-300 bp cDNA for PCR amplification. The 478 PCR product was purified again to obtain the final library, and after qualifying the library, finally 479 used Illumina Hiseq TM 4000 sequencing platform for sequencing. 480

Data filtering and splicing assembly 481
High-quality sequence data was obtained by filtering the raw reads to a small number of reads 482 containing sequencing adaptors, and which did not show the base information were considered low-483 sequencing reads. Through statistical Q20 (Phred = -10log10 (e)), Q30 (Phred = -10log10 (e)), GC 484 content and sequencing error (not more than 6%), and other indicators of sequencing quality control. 485 Trinity software [30] was used to obtain clean reads of transcripts for later analysis. The Tgicl 486 software was used to obtain Unigenes by homologous clustering for later analysis. Sequence the 487 splicing transcript according to its length from long to short, and add the length of the transcript to 488 the length of the splicing transcript that is not less than 50%/90% of the total length, namely 489 N50/N90, to measure the continuity of De Novo assembly, and its numerical value can be used to 490 evaluate the quality of the assembly effect. According to statistical analysis, the level of gene expression screening of differentially 502 expressed genes between different samples, DEG-Seq software [34] to the original read count for 503 standardization, through the negative binomial distribution hypothesis test probability statistical 504 model (P-value) calculation, and finally (BH) adjustment for multiple hypothesis testing, by FDR 505 value (False Discovery Rate, namely padj), |log2(FoldChange)| > 1 and padj < 0.05 was considered 506 as the standard for differential gene screening. 507 GO biological enrichment analysis was performed on the differential gene set by the GO-seq 508 method [35], and KEGG metabolism and signal transduction pathway enrichment analysis was 509 performed on the differential gene set by KOBAS method [36]. Padj < 0.05 was used as the 510 threshold of significant enrichment in both of the above two analyses. 511

Real-time PCR validated the transcriptome results 512
In order to verify the accuracy of transcriptome sequencing results, six genes were randomly 513 selected for real-time PCR detection. Primer5.0 was used to design differential gene sequence 514 primers (Tables. S9) Table 1 Statistics of gene annotation 548 Table 2 Transcription factors with large difference multiple under low temperature stress (after 12h 549 treatment) 550 Table 3 Transcription factors with large difference multiple under low temperature stress (after 36h 551 treatment) 552

Availability of data and materials 570
The datasets supporting the conclusions of this article are included within the article and its 571 additional files. The original data of RNA sequences will be uploaded to the NCBI database and 572 relevant login information will be provided. 573

Declarations 574
Ethics approval and consent to participate 575 The experimental materials for this study were cultivated by Sichuan Grass Industry Technology 576 Research and Promotion Center, Grassland Station of Ganzi Prefecture, Sichuan Province and 577 Sichuan Jinzong Liaoyuan Seed Industry Technology Co. Ltd and without any controversy, it 578 comply with relevant institutional, national, and international guidelines and legislation. 579

Consent for publication 580
Not applicable. 581  The abscissa is log10 of sample 1 (FPKM+1), the ordinate is log10 of sample 2 (FPKM+1), and R2 is the square of Pearson's correlation coe cient. c.d and e: Volcano map of differential genes. Note: In the gure, the abscissa is Log2 Fold Change value, and the ordinate is -Log10 padj or -Log10 p-value. The dashed blue line represents the threshold line of differential gene screening criteria. f: Venn diagram of different genes after 12h, 36h and 60h cryogenic treatment.    Fold change distribution of transcription factor genes Note: a, b and c are the expression multiples of transcription factor genes after 12h, 36h and 60h low temperature treatment, respectively. Photographs of Helictotrichon virescens