Genome-wide identification of the PEBP genes in pears and the putative role of PbFT in flower bud differentiation

Although Phosphatidylethanolamine-binding protein (PEBP) genes have been identified in several plants, little is known about PEBP genes in pears. In this study, a total of 24 PEBP genes were identified, in which 10, 5 and 9 were from Pyrus bretschneideri genome, Pyrus communis genome and Pyrus betuleafolia genome, respectively. Subsequently, gene structure, phylogenetic relationship, chromosomal localization, promoter regions, collinearity and expression were determined with these PEBP genes. It was found that only PbFT from PEBP genes of P. bretschneideri was relatively highly expressed in leaves during flower bud differentiation. Whereas, expression patterns of TFL1 homologues, gene23124 and gene16540, were different from PbFT in buds. The expression pattern and the treatment of reduction day-length indicated that the expression of PbFT in leaves were regulated by day-length and circadian clock. Additionally, the phenotype of transgenic Arabidopsis suggested that PbFT played a role in not only promoting flower bud differentiation, but also regulating the balance between vegetative and reproductive growth. These results may provide important information for further understanding of the evolution and function of PEBP genes in pears.


INTRODUCTION
Phosphatidylethanolamine-binding proteins (PEBPs) are named for their evolutionarily conserved phosphatidylethanolamine-binding domains and widely found in plants, animals and yeast (Chautard et al., 2004;Leeggangers et al., 2018;Sun et al., 2018). In plants, the PEBP genes family is mainly divided into three subfamilies, FLOWERING LOCUS T (FT)-like, TERMINAL FLOWER 1 (TFL1)-like and MOTHER OF FT AND TFL1 (MFT)-like (Li et al., 2015). MFT-like is the ancestor of the FT-like and TFL1-like subfamilies and is not found in moss and lycopodium. FT/TFL1 homologous genes appeared along with the evolution of seed plants. FT/TFL1-like gene plays an important role in the transformation process of seed plants from vegetative to reproductive growth (Karlgren et al., 2011).

Protein properties and sequence analysis
The pear PEBP sequences were uploaded to ExPASy (http://web.expasy.org/protparam/) to calculate the number of amino acids, molecular weights and isoelectric points (pI).

Phylogenetic tree construction, gene structure and protein motif analysis
The phylogenetic tree was constructed for the sequences of PEBP, using MEGA v6.0 (Tamura et al., 2013) by the Neighbor-Joining (NJ) method with 1,000 bootstrap replicates. The PEBP gene structures were obtained by alignment of open reading frames (ORFs) with corresponding genomic sequences along with the Gene Structure Display Server (Hu et al., 2014). Sequence logos of domain alignments were created using the MEME (http://meme-suite.org/tools/meme) with the default settings, except that the minimum and maximum motif widths were set to 10 and 60 amino acids, respectively. PEBP of both pears and Arabidopsis were analyzed.

Physical localization of PEBP genes on chromosomes and synteny analysis
To identify chromosomal location of PEBP genes, the genome annotation data was collected and mapped with TBtools (Chen et al., 2018). Synteny analysis of PEBP genes was performed by using TBtools.

Investigation on Cis-elements in the promoter region
The 2,500 bp genomic sequences upstream at the transcription start site (ATG) were extracted from Pyrus bretschneideri genomic database, which was used to analyze putative cis-elements in the promoter regions of the PEBP genes with PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/).

Plant materials, growth conditions and treatment
The 7-year-old pear trees (Pyrus bretschneideri, cv. Xueqing) were used as the trail samples that were grown in a natural environment in Hebei Agricultural University (N 38 , E 115 ).
To clarify the expression of PEBP genes in pears during the flower bud differentiation, leaf and apical bud samples at different developmental stages were collected at 30, 45, 60, 75, 90 and 105 days after full blooming (DAB). Meanwhile, the day-length data was obtained from the internet (https://richurimo.51240.com/). The treatment of reduction day-length on the pear leaves started from 30DAB to 120DAB in the next year. A branch of each pear tree was covered with an opaque bag from 6:30 PM until night (Fig. S1). The leaf samples were collected at 45, 75 and 105DAB in this experiment.
In order to verify the circadian rhythm of the genes, the leaf samples in the natural environment were collected from 9:00 AM to the next 6:00 AM, once for every three hours at 45, 75 and 105DAB.
Arabidopsis (Columbia-0) was used in this study. To generate PbFT overexpression transgenic plants, full-length cDNAs without the stop codon for PbFT were cloned into a PBI-121 vector in the sense orientation behind a cauliflower mosaic virus 35S promoter. Phenotypes of the transgenic lines were examined with T4 generation. Arabidopsis plants were cultured under long days (LD, light/dark, 16 h/8 h) in an artificially controlled growth room at 23 C.
All samples that were needed to further analyze were harvested and quickly frozen in liquid nitrogen and then stored at −80 C for further experiments. At least three independent replications were carried out for all assaying samples.

Meristem microscopy
Apical buds of current shoots harvest at 60DAB, 75DAB and 90DAB were fixed quickly in FAA solution (70% ethanol: formalin: acetic acid = 90: 5: 5, v/v) and dehydrated in a graded ethanol series (1 h in each of 70%, 83%, 95%, 100% ethanol and 100% ethanol). Then, samples were soaked in mixtures of ethanol and xylene (1:1 v/v) 1 h and twice in absolute xylene for 1 h. Subsequently, Samples were infiltrated at 37 C with mixtures of xylene and steedman's wax made up in the proportions of 1:1 (v/v) overnight and at 45 C 12 h, followed by two changes of pure wax (2 h each). Specimens were transferred to kraft paper molds filled in fresh wax at 50 C, then cooling to room temperature. Samples were sectioned into 8 mm slice and dewaxed before sealing. The slices were observed under a microscope and photographed.
RNA extraction, cDNA synthesis and quantitative real-time PCR Total RNA was extracted according to the manufacturer instructions of E.Z.N.A. plant RNA extraction kit (Omega Bio-tek, Norcross, GA, USA) including the gDNase. The first strand cDNA was prepared using FastQuant RT Kit (Tiangen, Beijing, China) according to the manufacturer's instructions together with gDNase.
Leaf samples harvested at 30DAB, 45DAB, 60DAB, 75DAB and 90DAB were using for transcriptome sequencing. Three independent replications were carried out for all assaying samples (since samples harvested at 30DAB and 90DAB have been confirmed as not critical samples in our research (unsubmitted), both were not set biological replicates for transcriptome sequencing). The subsequent library construction and RNA-sequencing (RNA-seq) were completed by the Shanghai Majorbio Bio-pharm Technology Co., Ltd. (Shanghai, China). The data was analyzed on the free online platform of Majorbio Cloud Platform (www.majorbio.com). The RNA-Seq data from this article can be found in the Sequence Read Archive (SRA) under accession number SRR10997902-SRR10997912.
The quantitative real-time PCR (qRT-PCR) reactions were executed using the TransStart Ò Top Green qPCR SuperMix (TransGen Biotech, Beijing, China). The qRT-PCR was performed using Eppendorf Mastercycler Ò ep realplex (Eppendorf, Hamburg, Germany). The thermal cycling conditions were 30 s at 94 C; 44 cycles for 5 s at 94 C, 15 s at 56 C and 10 s at 72 C and then performed the melting curve and cooling program. The data was calculated using the 2 −ΔΔCt (Shi, Zhang & Chen, 2019) method and the PbActin gene of P. bretschneideri was used as an internal control to normalize the expression levels of the target genes. Specific primers of the qRT-PCR were presented in Table S1. Three biological replications were used for the experiment.

Statistical analysis
Statistical Product and Service Solutions v. 17.0 (SPSS, Chicago, IL, USA) was used to analyze the experimental data. T-test and ANOVA Tukey's multiple comparison tests were used to exam significant differences (p < 0.05).

Genome-wide identification of PEBP genes
Based on a genome-wide analysis, 24 PEBP genes were identified in pears, of which 10 PEBP genes were from Pyrus bretschneideri genome, 5 PEBP genes were from P. communis genome and 9 PEBP genes were from the Pyrus betuleafolia genome. The detailed information of PEBP genes is given in Table 1. The PEBP proteins varied more or less with Mapping PEBP genes to the P. bretschneideri genome indicated that 10 PEBP genes were unevenly distributed on eight of the seventeen chromosomes, with six genes on Chr2, Chr3, Chr5, Chr7, Chr11 and Chr12, as well as two on Chr6 or Chr14. Mapping PEBP genes to Pyrus communis genome indicated that 5 PEBP genes were unevenly distributed on five of the seventeen chromosomes, that is, Chr1, Chr3, Chr4, Chr7, and Chr12.

Gene structure and conservative motifs analysis
Gene structure analysis showed that all PEBP genes from pears contained four exons and three introns. The sizes of the second exon and the third exon were exactly 62 bp and 41 bp

Selective pressure analysis
The gene pairs were obtained by synteny analysis. The ratio of non-synonymous substitutions (Ka)/synonymous substitutions (Ks) was evaluated to detect the modes of selection of PEBP genes. Generally, Ka/Ks > 1 indicated a positive selection with an accelerated evolution; Ka/Ks < 1 indicated functional constraints with a negative selection; and Ka/Ks = 1 signified a neutral selection. The results showed that the Ka/Ks ratios of 13 duplicated PEBP gene pairs ranged from 0.11 to 0.85 with an average of 0.27 (Table 2). Additionally, the majority of the ratios were less than 0.39, suggesting that the duplicated PEBP genes of pears had mainly undergone strong purifying selection, with the slowly evolving at the protein level of this gene family.

Cis-elements in the PEBP genes promoters of Pyrus bretschneideri
We analyzed 2,500 bp genomic sequences upstream at the transcription start site (ATG) of Pyrus bretschneideri PEBP genes. About 50 putative cis-elements were identified in each promoter of PEBP genes that were classified into three groups, including hormone responses, plant growth and development and biotic/abiotic stress responses (Fig. 3). There were numerous light-responsive elements on these promoters, such as AE-box, Box 4, G-box, I-box, and GT1-motif. The promoter of gene12374 possessed three Auxin-responsive elements and one GA-responsive element, and the promoter of gene23124 had one Auxin-responsive element, which was not found in the promoters of other Pyrus bretschneideri PEBP genes.

Expression of PEBP genes during the flower bud differentiation
Expression of pear PEBP genes obtained from the RNA-seq data of "Xueqing" leaves at different stages was analyzed (Table S2). It was found that just one gene (gene12374) was expressed in all five development stages and was relatively highly expressed in leaves (Fig. 4A). The other genes had low or no expression. Phylogenetic analysis revealed that the closest relationship between gene12374 and Arabidopsis PEBP gene was FT, so gene12374 was re-named PbFT. PbFT expression patterns with the qRT-PCR experiment were shown in Fig. 4B. The expression trends of PbFT were similar to RNA-seq data. The highest expression of PbFT was consistently observed at 75DAB. Then, the day-length of a period from 30DAB to 120DAB was recorded (Fig. 4C). According to the microscopic observation of current shoot apical meristems, it was not found abnormal at 60DAB. This period was maybe in the flower bud physiological differentiation stage. Whereas, it was found that the flower bud morphological differentiation stage had already entered at 75DAB (Figs. 4D-4F). Meanwhile, the expressions of PbFT was also the highest in apical buds (Fig. 5A). Those findings suggested that the period from 60DAB to 75DAB was critical for the flower bud differentiation, and the expression of PbFT was increasing strongly during this period. Additionally, the expressions of two TFL1 homologues (gene23124 and gene16540) were higher before 60DAB in apical buds. The lowest expressions of gene23124 and gene16540 occurred at 75DAB. Subsequently, the expression of gene23124 remained at a lower level, but the expression of gene16540 increased gradually (Figs. 5B and 5C). Interestingly, the time of the highest expression of PbFT was found around the day with the longest daylength. Subsequently, the treatment of day-length reduction on the pear leaves showed that if the day-length was reduced, the expression of PbFT also decreased (Fig. 6). The expression of PbFT under the natural condition was up to twice as much as under the treatment of day-length reduction at 75DAB.
On the other hand, the expression of PbFT was at a high level during the day, but it dropped suddenly as night fell and remained at a lower level until dawn (Fig. 7). Overexpression of PbFT in relation to acceleration of flowering time To investigate the role of PbFT in flower bud induction regulation, PbFT was overexpressed under the control of the 35S promoter in the wild-type Arabidopsis (Col-0). A total of nine independent lines of overexpressing PbFT (PbFT-OE) plants that showed obviously enhanced expression of PbFT were obtained. Two independent lines were randomly selected for further analysis. The transgenic lines were confirmed by PCR with the 35S forward primer (35S-F) and PbFT-specific reverse sequencing primer (PbFT-R) ( Table S1). The days from germination to flowering was 32.11 for WT plants, but 24.33 and 25.33 with PbFT-OE3 and PbFT-OE9 lines. The corresponding total rosette leaf number was 18.29, 8.57 and 10.17 in WT, PbFT-OE3 and PbFT-OE9 lines. The days from germination to flowering were significantly accelerated and the total rosette leaf number decreased significantly in the PbFT-OE lines compared with the WT plants, but two overexpression lines showed no significant difference in the relative parameters (Fig. 8).
The diameters of the transgenic lines PbFT-OE3 and PbFT-OE9 main stem were 0.76 and 0.82 mm, which were significantly larger than 0.5 mm of WT. In the absence of human interference, the transgenic plants were lodging, as opposing to the wild-type upright phenotype (Fig. 8).

Conservation and evolution of PEBP genes in pears
So far, the genome sequencing of "DangshanSuli" (Pyrus bretschneideri), "Bartlett" (Pyrus communis) and "Shanxi Duli" (Pyrus betuleafolia) has been completed. These existing information resources can be used to analyze the evolution and gene function of the PEBP family from the perspective of bioinformatics. The PEBPs were a kind of evolutionarily conserved proteins. As the PEBP genes of other species like Arabidopsis and apple (Malus domestica), the PEBP genes of pears consisted of four exons and three introns and meanwhile the second exon contained 62 bases and the third exon contained 41 bases. In addition, the Ka/Ks of PEBP gene pairs obtained from collinearity analysis in different genera were far less than 1. This result indicated that the PEBP genes of pears have undergone purifying selections, which also confirmed that the PEBP family was quite conservative. However, the number of PEBP genes and their distribution on chromosomes of Pyrus bretschneideri, Pyrus communis and Pyrus betuleafolia were different in this study. Ten and nine genes from Pyrus bretschneideri and Pyrus betuleafolia were identified, while only five genes were identified from Pyrus betuleafolia. This result suggested that different degrees of gene loss had occurred during the evolution of Pyrus, especially in the Pyrus communis genome, at least in the PEBP family of "Bartlett" pear. Interestingly, among the materials used for genome sequencing, "DangshanSuli" (Pyrus bretschneideri) and "Shanxi Duli" (Pyrus betuleafolia) belonged to Asian pears and "Bartlett" (Pyrus communis) was a European pear. This situation might be due to the independent domestication processes for either Asian or European pears .
In addition, more attention should be paid to the genes in the FT-like group. FT and TWIN SISTER OF FT (TSF) had redundant functions in promoting flowering (Jin et al., 2015;Yamaguchi et al., 2005), and they had five ortholog genes in pears. There were two genes belonging to both "Bartlett" and "Shanxi Duli", but just one, PbFT, was identified in "DangshanSuli" genome. Since PbFT was so closely related to the AtFT and AtTSF, we speculated that it may also be similar in functions.

Expression and function of PbFT
FT and TFL1 played important roles in flowering induction (Gao et al., 2017;Putterill & Varkonyi-Gasic, 2016). The expressions of PbFT and TFL1 homologues (gene23124 and gene16540) in pear buds showed the opposite expression patterns before and in the critical period of the flower bud differentiation. However, it was found that only PbFT among pear PEBP genes was relatively highly expressed in leaves. It suggested that PbFT played positive roles both in leaves and buds, but gene23124 and gene16540 played negative roles only in buds. In Arabidopsis, TFL1 expression was merely limited to shoots (Baumann et al., 2015). This was similar to our results. It suggested that gene23124 and gene16540 had the opposite function to PbFT. As a receptor of the plant, the leaves can directly sense the changes of the external environment, especially the time of day-length, so as to regulate the growth and development of the plant through corresponding signal transduction. FT protein is synthesized in the leaf vasculature and transported to the apical meristem of a shoot through the phloem (Corbesier et al., 2007;Lin et al., 2007;Tamaki et al., 2007;Zeevaart, 2008). By analyzing the promoter, it was found that there were typical light-responsive elements, AE-box, Box 4, G-box and TCT-motif, on the PbFT promoter. This suggested that the expression of PbFT may be regulated by light. According to our transcriptome data, only the PbFT gene, among the PEBP genes of Pyrus bretschneideri, had a relatively high expression level in the leaves. It was found that the expression of PbFT was increasing strongly during critical period, which was from the flower bud physiological differentiation stage to morphological differentiation stage and the highest expression of PbFT was observed at 75DAB, when was about the longest daytime of a year. Shortening the day-length of some branches indicated that the expression of PbFT was downregulated. This was a good indication that the expression of the PbFT gene was related to the flower bud induction and considerably regulated by the day-length in a year. Many scholars have concluded that FT is strongly expressed around dusk in LD in the incubator (Mouradov, Cremer & Coupland, 2002;Song, Ito & Imaizumi, 2013). However, we found that PbFT gene remained at a high expression level during the day and at a lower level in the evening until dawn. PbFT expression showed an unusual circadian rhythm. As the temperature and light in the incubator are usually constant but greater changes exist in the natural environment, these may result in a different pathway of regulation mechanism with PbFT.
At least six signal pathways (photoperiod, autonomous, age, gibberellin, vernalization, ambient temperature) that regulate the flowering process have been demonstrated (Blumel, Dally & Jung, 2015;Kim et al., 2009;Mutasa-Gottgens & Hedden, 2009;Srikanth & Schmid, 2011;Wang, 2014). As a mobile signal florigen gene, in several pathways, FT can initiate the floral induction in some species (Blumel, Dally & Jung, 2015;Cao et al., 2015;Corbesier et al., 2007;Liu et al., 2013;Oda et al., 2012;Pin & Nilsson, 2012;Tamaki et al., 2007). The research on FT genes has been progressed in many horticultural plants (Wilkie, Sedgley & Olesen, 2008). However, the role of the PbFT during flower bud differentiation remains unclear in pears. In this study, it was found the expression of PbFT increased during the period of floral initiation. Moreover, the overexpression of PbFT dramatically accelerated the starting time for flowering Arabidopsis. This implies that the PbFT might play a role in promoting flower bud differentiation in pears.
Additionally, homologs of FT genes in crops have pleiotropic functions that can mediate numerous developmental processes, such as growth, plant architecture control, fruit set and tuber formation (Pin & Nilsson, 2012). AcFT1 promotes bulb formation, while AcFT4 inhibits bulb production in onions (Lee et al., 2013). Tuberization is regulated by the FT homolog StSP6A through modulation of source-sink in potatoes (Abelenda et al., 2019;Lehretz et al., 2019). Our transgenic plants were thinner than WT, which was caused by overexpression of PbFT. This finding suggests that PbFT has a role in not only promoting flower bud differentiation, but also regulating the balance between vegetative and reproductive growth.

CONCLUSIONS
We identified 10, 5 and 9 PEBP genes from genomes of Pyrus bretschneideri, Pyrus communis and Pyrus betuleafolia, respectively. The PEBPs is a kind of evolutionarily conserved proteins. However, the number of PEBP genes and their distribution on chromosomes of Pyrus bretschneideri, Pyrus communis and Pyrus betuleafolia are different. This situation may be due to the independent domestication processes for either Asian or European pears. Moreover, we determined that PbFT, which can be regulated by day-length and circadian clock, possessed the function of promoting flower bud differentiation in pears. Additionally, PbFT may play a role in regulating the balance between vegetative and reproductive growth.