Evolutionary Understanding of Metacaspase Genes in Cultivated and Wild Oryza Species and Its Role in Disease Resistance Mechanism in Rice

Metacaspases (MCs), a class of cysteine-dependent proteases found in plants, fungi, and protozoa, are predominately involved in programmed cell death processes. In this study, we identified metacaspase genes in cultivated and wild rice species. Characterization of metacaspase genes identified both in cultivated subspecies of Oryza sativa, japonica, and indica and in nine wild rice species was performed. Extensive computational analysis was conducted to understand gene structures, phylogenetic relationships, cis-regulatory elements, expression patterns, and haplotypic variations. Further, the haplotyping study of metacaspase genes was conducted using the whole-genome resequencing data publicly available for 4726 diverse genotype and in-house resequencing data generated for north-east Indian rice lines. Sequence variations observed among wild and cultivated rice species for metacaspase genes were used to understand the duplication and neofunctionalization events. The expression profiles of metacaspase genes were analyzed using RNA-seq transcriptome profiling in rice during different developmental stages and stress conditions. Real-time quantitative PCR analysis of candidate metacaspase genes in rice cultivars Pusa Basmati-1 in response to Magnaporthe oryzae infection indicated a significant role in the disease resistance mechanism. The information provided here will help to understand the evolution of metacaspases and their role under stress conditions in rice.


Introduction
Plants are able to engage in programmed cell death, which can be defined as a schematic mode of action that mediates cell death both under stress and at developmental phases of life [1,2]. The phenomenon of programmed cell death is thoroughly characterized in higher animals where caspases are found to have the key role [3,4]. In contrast, very little is known about the programmed cell death in plants and fungi, which lack caspases [5]. Instead of caspases, metacaspases and paracaspases, which are structurally related to caspases, are present in plants and fungi [6]. The metacaspases and paracaspases belong to the C14 family of proteases. Interestingly, metacaspases can be found in all eukaryotic lineages except for a strict absence in animals [7]. Numerous studies have been performed genome. The O. nivara as a separate species is not yet established, but considering the whole genome resequencing of reference genome, here we considered it as a species. The retrieved protein data were used to create a local database in BioEdit [28]. Subsequently, BLASTp search was performed using nine metacaspase genes (AthMC1-AthMC9) reported in Arabidopsis [29]. To identify significant matches, BLASTp was performed with default parameters with an e-value lesser than 0.00001, and the bit score was set to be greater than 100. All the insignificant hits were excluded. The selected sequences were confirmed by screening for the peptidase domain (Peptidase C_14, pfam00656) by subjecting the sequences to Conserved Domains Database (CDD) [30].
The nomenclature of the selected genes was performed on the basis of initials of the rice species. The first letter "O" signifies Oryza, second and third letters stand for the species name, for example

Gene Structure, Motif Analysis, and Phylogenetic Tree Construction
The Gff3 (General feature format) files were used to retrieve the information regarding the exons and introns in the metacaspase genes. The gene structure was visualized by Gene Structure Display Server 2 [31]. Similarly, the MEME Suite was used to identify the conserved motifs in the metacaspases protein sequences [32]. Multiple sequence alignment was performed by ClustalW provided in the MEGA (version 10.0) software tool. The phylogenetic tree was constructed for protein sequences of 92 Oryza metacaspase genes by the Neighbor-Joining Tree Method with JTT + G + I model with bootstrap value 1000 replicates in Mega X software [33].

Chromosomal Localization, Gene Synteny, and Cis-Regulatory Element Analysis
On the basis of genomic annotations retrieved from EnsemblPlants, the chromosomal locations of metacaspases were visualized by Mapchart2.2 [34] and Circos tool [35]. Along with the chromosomal localization, the relationship between the genes was visualized using the Circos tool. The 2 kb upstream genomic region was retrieved for all the 92 metacaspase genes and subsequently analyzed by PlantCARE and PLACE databases to identify the transcription factor binding elements [36].

Localization and Tertiary (3D) Protein Structure Analyses
Various in silico tools, including CELLO [37], TMHMM 2.0 [38], and ProtParam [39] were used to study the localization of metacaspase proteins. A gene annotation study was conducted to annotate the roles of metacaspase proteins. The protein tertiary structure for the metacaspases was developed using I-TASSER (Iterative Threading ASSEmbly Refinement) server [40].

Haplotypic Evaluation
Haplotypic evaluation of the metacaspase genes was performed based on the sequence variations, including single-nucleotide polymorphisms (SNPs) and InDels. Whole-genome resequencing information available for 4500 diverse rice genotypes provided at RiceVarMap2.0 was used for the haplotyping of the metacaspase genes [41]. A haplotyping network was developed to understand the evolution of different haplotypes in respective rice groups such as O. sativa subsp. indica and japonica, and tropical.

Transcriptomic Evaluation and Co-Expression Analysis of Metacaspase Genes
Transcriptomic data (RNA-Seq) were retrieved from the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) to study the expression pattern of metacaspase genes. The data related to developmental, abiotic, and biotic stress conditions such as Magnaporthe oryzae infection, cold, drought, and salt stress conditions in susceptible and resistant cultivars were downloaded in the fastq format. CLC Genomic Workbench 12 (http://resources.qiagenbioinformatics.com/) was used to obtain the expression values in terms of reads per kilobase of transcripts, per million mapped reads (RPKM). The raw reads were processed and then mapped to the respective genomes with default parameters. The RPKM values were then used to develop the co-expression network by CoExpNetViz. The co-expression network was visualized using Cytoscape 3.7.1. Highly co-expressed genes were highlighted with a larger label size. Furthermore, log2-transformed RPKM expression values were used to generate a heatmap to show the expression pattern in different conditions by performing hierarchical clustering analysis using the Euclidean distance metric and average linkage clustering linkage method implemented in Multiple Experiment Viewer (MeV) v 4.9.0 software.

Plant Materials and Fungal Infection
To study the response of metacaspase genes family in compatible and incompatible interactions, a resistant rice genotype Tetep and a susceptible genotype Pusa Basmati-1 were subjected to infection with Magnaporthe oryzae race Pani NB. The plants were grown in a greenhouse under controlled environment conditions (16 h light/8 h dark at 25 • C) for three to four weeks. The M. oryzae infection assay was performed as earlier described by Valent et al. [42], with some modifications. At the 3-4 leaf stage, each plant was inoculated with 20-25 mL of blast spore suspension (10 5 spores/mL). To ensure uniform application of inoculum, an atomizer sprayer was used. The plants were maintained in complete darkness and 100% relative humidity for the initial 24 h. Later, the photoperiod was adjusted to the 16 h/8 h dark, keeping the relative humidity at 70-80%. Sampling was performed at five time-points, 0, 12, 24, 36, and 48 h post-infection in three biological and technical replicates each.

RNA Extraction and Expression Analysis
Total RNA isolation was performed by using Spectrum TM Plant Total RNA Isolation Kit (Sigma Aldrich, St. Louis, MO, USA) following the manufacturer's protocol. The quality and the quantity of the total isolated RNA was checked on agarose gel electrophoresis and Nanodrop Lite UV-Visible spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Total RNA was used to synthesize cDNA using iScript TM cDNA Synthesis kit (Bio-Rad, Hercules, CA, USA) following the manufacturer protocol. The primers for the real-time quantitative qRT-PCR were designed by QuantPrime online tool (https://quantprime.mpimp-golm.mpg.de/). The qRT-PCR was performed in 10 µL PCR mixture (5 µL SYBRGreen, 0.3 µL forward and reverse primer each, 1 µL cDNA (5 ng/µL), 3.4 µL ddH 2 O). The house-keeping gene Ubiquitin was used as a reference to normalize the expression values of target genes. The relative expression calculations were performed using the 2 −∆∆Ct formula [43] for three biological and technical replicates each.

Metacaspases in Wild and Cultivated Rice Genomes
A homology-based search performed in ten Oryza species initially identified 93 metacaspase genes. Furthermore, the peptidase domain (PF00656) was identified in all genes (Supplementary File S1). One of the candidate genes (OBART11G02150.1) was omitted from the analysis despite the presence of the peptidase domain, due to the presence of a redundant sequence. The number of metacaspase genes was roughly similar in all ten wild and cultivated Oryza species under study ( Figure 1  Blue bars represent Type I metacaspase genes, and red bars represent Type II metacaspase genes.

Predicted Structural and Functional Attributes of Metacaspase Genes
The proteins in this family were found to be highly variable (Supplementary File S2, Table 1). Mostly the metacaspase protein length falls in the range of approximately 200-400 amino acids. The smallest metacaspase protein was found to be 120 amino acids long (OLoMC1); whereas, the largest protein found was 1086 amino acid long (ONiMC8). The average size and molecular weight of all 92 metacaspase genes under study were 404 amino acids and 40 kDa. The pI ranges between 5.24 and 11.59, which indicates both the acidic and basic nature of the protein. Mostly, all the metacaspase proteins are acidic or slightly acidic in nature, but a few among them possess basic characteristics. Similar to the molecular weight, the gravy score fluctuates between −0.713 and −0.008. The studies regarding the subcellular localization of metacaspase proteins in the cell suggested that localization of metacaspase proteins in the cell is also highly variable. Few metacaspase proteins showed the presence of signals for the plasma membrane, whereas a single protein contains the signal sequence for mitochondrial localization (OBrMC1). Gene annotation studies revealed the role of metacaspase proteins majorly in biological processes and molecular functions (Supplementary Figure S1). The three-dimensional structure studied by the iTASSER server revealed the similar architecture of metacaspase proteins with a high level of structural conservations (Supplementary File S3).

Gene Architecture, Conserved Motifs, and Phylogenetic Distribution of Metacaspase Genes
The evolutionary relationship among the metacaspases identified in the ten Oryza species was studied by phylogenetic tree construction. Metacaspase proteins on the basis of their domains were categorized into Type I and Type II groups, which were also grouped into two separate clades in the phylogenetic tree ( Figure 2). Out of the entire 92 metacaspases in rice species and 9 in Arabidopsis metacaspases under study, 54 and 47 proteins belonged to Type I and Type II, respectively (Figure 2).

Figure 2.
Neighbor-Joining phylogenetic tree of metacaspase genes from wild and cultivated Oryza species. The phylogenetic tree comprising 92 Oryza and 9 Arabidopsis metacaspase genes was generated using Mega X software. Red square tiles represent Type II metacaspase genes whereas green square tiles indicate Type I metacaspases. Gene nomenclature used OJa for japonica, OIn for indica, ORu for rufipogon, ONi for nivara, barthii is named as OBa, brachyantha as OBr, glumaepatula as OGl, glaberrima as OGb, OPu stands for punctata, OMe stands for meridionalis, and OLo is for longistaminata, and MC denotes the metacaspase family, which is followed by a numeral representing the specific gene.
A thorough study was conducted to analyze the gene structure and conserved domains present in the metacaspase proteins in the Oryza species. An exon-intron organization for all the metacaspase genes in wild and cultivated rice species was constructed to understand the variations in the coding and non-coding metacaspases sequences (Figure 3). The number of introns in the metacaspase genes is highly variable, spanning from as low as one to as high as thirty introns. Exceptionally, ORuMC9, OJaMC9, OInMC8, OGbMC9, OJaMC6 genes were found to be intron-less, whereas OPuMC8, ONiMC8, OMeMC8, OLoMC9 genes had more than 25 introns. The genes coding for longer metacaspases possessed a higher number of introns. Interestingly, the protein length of all genes with no intron was 341 amino acids with the exception of OJaMC6 (305 amino acids). A total of twenty conserved motifs were identified in metacaspase proteins using the MEME software tool (Figure 4). Subsequently, a CDD search was performed to annotate the identified motifs. Interestingly, the domain architecture in Type II metacaspases was very much similar to that in the Type I metacaspases group. Exon-intron architecture of metacaspases identified in ten Oryza species. Exons, introns, and untranslated region are indicated by green, grey, and red, respectively. The coding DNA sequence (CDS) and genomic sequences of metacaspase from all ten Oryza species were used to develop pictorial presentation using the GSDS 2.0 online tool. Gene nomenclature used OJa for japonica, OIn for indica, ORu for rufipogon, ONi for nivara, barthii is named as OBa, brachyantha as OBr, glumaepatula as OGl, glaberrima as OGb, OPu stands for punctata, OMe stands for meridionalis, and OLo is for longistaminata, and MC denotes the metacaspase family, which is followed by a numeral representing the specific gene. . Conserved motifs identified in metacaspase protein from ten Oryza species. The motifs were identified by MEME online tool (http://meme-suite.org/) using default parameters. Gene nomenclature used OJa for O. sativa subsp. japonica, OIn for O. sativa subsp. indica, ORu for rufipogon, ONi for nivara, barthii is named as OBa, brachyantha as OBr, glumaepatula as OGl, glaberrima as OGb, OPu stands for punctata, OMe stands for meridionalis, and OLo is for longistaminata, and MC denotes the metacaspase family, which is followed by a numeral representing the specific gene.

Chromosomal Distribution, Gene Duplication, and Syntenic Relationship of Metacaspase Genes in Wild and Cultivated Oryza Species
Metacaspase genes were very asymmetrically distributed on eight different chromosomes (Chr.) in Oryza species ( Figure 5). No metacaspase genes were found to be located on Chr. 2, Chr. 4, Chr. 6, and Chr. 7 in all ten Oryza species. Interestingly, 42% of total metacaspase genes, i.e., 39 out of 92, were found to be located on Chr. 3. Gene duplication events were analyzed to understand the duplication event that possibly happened during the evolution of metacaspase genes. The results suggested that metacaspase genes evolved through the segmental and proximal duplication events over tandem and dispersal duplications (Supplementary Table S1). Exceptionally, O. meridionalis showed dispersal duplication in all of its genes. The syntenic relation was also studied in order to gain a better understanding of the relationship between the metacaspase genes. As expected, a high level of synteny between metacaspase genes from different Oryza species was observed (Supplementary Figure S2). Figure 5. Distribution of metacaspase genes in ten Oryza species genomes. Green represents Type I metacaspase genes, and red represents Type II metacaspase genes. The outer circle represents the chromosomes of ten Oryza species carrying the metacaspase genes. The inner circle denotes the chromosomal location of all 92 metacaspase genes.

Cis-Acting Promoter Analysis
Promoter analysis was carried out in metacaspase genes identified in Oryza species to gain a better understanding of the gene regulation across the species (Supplementary File S4). Many stress, development, and light-responsive cis-regulatory elements were found in 2 kb upstream regions of the metacaspase genes along with the most commonly occurring cis-regulatory elements which include CAAT, CTCC, TATA, and TATATA box ( Figure 6). Light responsive G-box was present majorly in all the metacaspase genes of Oryza species; in contrast, Sp1 (GGGCGG), which is also related to light responsiveness, was found in some metacaspase genes. Identified cis-elements known to be participating in abiotic and biotic stresses include WRKY-related (S000447, S000390, and S000457), MYB-related (S000176), and MYC-related (S000407) [44,45]. Abscisic responsive elements known as ABRELATED1 (S000414), ABRERATICAL (S000408), and ABREOSRAB21 (S000012) and a GCC box (S000430) were also found to be located in the upstream regions of the metacaspase genes. Along with this an auxin-responsive element, NtBBF1 (S000273), which was found in a few metacaspase genes [46] (Supplementary File S4). The most frequently occurring stress-related cis-elements in all the Oryza species were Myc-related MYCCONSENSUSAT and WRKY71OS, whereas in O. glumaepatula they were MYCCONSENSUSAT and WBOXNTERF3. Auxin-related cis-element, NtBBF1 (S000273), was found in two or three members of each Oryza species except for O. brachyantha where the NtBBF1 promoter element was present in six metacaspase genes. Figure 6. Frequency of transcription factor binding motifs identified in the 2 kb upstream region of the metacaspase genes in ten Oryza species. The variation in size represents the abundance of the respective sequences.

Haplotyping for Genes in Rice
Haplotyping of Oryza metacaspase genes was performed using whole-genome sequencing information data obtained from RiceVarMap2.0 ( Supplementary Files S5 and S6). Consequently, haplotype networks were developed to depict the distribution of haplotypes across Oryza subgroups. As shown in Table 2, gene Os03t0389100-01 depicted the maximum number of non-synonymous mutations, whereas the minimum number of non-synonymous mutations was observed in Os10t0565100-00. Similarly, haplotypes based on the number of SNPs were found to be maximum in number for Os03t0389100-01 and minimum for Os05t0496400-00. The effects of amino acid variations detected in the metacaspase genes were further studied using the PROVEAN tool (Supplementary File S7). PROVEAN results revealed no deleterious mutations present in Os01t0799900, Os03t0389501, Os03t0389000, and Os11t0134700, whereas maximum deleterious mutations were present in Os03t389100.

Expression Analysis of Metacaspase Genes during Different Developmental and Stress Conditions
For a complete expression atlas of metacaspase genes, the RNA-seq data were retrieved from PRJNA134239, which comprises 43 independent RNA-seq libraries for O. japonica (cv. Nipponbare) tissues. The expression pattern of metacaspase genes was studied for abiotic, biotic, and developmental stages using data from different studies. No expression was observed for OJaMC3 across different RNA-seq libraries, therefore hierarchical clustering of eight genes across abiotic stress in O. sativa was performed. Expression analysis revealed categorization of these genes into two groups, wherein four genes showed moderate-to-low gene expression under drought, salt, and cold conditions (Figure 7a,  Supplementary File S8). The OJaMC5 gene has higher expression in leaves under drought and salt stress compared to that in control conditions. Similarly, OJaMC8 and OJaMC7 has higher expression in leaves under control condition compared to that in cold stress (Figure 7a). Gene expression clustering across different developmental stages and tissues of O. sativa revealed a similar clustering into two groups (Figure 7b). However, among eight genes, five clustered under a single group and the rest into the second group. Based on expression patterns, the group with five genes showed low-to-moderate gene expression, whereas the remaining group showed moderate-to-high gene expression for several developmental stages. The OJaMC5 gene showed higher expression in 60 days old leaves, stem, and mature stigma and ovary (Figure 7b). The comparative expression analysis of metacaspase genes in Oryza species was studied using the RNA-seq data retrieved from the bioproject PRJNA393480. The study was conducted on root tissue for Oryza species O. rufipogon (Dongxiang) and O. sativa japonica (Dongdao-4) when exposed to M. oryzae (Guy 11) under controlled environmental conditions in triplicates. The expression pattern of all nine metacaspase genes was studied in all control and test samples. The ORuMC7 gene was found to have high expression in both O. sativa subsp. japonica and rufipogon samples as compared to control, but the RNA-seq expression of ORuMC7 was much higher in wild species as compared to that in cultivated ones (Figure 7c).
The expression profiling of metacaspase genes against blast disease in the O. sativa indica resistant and susceptible varieties, Tetep and Pusa Basmati-1, respectively, were carried out in controlled lab conditions. The expression of all eight OInMC(1-8) genes was studied by qRT-PCR (Figure 8). The relative expression values of all of the OInMC genes were found to be increasing at different time intervals. A high relative log2 fold change value was obtained for OInMC1 and OInMC3 genes 12 and 36 h post-infection (hpi), respectively. The house-keeping gene Ubiquitin was used as a reference to normalize the expression values of target genes. The relative expression calculations were performed using the 2 −∆∆Ct formula [43] for three biological and technical replicates each.

Discussion
The metacaspase genes have been studied in several plant species to understand their role in stress and developmental phases of the growth cycle [29,47]. The number of metacaspases in plant species varies considerably, for instance, 6 in grapes, 8 in rice and tomato, 9 in Arabidopsis, and 30 in cotton species [18]. Such wide range of occurrence may be due to the high level of variation in duplicated genome and ploidy level in plants. The genome-wide study of metacaspases in the Oryza clade will help to better understand its functional role and evolution in rice. The in-depth understanding of metacaspases in wild relatives will also provide an opportunity to introduce desirable alleles from wild genetic resources. Here, identification, characterization, and distribution of metacaspase genes in ten Oryza species were performed, which uniformly identified 8 or 9 homologs in each species despite having differences in the origin of the genomes (AA, BB, or FF) (Figure 1, Table 1). Recently, genome-wide identification of metacaspase genes performed in four Gossypium species identified a varying number of homologs: 16 in G. raimondii, 17 in G. arboreum, 26 in G. hirsutum, and the highest, 30, in G. barbadense [48]. Such a wide variation among Gossypium species as compared to the Oryza species is mostly due to the higher ploidy level and proportion of the duplicated genome.
The metacaspase genes were classified into Type I and Type II metacaspases on the basis of the linker region between p10 and p20 domains ( The zf domains are thought to be involved in protein-protein interactions. A study performed in yeast has shown active involvement of zf domains in the interaction of AtMC1 with LSD-1 in the case of the Oryza genus [49]. Therefore, the present/absent variation for the zf domains amongst species will help to study the significance of this domain. Information about the exon-intron organization of metacaspase genes helps to predict the origin, expression activity, as well as possibilities of alternate splicing [50]. The number of introns was higher in Type I metacaspase genes as compared to that in Type II MC (Figure 3). The number of introns was comparatively very high in O. longistaminata, O. punctata, O. nivara, O. meridionalis wild species, whereas three or four introns and intron-less metacaspase genes are common in O. sativa (japonica and indica), O. rufipogon, O. barthii, and O. brachyantha. The exon-intron organization study was consistent with the earlier reported studies about the number of introns higher in Type I than in Type II metacaspases [13,51]. Gene duplication plays an important role in the loss of introns. In general, recently evolved genes tend to have fewer introns. However, this is very controversial to say whether intron loss is dominant in higher plant or intron gain [50]. Nevertheless studies in rice clearly indicated that the newly evolved genes tend to have fewer introns and have more chances to lose their expression during the evolution [50]. Similarly, the organization of conserved motif also helps to understand the functional domains as well as the evolutionary relationship of the metacaspase gene family members. Here, a total of 20 conserved motifs identified with the MEME tool (Figure 4), including the signature p10 and p20 domains, and the overall organization looked highly conserved in the two types of metacaspases.
In the present study, whole-genome resequencing data from 4726 diverse genotypes were accessed to get the sequence variation in metacaspase genes. The sequence analysis helped to locate all the SNPs and InDels and subsequently predicted their effect on protein functionality ( Table 2). The deleterious mutations identified in the present study will be helpful in understanding the role of metacaspases. Earlier, such deleterious and neutral variation predicted with the PROVEAN tool was successfully validated using a site-directed mutagenesis experiment performed with SlNIP2-1, PtNIP2-1, and OsNIP2-1, respectively, from tomato, poplar, and rice [52]. Such predicted deleterious variation in the metacaspase genes need to be experimentally validated. The publicly available whole-genome resequencing data provide an opportunity to understand the functional role of metacaspase genes. The sequencing data were also used to define the haplotypes, which will be helpful to develop haplotype-specific markers and also to understand the haplotypic evolution (Supplementary File S5). The haplotypes representing allelic variation existing in the rice germplasm allow us to study the effect of each allele. The variations located in conserved domains can also help to elucidate the molecular mechanism involved in metacaspase activity.
Studying gene expression is essential to understand the tissue-specific or stress-specific function of metacaspase gene family members. Numerous studies have shown stress-induced expression of several metacaspase genes in different plants [47,48,53]. In the present study, available RNAseq transcriptome profiling was used to identify tissue-and stage-specific, stress-specific as well as differential expression among different Oryza species of metacaspases genes (Figure 7). The RNA-seq transcriptome data of control and different stress conditions clearly showed that the some metacaspase genes highly expressed in control also had high expression throughout all the stress conditions. Such abundant expression of some metacaspase genes is surprising. In this regard, a recent study where the crystal structure of Arabidopsis metacaspase 4 (AtMC4) was evaluated seems helpful to explain why some metacaspases have ubiquitous expression [54]. The structural analysis of AtMC4 revealed activation of the metacaspases only in the presence of Ca 2+ [54]. Such a Ca 2+ -dependency explains several phenomena that were elusive before. The study suggests that the highly expressed metacaspase genes are non-functional until an adequate supply of Ca 2+ occurs. This helps the plant to respond quickly to stress conditions. For instance, Ca 2+ stored in different cellular compartments get fluxed in the cytoplasm with stress signals, and then Ca 2+ -dependent activation of metacaspase takes place. The constitutive high expression of some of the metacaspases will be surely helpful in Ca 2+ -dependent stress response mechanism. Such a phenomenon may be dominantly observed with metacaspase, but many of the metacaspases family members may have a very specific mode of action. To understand such a tissue-and condition-specific role of metacaspases, more extensive studies using integrated omics approaches are needed.
In addition to extensive RNAseq transcriptome analysis, qRT-PCR analysis of metacaspase genes under M. oryzae infection was performed in the present study ( Figure 8). The metacaspase genes were found to be differentially expressed after M. oryzae infection, but the response time for each of the genes was unique. For instance, OInMC1 showed upregulation at 12 hpi, OInMC8 at 24 hpi, whereas OInMC4 showed upregulation at 48 hpi. This suggests the possible role of all the metacaspase genes in M. oryzae resistance mechanism, but they may have a different role in different cellular compartments or at different steps of the stress-response cascade. Similar results have been reported in rubber tree (Hevea brasiliensis) where most of the metacaspase genes showed differential expression patterns under various stresses [55]. For instance, eight out of nine metacaspase genes were found to be differentially expressed under drought conditions in rubber trees [55].

Conclusions
The role of metacaspase genes under stress and developmental processes in plants has been demonstrated in several studies. However, precise molecular mechanisms involving metacaspase genes are not yet well understood. In the present study, 92 metacaspase genes were identified in nine Oryza species and two subspecies of cultivated species O. sativa. Extensive gene characterization performed in this study will help to explain the functional role of metacaspases. The phylogenetic analysis performed here helped to define the distribution of Type 1 and Type 2 metacaspases in the Oryza clade. The distribution of Type 1 and Type 2 across the Oryza species is uniform despite having different genomes (AA, BB, or FF). The sequence variations identified in the metacaspase genes will be helpful for the understanding of structural features of the protein, role of conserved domains, allelic effects, and overall functionality of metacaspase genes. The defined haplotypes for metacaspase genes will be helpful to develop markers for breeding applications as well as for genetic studies. In addition, the present study has provided information about the metacaspase gene expression under biotic and abiotic stress conditions, in different tissues and developmental conditions. Metacaspase gene expression was also evaluated with a real-time qPCR assay, which suggested its role in different stages of disease development. The stepwise upregulation and downregulation of different metacaspases indicate the possibility of a metacaspase-regulated cascade of systematic stress response. The detailed information of evolution, distribution, expression, and allelic variation provided here will be helpful to better understand the metacaspases in rice and for the exploration of knowledge in crop improvement programs.