Comparative genomic analysis of the COBRA genes in six Rosaceae species and expression analysis in Chinese white pear (Pyrus bretschneideri)

Identification of COBRA genes in six Rosaceae species

In this work, the Pyrus bretschneideri genome was downloaded from GIGADB datasets (http://gigadb.org/dataset/100083), and five Rosaceae genomes (Fragaria vesca, Rubus occidentalis, Prunus avium, Malus domestica, Prunus mume) were obtained from the following website (https://www.rosaceae.org/). Arabidopsis COBRA gene family members amino acid sequences were used as the query sequence for BlastP (Protein BLAST: search protein databases using a protein query (nih.gov)) search (E = 0.001) from the local protein database. The SMART online software program (http://smart.embl-heidelberg.de/) was used to screen the genes (Letunic, Doerks & Bork, 2012). The protein sequences lacking a whole COBRA domain and redundant sequences were discarded. We used the ExPASY online website to predict the molecular weights, isoelectric points, and GRAVY values of COBRA protein (http://web.expasy.org/protparam/) (Artimo et al., 2012). Signal peptides were analyzed and predicted by SignalP 4.1 server (https://services.healthtech.dtu.dk/service.php?SignalP-5.0) software.

Evolutionary analysis of COBRA gene family

Sequence alignment of all COBRA proteins was done using the ClustalW tool in MEGA 7.0 software. The phylogenetic tree was constructed with MEGA 7.0 software using the NJ (Neighbor-Joining) (bootstrap = 1,000). The sequences of A. thaliana, O. sativa, Z. mays, and Populus were obtained from the article (Roudier et al., 2002; Li et al., 2003; Sindhu et al., 2007; Zhang et al., 2010).

COBRA gene structures and conserved motif prediction

The Gene Structure Display Server (http://gsds.gao-lab.org/) was used to investigate the gene structures (introns/exons) of the COBRA gene family members (Guo et al., 2007). The conserved sequences of the COBRA gene family were analyzed by MEME online (https://meme-suite.org/meme/) software (Bailey et al., 2015). The parameter settings are as follows: The number of identified motifs were 20 with Parameters for the conserved motif prediction were motif width greater than 6 and less than 200.

Promoter analysis and cis-acting element analysis of PbCOBL genes

In this study, we found the 1,500 bp–2,000 bp upstream of the start codon in the pear genome database, which was the promoter sequence of the gene. The online software PlantCare (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was used to examine cis-acting elements.

Chromosomal location, gene duplication, and Ka/Ks ratio analysis

Mapinspect software was used to map the position of COBRA genes on chromosomes (Zhu et al., 2015; Su et al., 2019a). The determination of gene replication events of COBRA genes in 6 species mainly followed the following principles: First, the similarity of the two genes was greater than 80%. The distance between two genes on the same chromosome was more than 200 kb, which was tandem-duplicated genes. Two genes located on different chromosomes were defined as segmentally duplicated genes. Finally, DnaSP v5.0 software calculated non-synonymous (Ka) and synonymous substitution (Ks) values and performed a sliding window analysis. The parameter was set to window size 150 bp and step size 9 bp (Zhao et al., 2021).

RNA extraction and qRT-PCR analysis

In this study, Dangshan su pear was used as the material, which growed in Dangshan County, Anhui Province. Samples of current-year flowers, buds, stems, mature leaves and fruits were obtained. The fruits were picked 15 days after pollination (DAP), 23 DAP 39 DAP, 47 DAP, 55 DAP, 63 DAP, 79 DAP, and 102 DAP. 39 DAP fruits were selected for different tissue expression analyses. We used the RNA extraction kit of Tiangen (Beijing, China) to extract RNA from different materials. Reverse transcription was performed using a PrimeScriptTM RT reagent kit with gDNA Eraser (TaKaRa, China). The qRT-PCR primers were designed using Beacon Designer 7 software (Table S1). The pear Tubulin gene (No. AB239680.1) was used as an internal reference (Wu et al., 2012). Each 20 µL qRT-PCR system included 10 µL of SYBR Premix Ex Taq TM II, 2 µL cDNA, 0.8 µL of Forward primer and Reverse primer, 6.4 µL of water. This study was conducted according to the procedures in the instruction manual, and three biological repetitions were run for each sample. The relative expression level of genes was calculated by the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Identification, characterization analysis of COBRA genes

Using the amino acid sequences of Arabidopsis COBRA gene family members as probes, we identified 87 COBRA proteins in six Rosaceae species. Including 16 COBRA proteins (PbCOBL1-PbCOBL16) in Pyrus bretschneideri, 22 in Malus domestica (MdCOBL1-MdCOBL22), 13 in Fragaria vesca (FvCOBL1-FvCOBL13), 11 in Prunus mume (PmCOBL1-PmCOBL11), 13 in Rubus occidentalis (RoCOBL1-RoCOBL13) and 12 in Prunus avium (PaCOBL1-PaCOBL12). The detailed information (gene name, gene identifiers, amino acid number, signal peptide, molecular weight, theoretical isoelectric point, Grand average of hydropathicity and subdivided subgroup) of each COBRA was presented in Table 1 and Table S2. These results showed that the largest molecular weight among the six Rosaceae species was MdCOBL12, which was 131.89 kDa. The smallest was PaCOBL4, which was 12.59 kDa. Except for MdCOBL17 and PaCOBL3, the grand average of hydropathicity of other genes were negative, indicating that most members of COBRA genes were hydrophilic proteins. In these Rosaceae species, the lowest pI value was 5.04 (MdCOBL16), whereas the highest pI value was 9.59 (MdCOBL1). Among the six Rosaceae species, about 74% of COBRA proteins contained signal peptides, and 26% of proteins did not contain a signal peptide.

Phylogenetic and hydrophobic analysis of COBRA genes

Phylogenetic analysis showed that the 87 COBRA genes were classified into two subclasses (Group A and Group B), similar to the Arabidopsis (Fig. 1, Table S5). Group A was structurally similar to AtCOBRA, and Group B had a higher similarity to AtCOBL7. We divided Group A into Class1–Class6 and Group B into Class7–Class8. MdCOBL17, RoCOBL7, RoCOBL5, and PbCOBL11 were independent branches in the evolutionary tree, and no genes related to them have been identified. Most PbCOBL genes were more tightly grouped with MdCOBLs. In Class1, AtCOBL4 and GhCOBL9A were clustered with MdCOBL11, MdCOBL13, MdCOBL2, MdCOBL3, PbCOBL13, PbCOBL5, PmCOBL3, FvCOBL5, RoCOBL12. OsBC1, ZmBK2, PtrCOBL4 with PmCOBL4, PaCOBL9 in one branch. AtCOBL2, OsBC1L4, and OsBC1L6 appeared in Class3, and OsBC1L5 was alone in Class8.

We compared the similarity of COBRA proteins of six Rosaceae species (Fig. 2, Figs. S1–S5). In Pyrus bretschneideri, comparisons among the Group B subgroup genes showed identity in the range of 47.11% to 94.36%. Within the Group A subgroup, the proteins are 15.22% to 95.61% identical. Protein similarity between the two subgroups ranged from 9.58% to 47.11% (Fig. 2). In Prunus mume, the similarities among the three members of the Group B subfamily were 53.3% for PmCOBL1-PmCOBL9, 50.66% for PmCOBL1-PmCOBL10, and 49.32% for PmCOBL9-PmCOBL10, respectively. The similarity between Group A and Group B subfamilies ranged from 12.28% to 19.25%, and the similarity between Group A subfamily members ranged from 15.27% to 76.75% (Fig. S1). In Rubus occidentalis, the similarity with other members of RoCOBL7 was low because RoCOBL7 was a separate branch in the evolutionary tree. The similarities between the three members of the Group B subclade, RoCOBL1, RoCOBL3, and RoCOBL6, were 49.78%, 56.63%, and 50.6%, respectively. The protein similarity between the two members of the two subclades ranged from 5.86% to 16.59%. The similarity of Group A subclade members ranged from 7.63%–85.40% (Fig. S2). In Fragaria vesca, Group B, Group A and two subclades individual member similarities ranged from 47.96% to 56.29%, 15.4%–93.21%, and 7.49%–19.45%, respectively (Fig. S3). In Prunus avium, the similarity of three members PaCOBL2, PaCOBL11, and PaCOBL12 was 23.45%, 41.84%, and 28.14%, respectively. The results of the comparison between subgroups were 8.05%–18.39%. Protein similarity among Group A members ranged from 6.64%–72.44% (Fig. S4). In Malus domestica, Group B, Group A and two subclades, individual member similarities ranged from 25.68%–79.05%, 5.98%–94.04%, 5.98–94.04% (Fig. S5).

We studied their hydrophobicity to determine if the 87 proteins were likely to have a GPI anchor similar to COBRA. We performed the hydrophobic analysis of 87 COBRA genes from six Rosaceae species. As shown in Fig. 2 and Figs. S1–S5, most amino acids showed a similar trend, with the middle part hydrophilic and the ends hydrophobic. The GPI modification sites and potential w-cleavage sites of 87 proteins were predicted using big-PI. Among the 87 proteins, the program found 34 significant potentials for GPI modification using the default parameters and proposed a possible GPI addition for the remaining three proteins (FvCOBL1, PbCOBL1, PbCOBL4, PmCOBL5, PmCOBL11, RoCOBL2, PaCOBL1). All proteins have potential w-cleavage sites.

Structural and conserved motif analysis of COBRA proteins

To gain a comprehensive understanding of the diversity of COBRA genes in the six Rosaceae families, we performed gene structure and conserved sequence analysis of 87 genes (Fig. 3). Group A had 66 members with the number of exons ranging from one to 12, of which 30 members had six exons, five members contained 12 exons, and one member contained 13 exons. Group B had 21 members, of which eight members contained only two exons, three members contained three exons, two members contained one exon, and one member contained four exons. In addition, PbCOBL8, PmCOBL9, PbCOBL2 contained six exons. MdCOBL16, PbCOBL15, MdCOBL5 contained seven exons and MdCOBL12 contained 13 exons. We performed a conserved structure analysis of 87 genes using MEME online software. We found that motif 8, motif 6, and motif 10 were relatively conservative, motif 15 and motif 13, and motif 14, were specific to Group B members (except for PaCOBL11), and motif 2 was specific to Group A.

We compared the sequences of six species (Fig. 4, Figs. S6–S11). Four conserved structural domains were identified in 87 genes, N-terminal signal peptide, Carbohydrate-binding module (CBM), central cysteine-rich domain(CCVS), and C-terminal hydrophobic domain, respectively. Sequence comparison revealed that GroupB members generally had more amino acids than Group A, and the N-terminal signal peptide was different in the two subgroups. Analysis of the conserved regions of the six species COBRA members revealed that the central cysteine-rich domain (CCVS) was the relatively conserved region in both subclades, and almost all members contained the CCVS structural domain. The Central cysteine-rich domain contained a consensus N-glycosylation site, this site was mainly associated with post-translational modifications of GPI-anchored proteins and more generally with extracellular proteins. There was also an N-glycosylation site in the C-terminal hydrophobic domain. N-terminal signal peptide and C-terminal hydrophobic domain played an important role in function. However, the sequence comparison results showed that the similarity between N-terminal signal peptide and C-terminal hydrophobic domain in each subgroup was not high, suggesting that there was no solid selective pressure on these areas as long as their hydrophobic nature was considered.

Chromosomal location and duplication events of COBRA family genes in six Rosaceae

Based on the genome-wide data of pear, strawberry, black raspberry, sweet cherry, japanese apricot, and apple, all COBRA genes exact chromosomal physical localization was determined, as shown in Fig. 5. In pear, 11 PbCOBL genes were distributed on six chromosomes (Chr3, Chr6, Chr8, Chr13, Chr14, Chr17), and five genes were not localized on any chromosome. In strawberry, 13 genes were located on five chromosomes (except Chr2 and Chr7). In Japanese apricot, four genes were distributed on chromosomes 7; two genes on chromosomes 2 and 4; and one gene on chromosomes 3 and 6, respectively. In black raspberry, there were no genes on chromosome 1; three genes on chromosomes 3, 5, and 6, and one gene on chromosomes 2, 4 and 7. In apple, 18 genes were distributed on chromosomes 3, 6, 8, 9, 10, 11, 13, 14, 16 and 17 (except Chr1, 2, 4, 5, 7, 12 and 15). In sweet cherry, 12 genes were distributed on six chromosomes (except Chr6 and 7), of which there were four genes on chromosome 3, only one gene on chromosome 1, 2 and 8, two genes on chromosome 4 and three genes on chromosome 5.

To understand the role of drivers of gene duplication in the evolutionary process, we calculated Ka, Ks, and Ka/Ks ratios of duplicated gene pairs in six Rosaceae species (Fig. 5, Table S3). Ka/Ks = 1 is the cut-off value that indicates neutral selection, Ka/Ks < 1 represents negative selection and Ka/Ks > 1 represents a positive selection. All 14 duplicated gene pairs were identified in pear, strawberry, apple, Japanese aprico, and black raspberry. There were four duplicated gene pairs in pear, two pairs in strawberry, six pairs in apple, one pair in plum and 1 pair in black raspberry. Among the 14 replicated gene pairs, only MdCOBL3-MdCOBL13 and RoCOBL10-RoCOBL11 Ka/Ks > 1 were 1.272 and 1.141, respectively, and the remaining replicated gene pairs had Ka/Ks < 1. Among the 14 replication gene pairs, 12 gene pairs experienced segmental duplication, and only two pairs of genes (RoCOBL10-RoCOBL11, FvCOBL1-FvCOBL3) experienced tandem duplication. To gain a comprehensive understanding of the selection pressure on COBRA genes, we performed a sliding window analysis (Fig. 6), which showed that in 14 pairs of gene replication events, most of the Ka/Ks loci were less than 1, and only very few loci had Ka/Ks values greater than 1.

Analysis of cis-acting elements of COBRA gene family promoter in Pyrus bretschneideri

In order to study the specific expression of the COBRA gene, We used the plant care database to study cis-acting elements of promoters of 16 COBRA genes in Pyrus bretschneideri (Fig. 7, Table S4). The promoter of COBRA genes contained many cis-acting elements related to hormones, including responses to Auxin-responsive element (TGA-box, AuxRE), MeJA-responsiveness (CGTCA-motif, TGACG-motif), Salicylic acid responsiveness (TCA-element), Gibberellin-responsiveness, Abscisic acid responsiveness. CGTCA-motif and GACG-motif cis-acting elements were identified in all 16 genes. A total of 11 cis-acting elements associated with auxin were found in 11 members. Among them, AuxRE was only present in PbCOBL14, TGA-box was not contained in PbCOBL1,3,4,7,13,14, and all other genes were contained. TCA-element, which are cis-acting elements involved in salicylic acid responsiveness, were identified in PbCOBL1, 2, 3, 5, 6, 8, 10, 11, 12, 13, 15 and 16. Gibberellin-responsiveness (P-box, GARE-motif, TATC-box) was identified only in PbCOBL4,9 and 13. Abscisic acid responsiveness (ABRE) were more prominentin numbers, 43 ABRE cis-acting elements were identified in 13 genes, and only PbCOBL3,7,9 was not identified with ABRE.

It is related to plant growth and development, including participating in light response elements (G-box, GT1 motif, GATA motif, Sp1, Box 4, TCT-motif, AE-box, Lamp element, CHS-cma1a, GA motif, I-box, 3-af1 binding site, ACE, TCCC motif) were identified 154 times in 16 genes. The cis-acting regulatory element O2-site involved in the regulation of maize alcohol-soluble protein metabolism was identified in 10 members (PbCOBL2,4,5,6,8,10,11,12,15 and 16). Meristem expression (CAT-box) was identified in PbCOBL3,4,7 and14. Endosperm specific negative expression (AACA motif). Wound-responsive element (WUN-motif) and Endosperm expression (GCN4-motif) were only identified in PbCOBL9. cis-acting elements associated with biotic and abiotic stress responses were also identified in COBRA genes. Six PbCOBL genes (PbCOBL1, 3, 4, 7, 13, and 14) had MBS, and two members (PbCOBL3,14) contained TC-rich repeats. Anaerobic induction (ARE) was identified in 15 genes (Except PbCOBL7). A total of 32 identifications in 13 members of the low-temperature responsiveness (LTR). COBRA family promoters were mostly engaged in hormone response and light response processes. Plant growth and development are facilitated by the reaction to various hormones and light responses.

Expression characteristics of Chinese white pear COBRA genes

To gain a deeper understanding of the COBRA gene family, we performed expression pattern analysis in the stems, leaves, fruits, flowers and buds of Dangshan su pear (Fig. 8). The results showed that PbCOBL1,3,10,12,13,14 were highly expressed only in fruits and almost not in other tissues. PbCOBL2 was highly expressed in flowers. PbCOBL5,15 were expressed at relatively high levels in several tissues. PbCOBL9 was highly expressed in stems. PbCOBL4,6,7,8,16 were expressed at low levels in leaves and at high levels in other tissues.

We examined the expression pattern of the COBRA gene family in eight periods of fruit development of Dangshan su pear (Fig. 9). As shown in the figure, PbCOBL1 and PbCOBL4 were substantially expressed at 15 and 23 DAP and non significant expression were found on remaining time interval. PbCOBL6,8,11,12,16,13 showed similar expression patterns, with higher expression in early fruit development (15 DAP, 23 DAP, 39 DAP) and a decreasing trend in late fruit development. PbCOBL10 expression showed two peaks at 23 DAP and 47 DAP during eight periods of fruit development, and PbCOBL2 showed a similar expression pattern with PbCOBL10, with two peaks at 39 DAP and 102 DAP. PbCOBL7 was barely expressed in 79 DAP. PbCOBL5 and PbCOBL9 showed a peak at 55 DAP, with lower expression in other periods. PbCOBL14,15 were expressed at each period of fruit development, with PbCOBL14 highly expressed at early fruit development (15,23,39 DAP) and 63 DAP, PbCOBL15 peaking at 47 DAP, and low expression at other periods. PbCOBL3 was mainly highly expressed in 39 DAP and 79 DAP.