Genome-wide identification, phylogeny and expression analysis of the PME and PMEI gene families in maize

Pectins, the major components of cell walls in plants, are synthesized and secreted to cell walls as highly methyl-esterified polymers and then demethyl-esterified by pectin methylesterases (PMEs). The PMEs are spatially regulated by pectin methylesterase inhibitors (PMEIs). In this study, 43 and 49 putative PME and PMEI genes were identified in maize, respectively. Gene structure and motif analysis revealed that members in the same paralogous pairs or in the same subgroup generally had common motif compositions and gene structure patterns, which indicates functional similarity between the closely related ZmPME/PMEI genes. Gene ontology annotation analysis showed that most of the ZmPME/PMEI genes are involved in cell wall modification and pectin catabolic process with molecular functions of pectinesterase or pectinesterase inhibitor activities. There are 35 ZmPME/PMEI genes expressed higher in anthers than in other tissues from the NimbleGen maize microarray data, and the semiq-RT-PCR assay revealed most of these ZmPME/PMEIs specially expressed in anthers and pollens, indicating they possibly had role in anther and pollen development. In addition, these ZmPME/PMEI genes were highly expressed in the fertile anthers, while lowly or no expressed in sterile anthers. This further indicated these genes might be involved in the development of anther and pollen.

For example, OsPMEI28 overexpression in rice had an effect on the growth process, which resulted in a dwarfed phenotype 31 , and overexpression of the PMEI5 resulted in a higher demethylesterification of seeds and reduced the PME activity, which was accompanied by an earlier and faster germination process compared to wildtype in Arabidopsis 32 .
In recent years, many reports have shown that some PME/PMEIs regulate plant stress resistance and pollen development. AtPMEI10, AtPMEI11 and AtPMEI12 were identified as upregulated in response to B. cinerea infection 33 . Expression profile of the genes TaPME21-2, TaPME21-1/2/4, TaPME58, TaPME63 and TaPME67 was induced in the susceptible cv. Bobwhite and repressed in the resistant cv. Sumai 3 34 . The transgenic rice overexpressing OsPME14 showed higher PME activity and Al content in root tip cell wall, and became more sensitive to Al stress 9 . In flax, 48 (77.4%) PME genes and 53 (80.3%) PMEI genes had higher expression level in the flowers 10 . In Arabidopsis, 15 PMEs were highly expressed in pollen and 10 of these contained PRO regions 35 . These suggest that the PME/PMEIs might play important roles in pollen development. Mutations of VANGUARD1 (VGD1), the type I PME gene with the highest expression levels in Arabidopsis pollen tubes, resulted in retarded growth in the style and transmitting tract and subsequent reduction in male fertility 36 . In maize, the ZmC5 of PMEs has a role in pollen tube elongation 37 and ZmGa1P, a pollen-expressed PME gene, can confer the male function in the maize unilateral cross-incompatibility (UCI) system 38 .
In this study, genome-wide identification of ZmPME/PMEI genes was firstly conducted in maize, and the phylogenetic tree, gene structure, conservative motif, expression, gene ontology annotations were also examined. In addition, semiq-RT-PCR assay was conducted to verify the gene expression pattern of some ZmPME/PMEI genes highly expressed in anthers, since more than half the genes highly expressed in anthers according to the NimbleGen maize microarray data. To further evaluate their possible roles on pollen development, gene expression of some ZmPME/PMEI genes in fertile and sterile anthers was also investigated. Results in this study would provide useful information for further investigate the function of maize PME/PMEIs, especially on the development of anther and pollen.

Identification of ZmPME/PMEI genes in maize.
To identify putative ZmPME/PMEI genes in maize genome, we searched the maize genome annotation data with known plant PME/PMEI domains (pfam01095/ pfam04043) as a query using HMMER 3.0 package 39 . In total, we obtained 43 putative PME genes and 49 putative PMEI genes in maize. These genes were designated as ZmPME1-43 and ZmPMEI1-49 ( Fig. 1 and Supplementary  Table 1), of them, 20 genes (PME1-20) had PRO region (which showed similarities with the PMEI domain) and the PME domain. Each ZmPME/PMEI gene model was selected by analyzing the similarity between the ZmPME/PMEI genes and homologous genes, as most of the ZmPME/PMEI genes had more than one transcript in the MaizeGDB database (https://www.maizegdb.org/). Then we randomly selected 15 ZmPME/PMEIs for reverse transcription polymerase chain reaction (RT-PCR) to assess the veracity of the ZmPME/PMEI genes models. The results indicated that the 15 ZmPME/PMEI genes were expressed in maize pollen and only a single amplicon was found ( Supplementary Fig. S1). The most identified ZmPME/PMEI genes encode proteins with 150-250 amino acids (aa). They are ranging from 149 (ZmPME23) to 1,360 (ZmPME28) aa, with an average of 346 aa ( Supplementary Fig. S4), and their isoelectric points (pI) are 4.28 to 10.23. These ZmPME/PMEIs are distributed on all the 10 maize chromosomes, and chromosomes 1, 2, 3, 7 and 8 have more ZmPME/PMEIs than others ( Supplementary Fig. S2).
Phylogenetic analysis. Phylogenetic trees were constructed by using MEGA 7.0 with the neighbor-joining model. In order to analyze the evolutionary relationships among the predicted ZmPMEs and ZmPMEIs, we aligned maize acid sequences with 101 and 106 predicted PMEs and PMEIs from rice and Arabidopsis. On the basis of phylogeny, the PMEs and PMEIs families in plants were subdivided into 5 and 12 groups, respectively ( Supplementary Fig. S3). PMEs in each group and PMEIs in groups I to III, V, VI and VIII are all from the three species ( Supplementary Fig. S3), indicating that these ZmPME/PMEIs might have the conserved function in evolution.
Meanwhile, according to cluster analysis, the ZmPME/PMEI families could be divided into 5 and 8 subfamilies, respectively (Fig. 1). The PMEI domains may be derived from duplication and divergence of the PRO domain and have rapidly evolved 12 . So, we constructed the PMEI phylogenetic tree used the protein sequences of all the ZmPMEIs and the 20 ZmPMEs containing PRO region. The ZmPMEI subfamily I includes the same genes in the ZmPME subfamily I (except ZmPME21 and ZmPME22), both of the subfamilies are the largest subfamily, and the genes had signal prediction or transmembrane region domain. For ZmPMEs, genes in the subfamilies II and III do not have signal prediction or transmembrane region domain (except ZmPME27); while genes in the subfamily IV have signal prediction domain (except ZmPME32 and -35); and the pI of subfamily V are higher than 8. For ZmPMEIs, genes in the subfamilies II and VII have signal prediction domain, and their pI are higher than 7 (except ZmPMEI26, -37 and -41); most genes in the subfamilies III, IV and VIII expressed higher in anthers than in other tissues (Fig. 2). Homologous ZmPME/PMEI genes were identified 24 paralogous pairs in maize (Supplementary Table 2). The value of the nonsynonymous substitution rate (Ka) to the synonymous substitution rate (Ks) substitutions (Ka/Ks) can be used as an indicator which could reflect selection pressure of a gene or a gene region during evolution. To infer the influence of selection on the evolution of the maize, we estimated Ka/Ks values for all of them (Supplementary Table 2). The Ka/Ks values of all the homologous genes are between 0.0033 and 0.3889, suggesting that most of the ZmPME/PMEI genes undergone negative selection and evolved slowly. The Ka/Ks values of maize PMEs paralogs are significantly lower than that of the PMEIs homologs (P < 0.005).
www.nature.com/scientificreports www.nature.com/scientificreports/ Gene structure and motif analysis of the ZmPME/PMEI families. Gene structures of the ZmPME/PMEI genes were constructed by aligning the extracted genomic sequences to predicted cDNA sequences of the maize PME/PMEI genes. As can be seen from Supplementary Fig. S5, most of the ZmPME genes in subfamily I have 2 exons, and the ZmPME members in the subfamily III have 5 introns (except ZmPME28). In addition, most of the ZmPME genes contain 1-10 introns, and most of the ZmPMEI genes do not have no intron.
As expected, most closely related members have a common motif composition and gene structure pattern, which indicates functional similarity between the ZmPME/PMEI proteins in paralogous pairs or in the same subfamily (Fig. 3). For the ZmPME genes, proteins in the subfamilies IV and V contain motifs 1-5 (except ZmPME32, -34 and -43) and the intron phase 2, 0, 0 and 2 separating the PME domain (except ZmPME32, Fig. 3a,b). Proteins in the subfamily VII of the ZmPMEIs have motifs 8 and 9 (except ZmPMEI37 and -41), while that in the subfamilies II, III and IV have motifs 7, 8 and 9 (except ZmPMEI35 and -40, Fig. 3c,d), and most of ZmPME genes in the subfamily IV have 5 exons. The ZmPMEI genes are clustered into 8 subfamilies. The neighbor-joining (NJ) phylogenetic tree of the ZmPME/PMEI genes was constructed based on the amino acid sequence alignment of the proteins. Bootstrap values > 50% are indicated at each node. The neighbor-joining (NJ) phylogenetic trees of the ZmPME/PMEIs were built by MEGA7 (https://www.megasoftware.net/). www.nature.com/scientificreports www.nature.com/scientificreports/ The PME and PMEI domains of the ZmPME/PMEIs, the PME domains in E. chrysanthemi (GenBank: Y00549), carrot (SwissProt Accession No. P83218) and A. aculeatus (Swiss-Prot code Q12535); and the PMEI domains in kiwi (SwissProt Accession number No. P83326) and Arabidopsis (AthPMEI-1, Accession Number NP_175236; AthPMEI-2, Accession Number NP_188348) were analyzed by T-Coffee (http://tcoffee.org/503/ index.html) and displayed by ESPript 3.0 (http://espript.ibcp.fr/ESPript/cgi-bin/ESPript.cgi, Supplementary Figs. S9 and S10). The ZmPMEs contain five characteristic sequence fragments (44_GxYxE, 113_QAVAL, 135_QDTL, 157_DFIFG, 223_LGRPW; carrot numbering), and several highly conserved aromatic residues ( Supplementary Fig. S9). The ZmPMEIs contain four conservative Cys residues, which were connected by two disulfide bridges (first to second and third to fourth) and do not have the fifth conservative Cys residue which has a free thiol group comparing to kiwi and Arabidopsis.

Gene ontology (GO) annotation and subcellular localization of the ZmPME/PMEI proteins. The 92
ZmPME/PMEI genes (except ZmPME28) were assigned a total of 37 GO terms ( Fig. 5 and Table 1). Among them, 175, 78 and 167 proteins were assigned terms under molecular function, cellular component and biological process, respectively. Under biological process, 41 ZmPME genes were predicted to be involved in cell wall modification, 42 ZmPME genes were related to pectin catabolic process and 73 genes (all of ZmPMEI genes and ZmPME21, -22, -23 and -24) were involved in negative regulation of catalytic activity. Under cellular component, 41 ZmPME genes were assigned to cell part. Under molecular function, most of the ZmPME genes and a few ZmPMEI genes had pectinesterase activity and aspartyl esterase activity; and most ZmPMEI genes and a few ZmPME genes had pectinesterase inhibitor activity or enzyme inhibitor activity. In addition, we analyzed the GO annotations for each subfamily. The same annotations exist in different genes of different subfamilies (e.g., the PMEI genes in the subfamilies II and III), and there are also different annotations for genes in the same subfamily (e.g., the PME genes in the subfamily II, Supplementary Table 4). These results suggested that different genes in the same subfamily may have different roles in the evolution process.
Expression assay of the ZmPME/PMEI genes. The NimbleGen maize microarray data 41 (ZM37) including 60 tissues representing 11 major organ systems and various developmental stages of the B73 maize inbred line was employed to analyze the expression pattern of the ZmPME/PMEI genes. All of the ZmPME/PMEI genes (except www.nature.com/scientificreports www.nature.com/scientificreports/ one ZmPME gene and 13 ZmPMEI genes) expression data was used to draw Heatmap. Of them, 35 ZmPME/PMEI genes had a much higher expression level in anthers than in other tissues (Fig. 2), these ZmPME/PMEI genes may be related to the development of anther or pollen. In general, expression pattern was similar for genes within the same paralogous gene pairs (e.g., ZmPMEI1/-17, Supplementary Table 1 and Fig. 2), indicating they might be formed by segmental duplication and retained their function. However, the expression profiles of the four paralogous gene pairs (ZmPME12/-13, ZmPME14/-15, ZmPME16/-17 and ZmPMEI29/-38) were fundamentally different in different tissues, suggesting that these genes may have differentiated with different roles.
To confirm the organ-specific expression of ZmPME/PMEI genes shown by the microarray data, 14 ZmPME/PMEI genes specifically expressed in anthers, ZmPME24 (not including in the maize microarray data) and ZmPME30 (expressed low in all tissues) were selected for conducting semiq-RT-PCR. Semiq-RT-PCR was  www.nature.com/scientificreports www.nature.com/scientificreports/ performed with total RNA isolated from the roots, leaves, ears, immature tassels, pollens, anthers, whole seed (after pollinated), endosperm and embryo of the B73 inbred line. Fourteen of them, ZmPME3, -5, -7, -11, -23, -31, -42, -43 and ZmPMEI2, -16, -25, -31, -32, -44 matched well with the microarray data, all of these ZmPME/PMEI genes expressed significantly higher in the anthers or pollens than in other organs. It is interesting to note that ZmPME23 was specifically expressed in the anthers. The gene not included in the microarray data, ZmPME24, was also specifically expressed in the anthers and pollens. However, expression of only one gene, ZmPME30, did not match with the microarray data, it had higher expression level in the anthers and pollens but showed little or no expression in roots, leaves, ears and seed according to the semiq-RT-PCR assay (Fig. 6a).
To further analyze the possible role of the ZmPME/PMEI genes on anther or pollen development, expression of the 13 genes (which had identified specifically expression in the anthers or pollens of B73 inbred line) was investigated in the anthers of three fertile and three sterile individuals of a maize backcrossing population derived from a cytoplasmic male sterility (CMS) line. The results showed that all of the selected ZmPME/PMEI genes were differentially expressed in the fertile and sterile anthers, although some ZmPMEI genes showed lower expression in the sterile anthers (Fig. 6b).

Discussion
Genome-wide analysis has identified the PMEs and PMEIs in nearly all vascular plants and in multiple gene members 18 . Up to now, the function of a number of the PME genes have been studied in Arabidopsis 42 , rice 9 , pea 43  www.nature.com/scientificreports www.nature.com/scientificreports/ and found to be involved in pollen tube elongation 37,38,51,52 . Thus genome-wide identification, evolution, and expression analysis of the PME/PMEI families in maize will facilitate to understanding of the function of the gene families.
In this study, 43 and 49 ZmPMEs and ZmPMEIs were identified in the maize genome, which were divided into 5 and 8 subfamilies (Supplementary Table 1; Fig. 1), respectively. The number of ZmPMEs is less than that identified in Arabidopsis 7 , and Gossypium raimondii 11 , while that of ZmPMEIs is much more than that identified in Sorghum bicolor 12 , and Brachypodium distachyon 53 . We identified 24 paralogous pairs in maize, but all Ka/Ks values of paralogous are less than 1 (Supplementary Table 2), implying that ZmPME/PMEI genes evolved mainly under the influence of stabilizing selection. The result of Ka/Ks analysis of PME, PRO and PMEI domains reveals that the homologous gene pairs of Arabidopsis, rice and sorghum experienced purifying selection 12 and the PME homologous gene pairs of G. arboreum, G. raimondii and G. hirsutum also experienced stabilizing selection 45 . Thus ZmPME/PMEI genes might play important role in growth and development of plants.
Intragroup ZmPME/PMEIs have conserved gene structure and motif composition, indicating that ZmPME/ PMEIs in the same group could have the same function and they might come from a common ancestor. For instance, the ZmPME subfamilies IV and V contain motifs 1-5 (except ZmPME32, -34 and -43) and most intron phases 2, 0, 0 and 2 separating the PME domain (Fig. 3a,b). Of the ZmPMEIs, most members in the subfamilies III, IV and VIII expressed higher in anthers than in other tissues (Fig. 2), and most members in the subfamily VII have motifs 8 and 9 (Fig. 3c,d). The PME and PMEI domains alignment implying that the ZmPME domains have five characteristic sequence fragments (44_GxYxE, 113_QAVAL, 135_QDTL, 157_DFIFG, 223_LGRPW; carrot numbering, Supplementary Fig. S9), which have all been shown to be functionally important in carrot 54 ; the ZmPMEI domains have four conservative Cys residues ( Supplementary Fig. S10), which are connected by two disulfide bridges, but do not have the fifth conservative Cys residue in comparison with that in kiwi 29 and Figure 6. Semiq-RT-PCR analysis of some ZmPME/PMEI genes in the nine different tissues of B73 inbred line and in the anthers of three fertile and three sterile individuals of a maize backcross population. (a) The total RNA of nine tissues including seedling roots, seedling leaves, 2-cm immature ears, non-emerged immature tassels, anthers, pollen, whole seeds after pollination for 18 days, endosperm and embryo after pollination for 20 days of B73 inbred line was isolated and used to perform the semiq-RT-PCR of the ZmPME/PMEI genes. (b) The total RNA of anthers of three sterile and three fertile individuals was isolated and used to perform the semiq-RT-PCR of the ZmPME/PMEI genes. Beta actin was used for internal controls to normalize the RNA contents in each sample. Primers used are shown in Supplementary Table 5. In the figure, the PCR products were separated with the same experiment condition and that gels were processed in parallel. The original gels were presented in Supplementary Fig. S11 www.nature.com/scientificreports www.nature.com/scientificreports/ Arabidopsis 55 . The structure of the carrot PME is almost completely superimposable to the structure of tomato 30 . In this study, the 3D structure of ZmPME3 and ZmPMEI2 are highly similar to that from Carrot 40 and kiwi 29 , respectively. This indicated that the PME and PMEI domains were highly conserved in different plant species.
In recent years, there are many reports on the function of PME/PMEI genes. Overexpression of PMEI5 in Arabidopsis thaliana caused aberrant growth morphology of the stems 56 ; VvPMEI1 expression negatively correlates with the PME activity during the early stage of grape berry development of Grapevine 57 . The expression pattern of the ZmPME/PMEI genes from the NimbleGen maize microarray data showed 35 ZmPME/PMEI genes had a much higher expression level in the anthers than in other tissues (Fig. 2), and semiq-RT-PCR analysis of different tissues from B73 inbred line verified that they had much higher expression in the anthers or pollens (Fig. 6a). In addition, semiq-RT-PCR analysis of some ZmPME/PMEI genes showed that 13 ZmPME/PMEI genes were differentially expressed in the anthers of fertile and sterile individuals derived from a maize S-type CMS line (Fig. 6b). Similar results also have been reported in other plant species. For example, antisense expression of a pollen-specific PMEI from broccoli (Brassica oleracea) in Arabidopsis triggered silencing of the orthologous Arabidopsis gene At1g10770 and resulted in male sterility 48 ; the expression of the pectin methylesterase gene (At3g06830) was significantly lower in male-sterile line than in male-fertile line at the <1 mm anther length stage of Brassica napus 58 ; in cotton, transcriptome analysis showed that many pectin methylesterase genes highly expressed in flowering buds of fertile plants compared to those of the CMS-D8 line 59 . This implied that a number of ZmPME/PMEI genes might play an important role in anther and pollen development, however, their detailed roles in male function of maize need to be further studied in future.

Analysis of gene structures and conserved motif of the ZmPME/PMEI genes. Several ZmPME/PMEI
genes had more than one gene model annotated in MaizeGDB (https://www.maizegdb.org/). To confirm the putative alternative splicing transcripts, transcript-specific primers (Supplementary Table 5) were designed to amplify corresponding DNA isolated from B73 seedlings and cDNA derived from B73 pollen RNA. Conserved PME/ PMEI domains and gene structures producing validated transcripts were drawn and displayed using the Gene Structure Display Server 60 (GSDS2.0; http://gsds.cbi.pku.edu.cn/index.php).
We used the MEME system (http://meme.sdsc.edu/meme/itro.html) to identify conserved motifs with parameters set as: number of repetitions, arbitrary; maximum number of patterns, 20; optimal width of the motif, between 6 and 50 residues 64 .

Calculation of synonymous (Ks) and non-synonymous (Ka) substitutions.
To identify homologous pairs of genes, the transcript sequences of the ZmPME/PMEIs were investigated by BLASTN searches 65 . Paralogous pairs within the genome of maize were defined as follows: the aligned sequences were longer than 300 bp and shared identities ≥40% 66 . If the amino acid shorter than 300 bp, the aligned region had an identity ≥60% and the alignment length covered ≥50% of the gene were defined paralogous pairs. Gene ontology (GO) annotation. The translated ZmPME/PMEIs protein sequences were annotated using the Blast2GO5.2.4 program to assign GO terms 67 (http://amigo.geneontology.org/amigo/term/). GO analysis e-value is 1.0E-6. GO terms are provided under three main categories, biological process, cellular component, and molecular function.