Novel function of a putative TaCOBL ortholog associated with cold response

The plant COBRA protein family plays an important role in secondary cell wall biosynthesis and the orientation of cell expansion. The COBRA gene family has been well studied in Arabidopsis thaliana, maize, rice, etc., but no systematic studies were conducted in wheat. In this study, the full-length sequence of TaCOBLs was obtained by homology cloning from wheat, and a conserved motif analysis confirmed that TaCOBLs belonged to the COBRA protein family. qRT-PCR results showed that the TaCOBL transcripts were induced by abiotic stresses, including cold, drought, salinity, and abscisic acid (ABA). Two haplotypes of TaCOBL-5B (Hap5B-a and Hap5B-b), harboring one indel (----/TATA) in the 5′ flanking region (− 550 bp), were found on chromosome 5BS. A co-dominant marker, Ta5BF/Ta5BR, was developed based on the polymorphism of the two TaCOBL-5B haplotypes. Significant correlations between the two TaCOBL-5B haplotypes and cold resistance were observed under four environmental conditions. Hap5B-a, a favored haplotype acquired during wheat polyploidization, may positively contribute to enhanced cold resistance in wheat. Based on the promoter activity analysis, the Hap5B-a promoter containing a TATA-box was more active than that of Hap5B-b without the TATA-box under low temperature. Our study provides valuable information indicating that the TaCOBL genes are associated with cold response in wheat.


Introduction
Low temperature stress is an important factor affecting wheat growth, development, and distribution. Cold injury leads to the decrease of wheat yield and quality. The cold resistance of plants mainly relies on the changes in morphological characteristics and internal physiological and biochemical traits. Wheat cold hardiness involves two important evolutionarily adaptive mechanisms: vernalization [1] and cold acclimation [2]. Vernalization is the process that allows winter-habit plants to synchronize flowering with favorable spring conditions after exposure to an extended time starting at low temperature [3]. There are four VRN loci in wheat: VRN1, VRN2, VRN3, and VRN4. A study suggests that VRN1 can downregulate some genes related to cold acclimation [4,5]. Under cold acclimation, VRN1 inhibits the expression of WCBF-2, TACBF-A22, and TACBF-B22 genes [3,6]. Expression analysis also has shown that VRN1 directly regulates CBF genes and represses their expression [7]. Cold acclimation enables plants to acquire freezing tolerance in response to a period of exposure to low temperatures prior to the onset of 0 °C [8]. In wheat, transcription factors involved in cold acclimation mainly include CBF/ DREB, ICE, BZIP, WRKY, MADS-box, and NAC/MYB. The cold signal is transduced downstream, and many pathways are activated in concert with the precise regulation of expression of cold-responsive genes, which participate in abscisic acid (ABA)-dependent or ABA-independent signaling pathways [9,10]. The ICE-CBF-COR pathway has been suggested to be functionally conserved in different plant species [11,12]. Cold resistance of wheat mainly involves two gene loci, Frost resistance locus 1 (Fr-1) [13] and Fr-2 [14], which are located on the long arm of homologous group 5 chromosomes.
Besides, under low temperature stress, plants will adjust their morphological structure to adapt to changes in the external environment. Research has indicated that the morphological structure of plants is closely related to their cold resistance capacity [15,16]. Plant cell wall form the first defense line upon different stresses. Many studies have shown that cell wall plays crucial roles in establishing the morphology of plant cells, including in plants' response to cold [17,18]. One of the genes identified to regulate secondary cell wall biosynthesis and determinant of the orientation of cell expansion is COBRA [19,20].
The COBRA (COB) and COBRA-like genes (COBLs) encode plant-specific glycosylphosphatidylinositol (GPI)anchored proteins, which harbor a hydrophobic signal peptide at the N-terminus, a hydrophilic middle, a highly hydrophobic C terminus, a cysteine-rich CCVS motif, and an ω-cleavage site for GPI processing [19]. The COBRA protein is predicted to be anchored to the external plasma membrane leaflet by a GPI moiety. There are 12 members of COBL genes in Arabidopsis [21], 11 in rice [22,23], and 9 in maize [24,25]. In Arabidopsis, AtCOBL4 is required for cellulose biosynthesis in the secondary wall [20], and AtCOBL9 plays a role in root hair development [26]. In rice, the OsBC1L family genes have been shown to be involved in cellulose biosynthesis and deposition of the secondary cell wall components [27][28][29][30]. In addition to the effects on cell wall biosynthesis, some of the OsBC1L members have various impacts on plant growth, including root hair development, plant height, young panicles, and pollen development [23]. ZmBK2 has been confirmed to affect the secondary cell wall formation [24,25]. However, the maize roothairless3 gene that encodes a putative GPIanchored COBRA-like protein plays an important role in root hair elongation and grain yield [31].
On the other hand, as a GPI-anchored protein located at the cell wall-plasma membrane interface, COBRA follows a typical GPI secretion pathway and can sense information from the cell wall and then transmit it to the cell membrane [32]. It is of high possibility that the COBRA family members play an important regulatory role in response to stresses and hormones. The objectives of this study were to: (i) identify the natural variations in TaCOBL, (ii) develop functional markers based on the variation, (iii) identify functional alleles and haplotypes of TaCOBL using association analysis, and (iv) analyze function of the TaCOBL. The flow diagram of different steps was showed in Supplementary Fig. S1.

Plant materials
Wheat (Triticum aestivum) cv. Annong 0711 was used for gene cloning and expression analyses. A set of Chinese Spring (CS) nulli-tetrasomic lines on homologous group 5 chromosomes were employed to verify the chromosomal location of TaCOBL. 20 accession with different cold-resistance capabilities were selected for polymorphism detection and functional marker development (Table S1, underlined varieties). The wheat natural population consisting of 336 wheat accessions from different ecological zones of wheat production were used for the validation of the newly developed sequence-tagged-site (STS) marker, as well as the association analysis between phenotypes and genotypes ( Table S1).

Strategy for cloning and sequence analysis of TaCOBL
The coding regions of rice BC1 (GenBank accession number: AY328910) [23] and maize BK2 (GenBank accession number: DQ139874) [24] were used for BLAST searches against the wheat EST database in GenBank and the wheat nucleotide database of the European bioinformatics institute (EBI) [33]. The constructed contigs were assembled to a putative full-length sequence of TaCOBL using the software DNAman. Primers TaCOBL5A, TaCOBL5B, and TaCOBL5D (Table S2) designed according to the assembled full-length TaCOBL sequence were tested using the genomic DNA of Annong 0711. The product was sequenced to confirm its specificity. Each primer combination was also tested with the CS nulli-tetrasomic lines N5A-T5B, N5A-T5D, N5B-T5A, and N5D-T5B to confirm a homologous group 5-specific location.
Genomic DNA (gDNA) was extracted from young leaves of 10-day-old seedlings using the phenol chloroform method. Primers were designed by the software Primer Premier v5.0 (Premier Biosoft International, Palo Alto, CA, USA), and all primers were synthesized by Sangon Biological Engineering Technology & Service Co., Ltd. (Shanghai, China). PCR reactions were performed in a volume of 20 μl, including 12.8 μl of ddH 2 O, 2.0 μl of 10 × PCR buffer (containing Mg 2+ ), 0.5 μl of 10 μM primer F, 0.5 μl of 10 μM primer R, 0.2 μl of Easy Taq DNA polymerase (5 U/μl), and 100 ng of gDNA. The PCR procedure was 95 ℃ for 3 min, followed by 36 cycles of 95 ℃ for 30 s, annealing (55-65 ℃) for 30 s, and extension at 72 ℃ (30 s to 3 min), with a final extension at 72 ℃ for 10 min. The annealing temperatures and extension times depended on the primer set and the length of the expected PCR product. To guarantee sequence accuracy, the PCR reactions and DNA sequencing were repeated at least three times. Multiple sequence alignment of the COBL family proteins was carried out using DNAman. Phylogenetic analysis was performed using the MEGA5 with adjacency method. The regulatory element analysis of promoter region by the PLACE tool (http:// www. dna. affrc. go. jp/ PLACE/). Prediction analysis of protein interaction network was performed by String (https:// cn. string-db. org/).

Abiotic stress treatment
Seeds of wheat cv. Annong 0711 were grown hydroponically in a light chamber at 25 ℃. Ten-day-old seedlings were subjected to various abiotic stresses. Seedlings were exposed to air on filter paper for induction of rapid drought conditions. Cold treatments were performed by transferring the seedlings to a pre-cooled medium in a growth chamber at 4 ℃. For drought, salt, and ABA treatments, seedlings were transferred to growth media supplemented with 20% PEG6000, 100 mM NaCl, and 100 μM ABA, respectively. Leave and root tissues were collected separately at 0, 1, 2, 6, 12, and 24 h after treatment; they were frozen immediately in liquid nitrogen and stored at − 80 ℃ for RNA extraction.
For real-time PCR analysis, Total RNAs were extracted with Plant RNA Reagent (Invitrogen) and used to synthesize cDNA with a Reverse Transcription System (Promega). The expression analysis of TaCOBL was performed on a cycler apparatus (Bio-Rad) with using the SYBR Green PCR master mix kit (TaKaRa Biotechnology Co., Ltd., Dalian, China; Product Code: DRR041A). The primer sets TaCOBL-2 and Actin-p (Table S2) were used for amplification of the TaCOBL-5B and actin genes, respectively. Three replicates were performed to obtain the average and standard deviation of the expression level. Quantification of the target gene expression was carried out by the comparative CT method [34].

Development of a functional marker for TaCOBL-5B
Twenty cultivars (Table S1) were used for screening the newly developed STS markers. Based on the sequence divergence of the two haplotypes at the TaCOBL-5B locus characterized in this study, two primer sets (TaCOBL-3, Table S2) were designed to identify allelic variants of the TaCOBL-5B gene in the 20 cultivars. The phenotypic variation was assumed to be associated with TaCOBL-5B, which was then validated using 336 wheat accessions. PCR products were separated using a high throughput tilling analysis system (Fragment Analyzer) to distinguish fragment lengths.

Nicotiana benthamiana transformation and promoter activity analysis
A 1608-bp fragment upstream of the start codon of TaCOBL-5B was amplified using the primers (TaCOBL-4, Table S2) containing the HindIII and BamHI restriction sites. The PCR product was confirmed by sequencing and fused to the upstream of the GUS gene in pCAMBIA1391 by replacing the CaMV 35S promoter to create new constructs, namely pCob1608a (Hap5B-a:GUS) and pCob1608b (Hap5B-b:GUS). The constructs were then transformed into Agrobacterium tumefaciens strain EHA105 through electrotransformation. Recombinant Agrobacterium was infiltrated into the young tobacco (Nicotiana benthamiana) leaves according to the described method [35].

Cloning and sequence analysis of TaCOBL
Ten wheat ESTs induced by abiotic stresses were detected through a BLAST search, which were assembled into three contigs based on their sequence similarities. Contig 1 comprised of ESTs BE500136, JZ886730, CA629775, HX137109, and HX137138 was located on chromosome 5A; contig 2 including CJ552171 and JZ885802 was mapped on chromosome 5B; contig 3 consisting of ESTs BJ258229, BU100163, and CJ691193 was mapped on chromosome 5D (Table S3). The three homoeologous TaCOBL sequences from wheat subgenomes A, B, and D, respectively. Sequence alignment of gDNA to their corresponding cDNA revealed that no intron presented in TaCOBLs, with an open reading frame (ORF) length of 1,350 bp for TaCOBL-5A, 1,353 bp for TaCOBL-5B, and 1,350 bp for TaCOBL-5D. These three genes have been submitted to GenBank under accession numbers MN082634, MN082635, and MN082636, respectively.
Prediction of functional conserved domains of TaCOBLs confirmed that they belonged to the COBRA protein family (Supplementary Fig. S2). In addition, multiple sequence alignments showed that all members had the CCVS conserved domain and the ω-sites, which was consistent with the conserved structure of COBRA protein (Supplementary Fig. S3). Based on protein structure, TaCOBLs showed 92% similarity with OsBC1 [23], 93% with ZmBK2 [24], 91% with AtCOBL4 [20], 89% with SlCOBRA (GenBank accession number: JN398667) [36], and 96% with BoCOBL (GenBank accession number: EU247930) [37] (Supplementary Fig. S3), implying they are orthology at the level of protein structure. The overall protein sequences including COBLs in the eudicots and monocots were utilized to construct an unrooted tree to clarify the phylogenetic relationship of COBL genes from different species. The COBL family members were clustered into two subgroups Group I and Group II (Fig. 2). The phylogenetic analysis revealed that the COBL genes in monocotyledons (TaCOBL in wheat, OsBC1 in rice, ZmBK2 in maize, and BoCOBL in bamboo) and those in dicotyledons (SlCOBRA in tomato and AtCOBL4 in Arabidopsis) were distinguished as separate sub-clades. String protein interactive network database was used to predict the the functional partners of TaCOBL. There were ten proteins may interact with TaCOBL protein (Supplementary Fig. S4). The protein partners Traes_3DL_940A89913.1, Traes_3B_29EA1B8AF1.1, Traes_3DL_B2FD2FBFA.1, Traes_1AL_F420A1BBE.1 and Traes_3B_53FF7A343.2 were mainly involved in cellulose biosynthesis. The five other partners were uncharacterized proteins.

Stress-induced expression of TaCOBL in wheat
The expression pattern of TaCOBL was analyzed in wheat individuals subjected to four abiotic stresses using quantitative real-time PCR (qRT-PCR). As shown in Fig. 3a, a significant increasing trend occurred with a peak seen at 6 h after cold treatment, followed by a steady decrease in leaves. By contrast, the expression level of TaCOBL did not change significantly in roots. However, upregulation was found in neither leaves nor roots after drought, salt, or ABA treatments (Fig. 3b-d). Conversely, a gradual decreasing trend was observed in TaCOBL transcripts short after the treatments started in both leaves and roots, indicating that TaCOBL-5B might play an important role in response to cold stress.

Genotypic analysis of TaCOBLs
Comparison of the full-length sequences of TaCOBLs among 20 wheat cultivars with different cold-resistance capabilities revealed no sequence divergence between TaCOBL-5A and TaCOBL-5D. However, a 4-bp insertion/deletion (InDel) was detected in the promoter region (− 550 bp) of TaCOBL-5B, forming two haplotypes: Hap5B-a and Hap5Bb (Fig. 4a). The regulatory element analysis showed that this InDel may consist of a classic TATA-box. This InDel provided an opportunity to develop an amplified polymorphism sequence marker to distinguish the TaCOBL-5B alleles. The genomic-specific primer set TaCOBL-3 was developed to screen the wheat natural population (Fig. 4b).
We analyzed the distribution of these two haplotypes across the 336 wheat accessions with different coldresistance capabilities. The results showed that Hap5B-a was more frequent in wheat accessions with stronger cold resistance, whereas Hap5B-b was more frequent in those with weaker cold resistance ( Supplementary Fig. S5). Significant associations between TaCOBL-5B haplotypes and cold resistance were observed in 4 environments. In particular, Hap5B-a was significantly positively correlated with cold resistance, and the correlation coefficients were 0.188, 0.205, 0.139, and 0.123 in E1, E2, E3, and E4, respectively, which were statistically significant (P < 0.05) ( Table 1). Thus, TaCOBL-5B Hap5B-a may represent a favorable haplotype that could make a positive contribution to cold resistance.

Distribution of TaCOBL-5B haplotypes
The wheat production area in China has been divided into ten agro-ecological zones based on environmental characteristics, variety types, and growing seasons [38]. The geographic distribution of the TaCOBL-5B haplotypes were evaluated in the 336 wheat accessions covering five out of ten agro-ecological zones, including I (northern winter wheat region), II (Yellow and Huai River valley winter wheat region), III (low and middle Yangtze River valley winter wheat region), IV (southwestern winter wheat region), and V (Introduced wheat variety of foreign). The distribution frequency of Hap5B-a was in an order of I (42%) > V (38%) > II (30%) > III (19%) > IV (5%) (Supplementary Fig.   S6). In China, Hap5B-a is mainly distributed in northern winter wheat region and Yellow and Huai River valley winter wheat region, indicating that wheat accessions harboring Hap5B-a are more adaptable to lower-temperature areas in China.

Promoter activity analysis of TaCOBL-5B haplotypes
To analyze the promoter activity of TaCOBL-5B haplotypes (Hap5B-a and Hap5B-b) in response to cold stress, a 1608bp fragment was fused to the upstream of GUS gene replacing the original 35S in the pCAMBIA1391 expression vector ( Supplementary Fig. S7). A. tumefaciens strain EHA105 was used as intermedia to express the construct in tobacco leaf cells using the infiltration method. We performed the quantitative analysis of the GUS activity produced by promoters of the two haplotypes using leaf samples of the transgenic tobacco plants. The results showed that the Hap5B-a promoter had stronger activity than the Hap5B-b promoter (Fig. 5a). To confirm the results of the GUS activity analysis, the leaf samples were subjected to histochemical staining. Strong GUS staining was observed in the leaf samples that were transformed with the Hap5B-a promoter, whereas light staining was observed in the leaf samples that were transformed with the Hap5B-b promoter (Fig. 5b), consistent with the GUS activity analysis results.

Discussion
The wheat whole-genome sequencing has already been completed [39][40][41][42][43], which provides a high-quality reference genome sequence for the cloning of wheat genes. In this study, we obtained partial sequences of the TaCOBL genes by homology cloning. Then, using the reference sequence of the wheat genome, we successfully obtained full-length sequences of three homoeologous TaCOBLs from the wheat chromosomes 5A, 5B, and 5D, respectively. TaCOBLs encode a putative GPI-anchored COBL protein. Arabidopsis AtCOBL4, rice OsBC1, and maize ZmBK2 specifically regulate cell wall biosynthesis and determine the orientation of cell expansion [23]. The phylogenetic analysis indicated that TaCOBLs and homologs from other plants formed a gene superfamily. In the phylogenetic tree, TaCOBLs and their orthologous genes OsBC1, ZmBK2, AtCOBL4, SlCOBRA, and BoCOBL were clustered into the same clade (Fig. 2). However, the COBL genes in monocotyledons (TaCOBLs in wheat, OsBC1 in rice, ZmBK2 in maize, and BoCOBL in bamboo) and those in dicotyledons (SlCOBRA in tomato and AtCOBL4 in Arabidopsis) were distinguished in the phylogenetic tree as separate sub-clades. This implied that the functions of homologous genes across genomes and species are generally conserved, but functional specificities may exist. As a GPI-anchored protein located at the cell wall-plasma membrane interface, COBRA follows a typical GPI secretion pathway and can sense information from the cell wall and then transmit it to the cell membrane [32]. It is of high possibility that the COBRA family members play an important regulatory role in response to stresses and hormones. The results of BLAST searches supported this hypothesis, in which the majority of TaCOBL ESTs could response to various abiotic stresses (Table S3). Expression pattern analyses suggested that TaCOBLs could response Fig. 3 Expression pattern of TaCOBL-5B in response to cold, drought, NaCl, and ABA treatments to cold stresses. Many studies have shown that String protein interaction network has high credibility in the application of gene function [44,45]. In this study, the protein interactive network showed TaCOBL could interacted with many genes involved in cellulose biosynthesis. Taking these aspects into consideration, we deduced that wheat TaCOBLs could play a role in cold response by regulating cellulose biosynthesis.
Correlation analyses indicated that there was a close positive relationship between Hap5B-a and cold resistance. Hap5B-a was more frequent in wheat accessions with stronger cold resistance, whereas Hap5B-b was more frequent in wheat accessions with weaker cold resistance ( Supplementary Fig. S5). In addition, the favored Hap5B-a was more frequent in lower-temperature ecological zones of wheat production in China (Supplementary Fig. S6). The above results proved that Hap5B-a was a superior allele in response to cold stress of wheat.
Promoter is a specific DNA sequence regulating the spatial and temporal expression of genes of interest through providing binding sites for other regulatory factors [46,47]. Many studies have reported sequence diversity in gene promoter regions affecting traits in plants. In barley, the insertion of a 1-kb sequence in the upstream region of HvAACT1 enabled barley to adapt to acidic soils [48]. A 23-bp motif found in the upstream regulatory region of the MYB10 gene  could be linked to color phenotypes in apple [49]. In this study, we found a TATA insertion in the 5' flanking region of TaCOBL-5B, viz Hap5B-a. Compared to the Hap5B-b haplotype, Hap5B-a conferred stronger cold resistance in wheat. Further analyses implied the insertion functioning as a classic TATA-box that may contribute to the higher transcript level of TaCOBL-5B. The transient expression of GUS in tobacco leaves under the promoters of the two haplotypes further confirmed the deduction (Fig. 5). Functional markers, developed from polymorphic sites within genes that causally affect phenotypic variation, are ideal tools for marker-assisted selection in wheat breeding [50]. In the present study, a co-dominant functional marker Ta5BF/Ta5BR associated with wheat cold resistance was developed and validated. The marker was shown to be highly relevant to cold resistance and thus can be used in wheat breeding programs aiming at improving the cold resistance of wheat. Nevertheless, some wheat cultivars with a 4-bp TATA-box showed weak cold resistance, while the others without the TATA-box possessed strong cold resistance. Considering wheat cold resistance is a very complex network that is regulated by some minor genes, this may be partly attributed to other genes responsible for wheat cold resistance, especially the Fr (freezing tolerance) and VRN (vernalization requirement) loci [13,14]. In common wheat, the most significant Fr loci have been localized on the long arms of 5A, 5B, and 5D chromosomes that carry Fr and VRN genes. Coincidentally, the three homoeologous TaCOBLs were also found on homologous group 5 chromosomes. However, whether

Conclusions
Marker-assisted selection (MAS) provides a strategy for accelerating the process of wheat breeding. Cold resistance is a very complex network that is regulated by some minor genes. Allele mining and functional marker development for these genes could be useful in breeding for improved cold tolerance. In this study, we report the isolation and characterization of TaCOBL. Expression pattern analyses suggested that TaCOBLs could response to cold stresses. TaCOBL-5B carried two haplotypes Hap5B-a and Hap5Bb. Association analysis indicated that Hap5B-a was positively associated with cold resistance. The co-dominant marker Ta5BF/Ta5BR, which was developed based on variants in the promoter region of TaCOBL-5B, could be used as functional marker for the selection of superior allelic variation in cold-resistance related traits. To improve the cold-resistance related traits in breeding, it will be better to perform marker-assisted backcrossing using a variety with the superior haplotype (e.g., Hap5B-a) as a donor parent.