Prediction of RNA editing sites and genome-wide characterization of PERK gene family in maize (Zea mays L.) in response to drought stress

2.1 PERK gene family identification and characterization in maize genome

Whole maize genome sequence, as well as the general feature format file (GFF3), was downloaded from the Maize Genetic and Genomic Database (Maize GDB)(https://gamma.maizegdb.org). For the purposes of finding the probable candidates of PERK family in maize, the online Pfam database (https://www.sanger.ac.uk/Software/Pfam/) was used to download the PERK domain HMM profile and then subjected as a query into Blastp (Finn et al., 2014). For all of the retrieved protein sequences, the SMART tools (https://smart.embl-heidelberg.de/) were used to verify the presence of the PERK domain (Letunic et al., 2015). Maize PERK gene family sequence were downloaded from maize genome database. TAIR 10 (http: /https://www.Arabidopsis.org) was used to retrieved the arabidopsis sequences while all other sequences of studied organisms were retrieved from online plant database Phytozome version 11 (https://phytozome.jgi.doe.gov/pz/portal.html). ExPASyProtParam, (https://us.expasy.org/tools/protparam.html) an online web tools, were used to retrieve the physiochemical properties.

2.2 Sequence logos and phylogenetic/evolutionary analysis

The MEGA 7 software was used to find out the conserved sequences for amino acids. Sequence are aligned using ClustalW and the structure was constructed using TBtool (https.//github.com/CJ-Chen/TBtools). Furthermore, using this software the Neighbor-Joining method was used to get the phylogenetic tree to deduce the evolutionary history (Chothia et al., 2003). The distances of the number of amino acid sites in units were measured using the poisson correction parameters (Yang et al., 2008). The Bootstrap algorithm employed with 1000 repetitions to estimate the stability of the nodes in the phylogenetic tree. Total 98 amino acid sequences were used for this analysis.

2.3 Predicted protein motifs, structure of exon/intron and conserved domain analysis

To find preserved motif of PERK protein online web server Multiple Em for Motif Elicitation (MEME) is used (https://meme-suite.org/tools/meme). The TBtool was used to construct the motif structure using the MEME.xml file which is obtained through MEME suite. The default parameters were as follows: motif recurrence was set to 1 per sequence; frequency of motifs was set to 10; motif width was set to 5–50 residues; and the minimum number of motif sites was set to 5. Arrangement of Exon and Intron of PERK genes was investigated by using gff3 file downloaded from maize GDB. Structure is constructed using TBtool software. Afterward, the NCBI CDD tool (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) was used to perform the conserved domain analysis.

2.4 Chromosomal localization, gene duplication events and synteny analysis

The start and end position for each identified ZmPERK gene were procure from the Maize Genetic and Genomic Database (Maize GDB) and validated against the GFF3 file. Total chromosomal length is retrieved through using FASTA stat in TBtool. Finally, using MapChart v2.32 (Voorrips, 2002), the ZmPERK genes were spatially mapped onto the maize chromosomes. The phylogenetic tree was used to identify the putative paralogous PERK gene pairs. The resulting pairs were subjected to TBtool software to determine synonymous and non-synonymous substitution rates (Ka). To determine the nature of codon selection the Ka/Ks ratio was also calculated that allegedly occurred during evolution. Further, using the formula T = Ks/2 and assuming a clock rate of 6.05 X10 9 substitutions/synonymous site/year for maize, the approximate period of duplication event was calculated (Kong et al., 2013). The genome sequence files, and gene annotation files (GFF3) of sorghum, rice and maize were used for the collinearity analysis. Required files generated using one step MC scan. TBtool software was used to visualize the results, and the parameter filtering genes in the collinearity block was set to 40.

2.5 RNA editing sites prediction, subcellular localization, and promoter analysis

RNA editing is a method in which certain cytidines in mitochondrial and chloroplast transcripts of plants are converted to uridines. The online web server like PREP-Cp (for chloroplast genes) and PREP-Mt (for mitochondrial genes) software (https://prep.unl.edu/) with the cutoff value to 0.8 were used in predicting the RNA editing sites (Mower, 2009). Location of genes at cellular level were also predicted using online web server softberry (https://www.softberry.com). In order to perform promoter analysis, the 5′ upstream region of each gene of the ZmPERK was downloaded from NCBI (https://www.ncbi.nlm.nih.gov/) and the resulting file was submitted to the online database plantCARE (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) for a cis-element scan.

2.6 In-silco expression analysis

The NCBI Geodataset (https://www.ncbi.nlm.nih.gov/gds) was used to obtain transcriptome data (GSE136087, GSE40070). RNA-sequencing analysis was performed on seed embryos from 15 days after pollination (DAP) to 40 days after pollination (DAP) (Sekhon et al., 2013). We also retrieved drought stress expression data. Inbred line B73 maize (Zea mays) plants grown in the green house under well-watered and drought stress circumstances until they reached the reproductive stage (at the onset of silk emergence). Plants were hand pollinated for drought stress two to three days after irrigation was stopped, and measurements and samples were collected 24 h later for transcriptome analysis. This gene expression data was used to construct heat map.

3.1 Identification and characterization of PERK gene family in maize genome

A systematic approach was followed to identify the PERK protein encoding genes. Taking the advantage of publically available genome of maize, after removing redundant genes we excavated 23 PERK genes in maize genome and named them as ZmPERK1 to ZmPERK23. The PERK gene family comprises 15 Arabidopsis thaliana genes, 8 Oryza sativa genes, 14 Populus trichocarpa genes, 14 Sorghum bicolor genes, 9 Theobroma cacao genes, 23 Z. mays genes, 10 Ananas comosus genes, and 5 Physcomiterlla patens genes. The biophysical properties of the ZmPERK family genes were then determined, including genes ID, start and end positions of genes on chromosomes, polarity of strand, length of CDS sequence (bp), length of protein sequence (aa), protein molecular weights (MW), isoelectric points (pl), and predicted subcellular localization of ZmPERK genes. Table 1 presents all additional estimated biophysical properties. ZmPERK proteins had a peptide length ranged from 432 to 958 amino acids, with an average of 695 A.A (Table 1). The PI (Isoelectric point) of maize PERKs varied between 6.2 and 9.134, while the molecular weight ranged between 36.45 and 100.93 kDa, with an average of 68.69 kDa. The length of nucleotide, amino acid sequences varied greatly, indicating that the ZmPERK genes are highly complex, implying a high level of complexity.

3.2 Sequence logos and phylogenetic/evolutionary analysis

Sequence logos analysis provide more comprehensive information for sequence similarities, significant alignment aspects, and sequence conservation patterns. To check PERK family evolution, we generate sequence logos and results showed that this family remained conserved throughout evolution. For comparison study, the protein sequences of maize, rice and Arabidopsis were used. The results reveal that consensus sequence residues were highly conserved, and there was no compositional bias seen across any specie. These results help in discover and analyze and evaluate PERK gene family protein sequence across the species.

Phylogenetic tree serves an important way to understand the evolutionary relationships pathways. In our study we created phylogenetic tree of PERK genes to depict the evolutionary relationships. The phylogenetic or evolutionary analysis revealed the oldest plant lineage, of the PERK gene family as its members were found in Ananas comosus (angiosperm), Physcomitrella patens (bryophytes), dicots (Arabidopsis thaliana, Theobroma cacao, Populus trichocarpa), and monocots (Oryza sativa, Sorghum bicolor and Z. mays). These findings suggested that these genes evolved in ancient land plants, and that probable orthologous genes can be found across the plant kingdom. The PERK genes were characterized by 29 members in the PERK-A clade, 26 members in PERK-B, 19 members in PERK-C, and 22 members in PERK-D in the phylogenetic study. PERK genes were randomly distributed in all four clades from dicot, monocot, and bryophytes plant species, indicating that these genes evolved after the split of bryophytes. This finding showed that PERK genes possibly expanded and diversified after the radiation of these different species. These evolutionary linkages can facilitate the identification of orthologous genes and help to accelerate their functional characterization.

3.3 Predicted protein motifs, structure of exon/intron and conserved domain analysis

The 23 ZmPERK protein sequences were classified into four subfamilies using a rectangular phylogenetic tree (subfamily I, II, III, IV). 10 members were found in the Subfamily I, followed by subfamily II (5), subfamily III (4), and subfamily IV (4) (Fig. 3A). In addition, we examined the conserved motifs using MEME software to further investigate the diversity of ZmPERK protein family (Fig. 3). In this study, total 10 motifs were found (Table S3). All the gene exhibits same motif pattern. The type, order, and number of motifs were consistent within a subfamily, but varied across subfamilies. The patterns of ZmPERK protein motif distribution revealed that conserved distribution patterns existed for similar motifs. Domains 1, 2 and 3 represent the distinctive protein kinase-binding domain that is found in all 23 ZmPERK proteins (Fig. 3C). Similarly, Fig. 3D depicts the relative lengths of introns and exon sequence conservation within each ZmPERK gene in maize. A gene's biological function is linked to the distribution of exons and introns. All these genes contain exons ranged between 2 and 10. The findings demonstrated obvious conservation, laying the groundwork for functional conservatism and guiding future functional research.

Close

Fig. 3

A: Phylogenetic tree-based categorization of ZmPERK genes. An un-rooted phylogenetic tree an un-rooted phylogenetic tree based on full-length peptide sequences (ZmPERK) was generated. Classification is shown based on a phylogenetic tree using differences into groups. 3B: Motif pattern of ZmPERK genes 3C: Conserved domains of maize PERK protein 3D: Exon–intron structure analyses of ZmPERK genes. The purple line represents introns, while the purple boxes represent exons.

Export to PPT

3.4 Chromosomal localization, gene duplication events and synteny analysis

Each ZmPERK gene's genomic DNA sequence was analyzed in the maize genome database using BLASTn to establish its location, and MapChart was used to visualize the position of identified ZmPERK members on their respective chromosomes. The chromosome map revealed that 23 PERK genes were dispersed out over 8 of the 10 chromosomes. (Fig. 4A). The most ZmPERK genes were found on chr06, which had six members, followed by chr01, 03, and 08, which had 4, 5, and 4 members, respectively. While chromosomes 2, 4, 7, and 10 each had only one gene. Zmchr06 had the most PERK genes (26.08%), followed by Zmchr3 (21.73%), Zmchr1, and Zmchr8 (17.39%), while Zmchr02, Zmchr04, Zmchr06, Zmchr07, and Zmchr10 had the lowest percentage (4.34%) (Fig. 4B).

Close

Fig. 4

A: Distribution of 23 ZmPERK genes on their respective chromosomes.4B. Pie chart representing percentage of genes present on chromosome.4C: Pictorial representation of paralog gene pairs on chromosome indicating the type of duplication either tandem or segmental.

Export to PPT

Gene duplications, either whole-genome or segmental, as well as tandem duplications, are critical for gene family evolution. Although it has been proven that segmental and tandem duplications play a key role in gene family evolution in all plants specie (Cannon et al., 2004).To investigate ZmPERK gene duplications and evolutionary processes in maize, we identified 09 pairs of probable paralogous genes using the maize PERK phylogenetic tree. It is well documented facts that tandem duplication observed when paralogous genes are present on the same chromosome, whereas segmental duplication arise when paralogous genes are located on distinct chromosomes (Panchy et al., 2016). All the paralogous gene pair appeared to have evolved by segmental duplication except one (ZmPERK5-ZmPERK9) which evolved through tandem duplication indicating that the evolution of PERK genes appears to have been dominated by segmental duplications in maize. (Fig. 4C) Segmental duplication is the primary force that drives the evolution of a gene family. The estimated time of divergence for paralogous gene pairs was determined using synonymous (Ks) and non-synonymous (Ka) substitution rates. The Ka/Ks ratios for all paralog ZmPERK varied between 0.10 and 0.65 (Table 2). It indicates that purifying selection may have been performed on codons in the development and proliferation of parallel PERK genes in maize.

Multiple collinearity scan tool was used to find orthologous genes among genomes of maize, Sorghum, and rice to further understand the Synteny links of ZmPERK genes with these plant species. (Fig. 5). 18 pairs of collinearity genes of PERK gene family between maize and rice whereas sixteen pairs in maize and sorghum were observed in the synteny analysis (Table S2). Gene IDs of all collinear genes is given supplemental file. According to these findings, the collinearity between maize and sorghum is significant as compared to the collinearity values between maize and rice furthermore, these PERK genes in maize derived from a common ancestor.

Close

Fig. 5

Collinearity analysis of maize, rice, and sorghum. (A) Collinearity analysis of all chromosomes reveals duplicated PERK genes in maize and sorghum. The lines connect the pairs of duplicated genes. (B) The collinearity study of maize and rice chromosomes. The PERK genes are represented by the red flags on distinct chromosomes.

Export to PPT

3.5 Prediction of RNA editing sites, subcellular localization, and promoter analysis

The Prep-CP and Prep-Mt prediction tools were used to find the RNA editing sites of ZmPERK chloroplast and mitochondrial genes, respectively. In chloroplast genes, 196 RNA editing sites were predicted (Table 3A) and 268 in mitochondrial genes (Table 3B). All predicted editing sites in the chloroplast and mitochondrial genomes were distributed among 23 genes, with an average of 8.5 and 11.26 editing sites per gene, respectively. The chloroplast gene ZmPERK1 contains maximum RNA editing sites (12) while minimum editing sites (5) were predicted in ZmPERK14. Similarly, mitochondrial gene ZmPERK12 contain maximum editing sites (21) while minimum sites (7) were predicted in ZmPERK14. The position of RNA editing sites was further explored, and it was observed that all the predicted sites were based on first and second codon base changes. At the third codon base, we couldn't find any site for RNA editing. The transition of cytosine → uracil (C-U) seems to be present in all editing sites, resulting in amino acid substitutions. Eleven type of amino acid change found in chloroplast and mitochondrial genes (Fig. 6A). Amino acid conservation caused by RNA editing including A (Alanine) → V (Valine), T (Threonine) → I (Isoleucine), H (Histidine) → Y (Tyrosine), P (Proline) → S (Serine, P (Proline) → L (Leucine), R (Arginine) → C (Cysteine), S (Serine) → F (Phenylalanine), R (Arginine) → W (Tryptophan), P (Proline) → F (Phenylalanine), S (Serine) → L (Leucine), T (Threonine) → M (Methionine), L (Leucine)F (Phenylalanine).

Close

Fig. 6

(A) RNA editing of the PERK genes results in amino acid conservation. (B) Identified cis-acting elements in ZmPERK gene family promoters.

Export to PPT

Location of gene at cellular level was also determined. Results indicated that 21 of the 23 ZmPERK proteins were localized to the plasma membrane, while two (ZmPERK18 and ZmPERK23) were localized to the cytoplasm, Table1 contains the details of these parameters. The promoter region, which is located upstream of the start codon area, controls gene transcription. Understanding gene regulation and function requires a thorough examination of cis-elements (Higo et al., 1999). We discovered and classified cis-acting factors in the upstream region of the ZmPERK genes. (See Table S4) The cis-elements were categorized based on their roles in growth and development, as well as light and stress actions. The upstream region of ZmPERK genes contained cis-acting factors like MeJA responsive, MYB-binding sites associated with light responsive elements, ABA responsive elements, defense, stress, low temperature, gibberellin acid (GA), and salicylic acid (SA)responsive elements, as per the promoter analysis results (Fig. 6B). The promoters of the ZmPERK gene have the most MeJA responsiveness elements. where they were found in 296 promoters. There were 166 light responsive elements, 88 ABA responsive elements, 27 GA responsive elements, 25 MYB light responsive elements, and 4 auxin responsive elements. The cis-element analysis showed that during abiotic stress and plant development phase the ZmPERK genes could respond.

3.6 In-silico expressions analysis

Expression patterns give information regarding the biological activities of genes because gene expression is required for optimal regulation of plant growth and development. We looked examined the expression patterns of the ZmPERK under drought stress and at different stages of seed embryo development from 15 days after pollination (DAP) to 40 days after pollination (DAP) in two distinct varieties (High oil content and low oil content) (Fig. 7A).Under drought stress 9 genes were upregulated in both tissues leave and cob (ZmPERK2, ZmPERK 3, ZmPERK 8, ZmPERK 10, ZmPERK 14, ZmPERK 16, ZmPERK 19, ZmPERK 20) 4 genes (ZmPERK 1, ZmPERK 7, ZmPERK 18, ZmPERK 21) only upregulated in leave tissue and 1 gene ZmPERK 18 upregulated in cob (Fig. 7B). While Seven genes downregulated (ZmPERK 4, ZmPERK 5, ZmPERK 11, ZmPERK 12, ZmPERK 15, ZmPERK 17, ZmPERK 22) under drought stress.For oil content accumulation in embryo 9 genes shows upregulated expression pattern (ZmPERK 1, ZmPERK 3, ZmPERK 6, ZmPERK 8, ZmPERK 13, ZmPERK 14, ZmPERK 16, ZmPERK 19 ZmPERK,20) while 12 genes show downregualted trend (ZmPERK 2, ZmPERK 5, ZmPERK 6, ZmPERK 7, ZmPERK 9, ZmPERK 10, ZmPERK 11, ZmPERK 12, ZmPERK 17, ZmPERK 18, ZmPERK 22, ZmPERK 23). Interestingly, some genes show similar expression pattern under both conditions. For example, ZmPERK 3 ZmPERK 8 ZmPERK 14 ZmPERK 16 ZmPERK 19 ZmPERK 20 regardless of the tissues or stresses applied, they were always upregulated. We may conclude from these findings that ZmPERK gene expression is involved in drought stress and oil content accumulation in the embryo.

Close

Fig. 7

A: ZmPERK gene expression patterns at different developmental stages of seed embryos from 15 days after pollination (DAP) to 40 days after pollination (DAP). H-represent expression in high oil content varaiety while L represent low oil content varaiety B: Expression pattern of ZmPERKgene under drought stress DS-L(Drought stress leave sample) DS-C (Drought stress cob sample) CL(Controled leave sample) CC (Control Cob sample.

Export to PPT