Identification Of Arabidopsis genes associated with cold tolerance based on integrated bioinformatics analysis

ABSTRACT Cold stress is a major environmental factor that limits plant growth and productivity. Plants have evolved various strategies to adapt to these environmental conditions. To better explain the mechanisms used to survive environmental challenges, we retrieved the cold-responsive genes of Arabidopsis thaliana from the Gene Expression Omnibus (GEO) database. The GEO raw data were normalized by the quantile method, and then the differentially expressed genes (DEGs) under cold stress were screened using the robust rank aggregation (RRA) algorithm, including 261 upregulated and 177 downregulated genes out of more than 20,000 genes. Further, the integrated bioinformatics analyses of PUBMED, PANTHER, DAVID, and STRING indicated that the upregulated DEGs were involved in cellular response to red light, negative regulation of circadian rhythm, photoprotection, monosaccharide transport, cold acclimation, and phosphate ion homeostasis, while the downregulated DEGs were associated with the indole glucosinolate biosynthetic process, regulation of RNA splicing, water transport, cell wall modification, cell wall loosening, cellular water homeostasis, and cell wall homeostasis. Furthermore, the up-regulated DEGs had about four times protein-protein-interactions (PPIs) than the down-regulated DEGs, and the cold-responsive genes were identified using Cytoscape software. Furthermore, qRT-PCR of low-temperature-responsive protein 78 (LTI78), transducin family protein (SWA1), and arginine methyltransferase 11 ( PRMT11) were performed to validate the outcome of integrated bioinformatics analysis. Our work will improve our knowledge of cold-responsive mechanisms and these DEGs might be targets for plant cold stress-resistance research.


Introduction
Extreme low temperatures negatively affect crop productivity and threaten food security. Plants can perceive extreme temperatures and adjust their growth, reproduction, and development to adapt them. To better understand plant cold-resistance mechanisms and increase crop yield, Arabidopsis thaliana has been used as an ideal model plant to explore the response of differentially expressed genes (DEGs) to cold stress (Gong et al. 2020).
In previous studies, several important Arabidopsis genes and pathways, including calcium signaling, mitogen-activated protein kinase (MAPK) signaling, jasmonic acid (JA) signaling, and C-repeat/DREB binding factor (CBF) signaling, have been reported to be involved in the cold response (Shi et al. 2018;Wu et al. 2019;Yuan et al. 2018). For example, glutamate receptor family proteins (ATGLR1.2 and ATGLR1.3) are glutamate-like receptors related to cold stress. Overexpression of ATGLR1.2 and ATGLR1.3 could enhance the gene expression of the CBF signaling pathway, which positively improves cold adaption in Arabidopsis (Zheng et al. 2018). The β-expansin gene (TaEXPB7-B) is another cold-responsive gene located in the Arabidopsis cell wall. Overexpression of TaEXPB7-B improved cellulose and lignin content, and increased antioxidant activity to survive at low-temperatures (Feng et al. 2019). The gene open stomata 1 (OST1) also plays an important role in enhancing low-temperature tolerance in Arabidopsis, and interacts with a plasma membranelocalized clade-E growth-regulating 2 (EGR2) phosphatase to facilitate better adaption of plants under cold stress (Ding et al. 2019). Abscisic acid (ABA) is associated with plant adaptation to cold stress. Overexpression of rice pyrabactin resistance-like gene 3 (OsPYL3) in Arabidopsis could increase the sensitivity of the ABA signaling pathway. As a result, the plant cold stress was enhanced (Lenka et al. 2018). Arabidopsis thaliana DEAD-box RNA helicase 7 (AtRH7) is a DEAD-box RNA helicase that interacts with cold shock domain protein 3 (AtCSP3), which plays a considerable role in cold tolerance. Mutants of AtRH7 negatively affect pre-rRNA processing and cause delay in first leaf emergence in Arabidopsis . Ubiquitin-conjugating enzyme 13 (UBC13) is an important cold-responsive gene that participates in programmed cell death pathways in Arabidopsis. The UBC 13 mutant was involved in the lesion mimic phenotype and interacted with F-box and associated interaction domains-containing protein 1 (CPR1), which is an F-box protein that regulates TIR-NBS-LRR class disease resistance protein 1 (SNC1) degradation under cold stress . Arabidopsis cystatin A (AtCYSa) and cystatin B (AtCYSb) are cysteine proteinase inhibitors that can be induced by cold stress. Their promoter regions included a dehydration-responsive element (DRE) and abscisic acid-responsive element (ABRE). Therefore, AtCYSa and AtCYSb are target genes of the dehydration response element B1A (DREB1A) and AREB. Overexpression of AtCYSa and AtCYSb enhances plant tolerance to environmental stress (Zhang et al. 2008). By regulating the expression of cold-responsive salicylic acid and bZIP transcription factor family protein (TGA), the regulatory protein (NPR1) also plays a vital role in enhancing the cold acclimation of Arabidopsis. Interacting with HSFA1 significantly promoted cold tolerance in plants, and NPR1 might be an essential gene during plant cold acclimation (Olate et al. 2018).
Although cold-responsive genes have been reported by different laboratories, more comprehensive DEGs will better disclose the mechanisms of cold acclimation in plants.
Recently, 5742 genes were differentially expressed to reveal the gene regulators and pathways involved in cold tolerance in Brassica napus. These DEGs were related to the inhibition of photosynthesis and the primary biological processes (Ke et al. 2020). To identify early responsive events in Oryza sativa under cold stress, a transcriptome profile was performed to identify 516 DEGs that were involved in Ca 2+ and ROS-mediated signaling, and the DREB/CBF pathway (Dasgupta et al. 2020).
In the present study, to disclose the details of cold-responsive DEGs in Arabidopsis, the Gene Expression Omnibus database (GEO) was used to retrieve the DEGs of plant cold tolerance. Then, 261 upregulated and 177 downregulated genes from the GEO database were identified using a rank aggregation method. Integrated bioinformatics strategies have been applied to demonstrate the molecular functions, signaling pathways, and PPI interaction network of cold-responsive DEGs. These DEGs will greatly improve our knowledge of plant cold-response mechanisms and might become useful targets for crop plant growth, development, and crop output.

Retrieval of cold-responsive genes in Arabidopsis from GEO database
The eight responsive gene expression profiles of the microarray and RNA-seq data related to cold stress in Arabidopsis (Table 1) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Microarray data normalization and robust rank aggregation (RRA) algorithm in Arabidopsis
Microarray data were downloaded from the GEO database in TXT format (https://www.ncbi.nlm.nih.gov/geo/). The R software package was used to process the matrix files and filter low-quality data. The resultant data were log2 transformed and processed with the limma package (http:// www.bioconductor.org/) to retain the DEGs that had a pvalue < 0.05, and a |log2fold change (FC)|> 1. Furthermore, the DEGs identified above were integrated into the robust rank aggregation (RRA) software package (https://cran.rproject.org/web/packages /RobustRankAggreg/index.html). Using a null hypothesis of uncorrelated inputs, the RRA algorithm defines genes ranked consistently better than expected, and then assigns a significance score to each gene. The hypothesis of the RRA method is that each gene is in a random order in each experiment. If a gene highly ranked in all experiments, its p-value is smaller (Niu et al. 2019).

Cluster analysis of cold-responsive DEGs in Arabidopsis
The resultant cold-responsive DEGs were clustered by Gene Cluster 3.0, using the C Clustering Library version 1.24 created by Michael de Hoon, Seiya Imoto and Satoru Miyano, and the results were viewed using Java TreeView software created by Alok (Version:1.0.4; http://jtreeview.sourceforge. net).

Gene ontology classification and KEGG Pathway in Arabidopsis
The bioinformatics tools of the PANTHER (Protein Analysis THrough Evolutionary Relationships) classification system (Version 15.0 released 2019_04) (http://pantherdb.org/) and DAVID (Database for annotation, visualization and integrated discovery, Version 6.8) (https://david.ncifcrf. gov/) were applied for gene ontology enrichment analysis of the DEGs. Each gene was classified into a single category.

Protein-protein-interaction (PPI) networks in Arabidopsis
PPI networks were analyzed using the STRING (search tool for recurring instances of neighbouring genes) database (Version 11.0, released January 19, 2019) (http://string-db. org/). Subsequently, the maximal clique centrality (MCC) app in Cytoscape software was used to screen modules within the PPI network with the default parameters.

Growth conditions and harvest of Arabidopsis
According to a previous publication (Guo et al. 2020b), ecotype Col-0 Arabidopsis seeds were germinated on a normal medium plate at 22/20°C day/night, an 8/16 h light/ dark cycle, and 60 μmol·m −2 s −1 of light intensity. After 7 days, the A. thaliana seedlings individually underwent cold stress by exposure to a temperature of 4°C for 12 h.
Then, the Arabidopsis roots and leaves were harvested and stored at -80°C.

qRT-PCR analysis
According to the previous method of RNA isolation and qRT-PCR analysis (Guo et al. 2020b), the relative expression levels of Arabidopsis gene low-temperature-responsive protein 78 (LTI78), transducin family protein (SWA1), and arginine methyltransferase 11 (PRMT11) were determined and Arabidopsis gene actin 2 (AT3G18780) was used as an endogenous reference in Arabidopsis roots and leaves. Briefly, total RNA was extracted using the RNAiso Plus reagent (TaKaRa). RNA concentrations were determined by spectrophotometry (NanoDrop 2000/2000C, Thermo Scientific). Then, 2 μg of total RNA was reverse transcribed using ReverTra Ace (TOYOBO). Real-time RT-PCR analysis was performed in a Roter-Gene Q (QIAGEN) using the Platinum SYBR Green qPCR SuperMix-UDG kit (Life Technologies Corporation). Cycling conditions were as follows: 50°C for 2 min, 95°C for 5 min, followed by 40 cycles of 95°C for 10 s and 60°C for 45 s. The 2-ΔΔCt calculation was used to determine the differences in the cold-responsive genes. All experiments were performed in triplicate (Table 2).

Statistical analysis
The data of qRT-PCR analysis were statistically analyzed by means of an unpaired t-test using GraphPad Prism 7 (GraphPad Prism, La Jolla, CA, USA). Statistical significance was set at p < 0.05.

Quantification of cold stress-responsive DEGs in Arabidopsis
These cold-responsive gene datasets were screened using the limma package in R with the following condition: p-value < 0.05 and |log 2 FC| > 1 (Figure 1). In GSE3326, 2193 downregulated and 2452 upregulated DEGs were screened in Arabidopsis. In GSE31837, 1872 downregulated and 1959 upregulated DEGs were screened in Arabidopsis. In GSE39090, 402 downregulated and 466 upregulated DEGs were screened in Arabidopsis. In GSE43818, 875 downregulated and 843 upregulated DEGs were identified in Arabidopsis. In GSE55907, 608 downregulated and 689 upregulated DEGs were screened in Arabidopsis. In GSE106635, 367 downregulated and 485 upregulated DEGs were identified in Arabidopsis. In GSE112389, 677 downregulated and 906 upregulated DEGs were identifie in Arabidopsis. In GSE113547, 1003 downregulated and 962 upregulated DEGs were screened in Arabidopsis. These cold-responsive DEGs were further clustered using Cluster 3.0 software, and viewed using Java TreeView software (Figure 2).

Cold-responsive DEGs screening by RRA algorithm in Arabidopsis
After the RRA method was used to integrate the above DEGs, 261 upregulated and 177 downregulated genes were

qRT-PCR analysis of cold-responsive genes in Arabidopsis
Our results showed that cold-responsive Arabidopsis genes LTI78, SWA1, and PRMT11 were upregulated under cold stress. Therefore, the relative mRNA expression of Arabidopsis genes LTI78, SWA1, and PRMT11 in both Arabidopsis roots and leaves were assayed under cold stress to validate the confidence of the bioinformatic outcome. The results indicated that the relative mRNA expression levels of Arabidopsis genes LTI78, SWA1, and PRMT11 in Arabidopsis roots and leaves were all increased under cold stress (Figure 7).

Discussion
Cold stress negatively affects plant growth and crop yield. Arabidopsis has been used as a model plant for research on cold-responsive mechanisms. Our previous work indicated that there are complex signaling pathways and PPI interaction networks under salt stress in Arabidopsis (Guo et al. 2014;Guo et al. 2019a;Guo et al. 2020a;Guo et al. 2020b). However, it is still not clear how plants respond to other stresses, such as low temperatures, and further research is still needed.
Here, we retrieved the cold-responsive DEGs from the GEO database. To integrate the DEGs from the eight datasets, an RRA method was executed according to a previous publication (Niu et al. 2019). As a result, 261 upregulated and 177 downregulated genes were screened using the RRA algorithm. The integrated bioinformatics methods indicated that the DEGs were involved in biological processes, KEGG pathways, and PPI interaction networks. The upregulated DEGs were involved in the negative regulation of circadian rhythm, photoprotection, monosaccharide transport, cold acclimation, and phosphate ion homeostasis, while the downregulated DEGs were mainly involved in regulation of RNA splicing, water transport, cell wall modification, cell wall loosening, cellular water homeostasis, and cell wall homeostasis. Among the upregulated DEGs under cold stress, zinc finger protein ZAT12 (ZAT12) is a protein that is necessary for the expression of ascorbate peroxidases (APXs), maintaining the balance of ROS (Rizhsky et al. 2004). In the ABA signaling pathway, ABI3/VP1 1 (RAV1) is activated by serine/threonine-protein kinase (SRK2), which then downregulats the expression of ABI3, ABI4 and ABI5. In the present study, RAV1 was upregulated, which was consistent with previous findings (Feng et al. 2014) and might play an important role in response to cold stress.
The ICE-CBF-COR pathway is an important signaling pathway associated with low-temperature stress (Guo et al. 2019b;Jin et al. 2018). In the present study, the CBFs, inducer of CBF expression (ICEs), and cold-responsive (COR) genes were differentially expressed under cold stress. Jasmonate zim-domain protein 1 (JAZ1) also increased in response to low temperature. It can inhibit the JA signaling pathway and regulate ICE1 expression (Hu et al. 2013). Furthermore, JAZ1 interacted with ICE1 and regulated the expression of the CBF. CBF is regulated by the circadian clock late elongated hypocotyl (LHY) gene, related to plant circadian rhythm. LYH positively binds to the promoter of CBF. Furthermore, under cold stress, osmotically responsive gene 1 (HOS1) ubiquitinated ICE1, and SUMO E3 ligase (SIZ1) sumoylated ICE1. Then, the resultant ICE1 was degraded by the 26S proteasome pathway (Dong et al. 2006;Dong et al. 2011;Miura et al. 2007). Additionally, cold-regulated protein 28 (COR28) interacts with protein CCA1 (CCA1) and negatively regulates the expression of CBF. In the present study, JAZ1, ICE1, dehydration-responsive elementbinding protein 1A/1B/1C (CBF1/2/3), LHY, CCA1, and COR27 were all overexpressed and interacted with each other, which played an important role in adapting to low temperature (Li et al. 2016). Furthermore, to compare with a previous publication, the DEGs involved in the CBF pathway under cold stress could also be found in rice through comparative transcriptome analysis (Dasgupta et al. 2020), which indicated the confidence and applicability of our results.

Conclusion
This work provides a new understanding of the details involved in tolerance of Arabidopsis under cold stress, which showed new signaling pathways, more cold-responsive DEGs, and more comprehensive interaction networks. This study has been helpful in demonstrating how plants survive under low temperature, and the mechanisms involved in cold tolerance might be potential targets for the research on cold-response in plants.