In silico identification and validation of miRNA and their DIR specific targets in Oryza sativa Indica under abiotic stress

Several biotic (bacterial and viral pathogenesis) and abiotic stress factors like salt, drought, cold, and extreme temperatures significantly reduce crop productivity and grain quality throughout the world. MicroRNAs (miRNAs) are small (~22 nucleotides) non-coding endogenous RNA molecules which negatively regulate gene expression at the post-transcriptional level either by degrading the target protein-coding mRNA genes or suppressing translation in plants. Dirigent (DIR) gene protein plays a crucial role as they are involved to dictate the stereochemistry of a compound synthesized by other enzymes as well as in lignifications against biotic and abiotic stress. In plants, several miRNAs, as well as their targets, are known to regulate stress response but systematic identification of the same is limited. The present work has been designed for in silico identification of miRNAs against a total of sixty-one DIR genes in Oryza sativa Indica followed by target prediction of identified miRNAs through the computational approach and thereafter validation of potential miRNAs in rice genotypes. We systematically identified 3 miRNA and their respective DIR specific target gene in Oryza sativa Indica. The expression of these three miRNAs and their respective DIR specific targets were validated in rice seedlings subjected to five different abiotic stress conditions (heavy metal, high temperature, low temperature, salinity and drought) by quantitative Real-Time PCR (qRT-PCR). Expression analysis indicated that miRNA under stress conditions regulates the gene expression of the DIR gene in rice. To the best of our knowledge this is this is the first report in any organism showing the expression of ath-miRf10317-akr, and osamiRf10761-akr miRNAs in response to various abiotic stresses.


Introduction
Rice (Oryza sativa L.) is one of the major cereal grain crops of the developing countries and the staple food of around 78% of the world's population. Productivity of rice is greatly affected due to various abiotic stresses which includes salinity, drought, heavy metal and extreme temperature condition. In response to various abiotic and biotic stresses, plants synthesize several proteins including the dirigent proteins (DIRs) [1].
Dirigent proteins are extracellular glycoproteins with high β-strand content and have been found in all vascular plants, including lichens, ferns, gymnosperms and angiosperms [2][3][4]. DIRs involved in secondary metabolism, lignan and lignin biosynthesis [5,6]. DIRs were found to mediate regio-and stereoselectivity of bimolecular phenoxy radical coupling during lignin biosynthesis [7][8][9]. Exposure to abiotic stress lead to the modulation of lignifications levels which is an implication of DIRs and peroxidases [1]. The expression of the most responsive DIRs was found to be correlated with increased lignifications https://doi.org/10.1016/j.ncrna.2020.09.002 Received 16 July 2020; Received in revised form 13 September 2020; Accepted 15 September 2020 when studied for cold stress [10,11]. In soybean roots, the peroxidase activity and wall lignifications get enhanced under stress due to manganese (Mn) toxicity [12]. Application of drought, salt, and oxidative stress resulted in the exhibition of stem specific expression by a Saccharum spp. DIR gene (ScDIR) [13]. In Medicago sativa, transcriptional up-regulation of one DIR gene was observed under heat stress whereas transcriptional down-regulation of two peroxidases and another DIR was observed under cold stress [14].
In plants, many studies have revealed that microRNAs (miRNAs) play a vital post-transcriptional regulatory role in gene expression by target mRNA cleavage or translational inhibition [15]. In plants, mature miRNAs are generated from the long stem-loop precursor (pre-miRNAs) by a DICER-like RNA endonuclease and then the RNA-Induced Silencing Complex (RISC) guided by ARGONAUTE 1 (AGO1) protein directs the miRNA to the complementary target mRNA sequence [15][16][17][18]. Plant miRNAs are reported to possess important functions in several metabolic and biological pathways such as tissue development and differentiation, biotic and abiotic stress responses, phytohormones signaling, and secondary metabolite production [19,20]. Nonetheless, the evolutionary highly conserved nature of an extensive number of miRNAs simplified the process of characterization of novel miRNA orthologs in new plant species through homologs identification [21]. Several abiotic stress-sensitive miRNAs have been reported over a period of time in various studies, e.g. In Arabidopsis thaliana, miRNA398 is involved in oxidative stress tolerance [22], and gene expression of 21 miRNAs are up-regulated in response to UV-B exposure [23].
Knowing the importance of miRNA and their roles in gene regulation, in the present investigation, an experiment has been designed for in silico identification of miRNAs and their potential DIR targets in rice through computational approach and validation of putative miRNAs using quantitative real-time PCR (qRT-PCR) under different abiotic stress condition.

Identification and domain analysis of DIR family genes in rice
We identified candidate DIR family genes by using their respective Pfam ID i.e., "PF03018" against the rice genome database in the Ensembl Plant (http://plants.ensembl.org/Oryza_indica/Info/Index). The amino acids, as well as the cDNA sequence of all the selected DIR proteins, were then retrieved from the Ensembl Plant database for further analysis. The amino acid sequences obtained were used for domain analysis using MEME Suite [24], Pfam [25] and NCBI's Conserved Domains Database [26]. Multiple sequence alignment was performed in MEGA X program applying MUSCLE algorithm using default parameters [27]. The aligned DIR protein sequences were used for construction of phylogenetic trees using default parameters in MEGA X program applying neighbor joining algorithm [28].

Identification of potential miRNAs and their target DIR gene
Workflow of the identification and characterization of potential miRNAs, and target genes is depicted in Fig. 1. A total of 10898 mature miRNA sequences were retrieved from PMRD: Plant micro RNA Database (http://bioinformatics.cau.edu.cn/PMRD/) [29]. With identity value 90, CD-HIT-v4.5.4 was used to remove the redundancy in miRNA sequences [30]. In order to identify miRNA-targeted DIR genes of Indica rice, Local BLAST was performed using Blast2GO version 5.2 [31]. BLASTx analysis (E-value ≤ 1e −10 ) was performed in order to remove protein-coding sequences from precursor sequences.

Prediction of the secondary structure of pre-miRNAs
Prediction of the secondary structure was done by using the software MFOLD 3.1 [32] available at (http://www.bioinfo.rpi.edu/ applications/mfold/rna/form1.cgi). The following criteria were used for screening the candidates of potential miRNAs: minimum length of the pre-miRNA to be 60 nt; pre-miRNA should be folded into appropriate stem-loop hairpin secondary structure; mature miRNA sequence should be placed in one arm of the hairpin structure; not > 6 nt mismatches in miRNA/miRNA duplex; No loops or breaks between the miRNA/miRNA duplex; A+U content within 30-70%; Predicted secondary structure should have higher minimal folding free energy index (MFEI) and negative minimal folding free energy (MFE) values [33]. The MFE or ΔG (-kcal/mol) values generated from the MFOLD web server of the stem-loop structures were used for calculating the MFE index values using the following formula: (MFE/length of precursor miRNA sequence) 100 % GC content

Plant growth and stress treatment
Seeds of Oryza sativa were surface sterilized with 70% ethanol for 1 min followed with 0.1% Mercuric chloride (HgCl 2 ) for 5 min and then 0.2% Bavastin for 10 min. All seeds were placed in the dark for 2 days then allowed to germinate for 15 days under control condition maintained at 26 ± 2 °C with 16/8 h light/dark photoperiods cycle. Thereafter the seedlings were subjected to different abiotic stress for 48 h, which includes high temperature (48 °C), low temperature (4 °C), heavy metal (6 mM CdCl 2 ·H 2 O), salinity (200 mM NaCl) and drought (15% polyethylene glycol). Total RNA was isolated from various abiotic stress treated as well as untreated seedlings and were used as a template for the cDNA synthesis following the manufacturer's protocol (Fermentas, EU). Two sets of cDNA were synthesized and used for expression analysis. For DIR gene expression analysis, 2 μg total RNA was used for first-strand cDNA synthesis. However, for miRNA expression analysis, the total RNA was initially polyadenylated with Poly (A) Polymerase Tailing Kit following the manufacturer's protocol (Epicentre, USA; Cat. No. PAP5104H). The polyadenylated RNA was then used for the synthesis of cDNA using a reverse primer of specific miRNA. The cDNA synthesized was later used as a template for expression studies using quantitative real-time PCR (qRT-PCR).

Validation using quantitative real-time PCR
Gene-specific primers for all miRNAs and their DIR specific target gene were designed using miRPrimer [34] and OligoAnalyzer Tool 3.1 tool (Integrated DNA Technologies, Inc) respectively. The rice actin 1 gene was used as an internal control to normalize the gene expression level. The qRT-PCR was performed on an AriaMx Real-Time PCR System (Agilent Technologies). The total reaction volume was 10 μl which contained 5 μl of 2X KAPA SYBR FAST qPCR Master Mix Universal, 200 nM gene-specific primers and 0.5 μl of cDNA. The thermal cycle reaction conditions were 95 °C for 3 min, followed by 40 cycles of 95 °C at 10 s and then 57 °C for 30 s. A melting curve was generated at the end of 40 cycles for analyzing the specificity of each gene. The experiment was conducted with two independent biological replicates and three technical replicates for each sample. The relative gene expression of the individual gene was calculated via 2 ∧ −ΔΔCT [35].

Statistical analysis
All the experimental data are means of triplicates and represented as means ± standard deviation (SD). The significance was tested using SPSS (Statistical Package for the Social Sciences) software (version 21 for Windows; IBM Ltd., Japan) for calculating Students's t-test at significance level p ≤ 0.05.

Identification and domain analysis of DIR family genes in rice
In the present study, a total of 61 potential DIR family genes in rice were identified from Ensemble Plant database and thereafter their amino acid sequences were retrieved for further analysis ( Table 1). The amino acid sequence of all theses 61 proteins showed the presence of the DIR domain (Accession No: PF03018) ( Table 2). These 61 DIR proteins of Indica were aligned with 49 DIR proteins of Japonica rice as reported by Liao et al., 2016 [36]. The sequence homology was evident between the DIR proteins of Indica and Japonica rice as observed through the phylogenetic tree. The phylogenetic classification of DIR proteins revealed that the DIR proteins of Indica and Japonica rice can be divided into six major groups and in some cases into further 2-3 subgroups with characteristic motifs (Fig. 2).

Identification of miRNAs and their DIR specific target genes in rice
For the in silico prediction of potential rice miRNAs, a reference set of 10898 mature plant miRNAs was retrieved from the PMRD database (http://bioinformatics.cau.edu.cn/PMRD/). With identity value 90, CD-HIT-v4.5.4 was used to remove the redundancy in miRNA sequences. After removing redundant sequences, a set of 5025 miRNA sequences (reference set of miRNA sequences) were blasted (Local BLAST by using Blast2GO-v5.2) to the DIR genes assembly (Fig. 1).
Further, BLASTx analysis (E-value ≤1e −10 ) showed that out of 6 miRNA only 3 miRNA sequences were found to be non-coding indicating that they may play the role in Oryza sativa while the other 3 are coding for some protein. The putative miRNAs obtained had lengths of 21, 21 and 27 nucleotides for ath-miRf10317-akr, cre-miR910, and osa-miRf10761-akr, respectively (Table 3). Further, in order to identify miRNA-targeted DIR genes of rice, Local BLAST was performed using Blast2GO. The potential DIR specific targets of ath-miRf10317-akr, cre-mir910 and osa-miRf10761-akr were found to be BGIOSGA034397, BGIOSGA024969 and BGIOSGA036979, respectively (Table 3).

Prediction of the secondary structure of potential miRNAs
The three non-coding miRNA sequences i.e., ath-mirf10317-akr, cre-mir910 and osa-mirf10761-akr, were further used for secondary structure analysis including hairpin stem-loop structure using MFOLD version 3.1. The negative MFE (-ΔG) of the miRNA precursors were also generated to study the stability of the hairpin stem-loop structure ( Table 4). In comparison to the length of miRNAs, the length of putative precursor miRNAs of rice also varied significantly. It was observed as 227 nt, 129 nt and 147 nt for cre-miR910, osa-miRf10761-akr, and ath-miRf10317-akr, respectively ( Fig. 3A-C). The stability of the secondary hairpin structure of pre-miRNA was determined by MFE (-ΔG). It was 211.20, 62.10 and 39.60 kcal/mol at 37 °C for cre-miR910, osa-miRf10761-akr and ath-miRf10317-akr, respectively. The distribution of G, C, A, and U nucleotides in the pre-miRNA were found to be different, where it ranged from 34.88 to 11.89% for A, 37.21-12.77% for U, 37.44-12.77% for G and 37.89-13.8% for C, respectively (Table 5). In the present investigation, AU content of miRNA cre-miR910 was found to be 25% which is below the 30-70% set range. Thus, miRNA cre-miR910 failed to qualify one of the eight criteria that we used to identify potential miRNA.

Seed germination and stress treatment
Rice seedlings were germinated for 15 days followed by 48 h abiotic stress treatment which includes salinity, drought, heavy metal, high and low temperature. After 48 h, the seedlings which were exposed to stress treatment showed significant morphological changes (Fig. 4). Seedlings showed chlorosis of the shoot region, curling of leaves, necrosis in shoot and root tissue, and stunted growth of both shoot and root as compared to the untreated control seedlings.

Quantitative RT-PCR analysis
The expression of miRNA and their DIR specific target genes under different abiotic stress conditions was analyzed using quantitative Real Time PCR. Expression of osa-miRf10761-akr miRNA and it targets DIR gene showed inverse expression pattern. Expression of osa-miRf10761akr miRNA showed significant down-regulation (p ≤ 0.05) whereas its target DIR gene (BGIOSGA036979) showed significant up-regulation in response to drought, low temperature and salt stress. Under high temperature stress, expression of DIR gene (BGIOSGA036979) got significantly down-regulated whereas osa-miRf10761-akr miRNA showed significant up-regulation. However, the expression in response to heavy metal treatment was insignificant (Fig. 5A). Except heavy metal, expression of ath-miRf10317-akr miRNA showed significant up-regulation in response to drought, salinity, high and low temperature stress whereas it target DIR gene (BGIOSGA034397) showed significant down-regulation. Under heavy metal stress, the DIR gene and ath-miRf10317-akr miRNA showed inverse pattern (Fig. 5B). Expression of cre-miR910 miRNA showed significant down-regulation and its target DIR gene BGIOSGA024969 remained up-regulated in response to heavy metal, low temperature and salt treatment. Under high temperature and drought treatment, the expression of cre-miR910 miRNA was up-regulated while its target DIR gene was observed as down-regulated (Fig. 5C).

Discussion
Rice (Oryza sativa L.) is a model plant species, which is ranked second after maize in production globally. It is mostly cultivated in Asian countries like China, India, Indonesia, Bangladesh, Vietnam and Thailand. Rice is sensitive to various abiotic stresses which includes drought, salinity, heavy metals, high and low temperature. Rice in response to different abiotic stress undergoes several changes at the morphological, physiological, biochemical and molecular level [37,38]. Upon exposure to abiotic stress several genes/proteins gets up-or down regulated [39][40][41]. These proteins are believed to play vital role against biotic and abiotic stresses in plants. One such protein is the dirigent protein (DIRs) which is involved in lignifications. Lignin is mainly deposited in the vascular tissues during plant development and provides additional strength and protection to the cell wall. Liao et al., 2016 using the NCBI database reported genome-wide analysis of 49 DIR or DIR-likes genes in rice (Oryza sativa Japonica Group) [36]. However, in the present study, we have used Ensemble Plant database for the identification of DIR family genes in Indica rice. Using Pfam ID "PF03018" as keyword, we have identified a total of 61 DIR family genes in Indica rice. All the DIR proteins showed the presence of conserved dirigent domain (Accession No: PF03018) which is the characteristic feature of DIR family protein [2]. In pepper (Capsicum annuum L.), Pfam ID "PF03018" was used as keyword for the identification of dirigent gene family [42]. DIR proteins were also reported in numerous plant species, including lichens, ferns, gymnosperms, and angiosperms [1][2][3]. The phylogenetic classification along with the identified motifs, as shown in Fig. 2, distinctively identified six major groups of DIR proteins in both Indica and Japonica group of Oryza sativa named as DIR Group I to VI. All the groups showed the presence of 6 unique motifs except Group II and Group VIB which excluded the motif "MNLVFTDGPYNGSTL" marked in yellow. Some unique features could also be observed in both Indica as well as Japonica were there was loss of the characteristic motif(s) in the individual group members like BGIOSGA036979, BGIOSGA037221, BGIOSGA008145, BGIOSGA026223, LOC_Os07g44380.1, BGIOSGA023849, LOC_Os07g44260.1, BGIOSGA000 530, BGIOSGA031974, LOC_Os10g25900.1, BGIOSGA004916, BGIOS GA004917, BGIOSGA021931, LOC_Os11g42500.1 and double occurance of all the 6 motifs in LOC_Os07g01620.1. The results in Table 2 where four dirigent superfamily members in Indica showed absence of the unique dirigent motif (PF03018) could also be corroborated with Fig. 2 where the dirigent motif "DIAAEVRELSVVGGTGKFRMARGYALLRT" marked in red was found absent in BGIOSGA036979, BGIOSGA037221, BGIO SGA008145 and BGIOSGA004917.
Earlier studies have shown that dirigent genes were expressed in different patterns in response to various abiotic stresses. In Japonica rice, expression of 13 OsDIR genes in response to dehydration, salinity and cold stresses were analyzed in rice seedling [36]. Expression analyses showed up or down regulation of OsDIR genes indicated that OsDIR genes are involved in the response process of abiotic stresses. In the present study, expression of three DIR genes i.e., BGIOSGA036979, BGIOSGA034397 and BGIOSGA024969 were analyzed in response to salinity, drought, high and low temperature and cadmium stresses in rice seedling. Differential expression showed low or high expression of DIR genes in response to various abiotic stresses. Expression of all the three DIR genes showed more than two-fold increase under cadmium, salinity and drought stress. Similar observation with two-fold down or up expression for selected sixteen OsDIR genes was also noted under at least one of the stress conditions studied in Japonica rice [36].
We identified and characterized miRNA and their DIR specific targets in Oryza sativa var Indica using bioinformatics approaches. Using the Plant micro RNA Database, we obtained 3 miRNAs i.e., ath-mirf10317-akr, cre-mir910, osa-mirf10761-akr, and their respective DIR specific targets were BGIOSGA034397, BGIOSGA024969, BGIOSGA036979. In Chlamydomonas reinhardtii, cre-mir910 miRNA targets the NCR2 gene which gets up-regulated in stress to protect the cells from damage induced by stress [59]. In the present investigation, differential expression of ath-mirf10317-akr, cre-mir910, osa-mirf10761-akr, and their respective DIR specific targets genes in response to salinity, drought, cadmium, high and low temperature stresses in rice seedling were analyzed. In response to salinity, cadmium and low temperature stresses the expression of cre-mir910 and osa-mirf10761-akr miRNA were down-regulated compared to its respective DIR specific targets genes. Similar result was observed for the cre-mir910 miRNA where its expression was found down-regulated under multiple stress conditions in C. reinhardtii [59]. Role of miRNA ath-mirf10317-akr and osa-mirf10761-akr in various abiotic stresses is not reported till date. However, few reports showed the targets for miRNAs ath-mirf10317-akr codes for Synaptobrevin/vesicle-associated membrane protein (VAMP) (AT1G08820.2) [60] and several targets of miRNA osa-miRf10761-akr has been identified which includes inorganic H+ pyrophosphatase (LOC_Os06g08080.1), DEAD-box ATPdependent RNA helicase (LOC_Os02g54020.1), PPR repeat containing protein (LOC_Os10g33700.2) and expressed protein (LOC_Os08g38620; LOC_Os09g15639) [61]. Thus, the identification and functional elucidation of miRNA targets are crucial to uncover the roles of miRNAs under abiotic and biotic stress [62][63][64][65].

Conclusion
In the present study, we identified 3 miRNAs i.e., ath-mirf10317akr, cre-mir910, osa-mirf10761-akr and their respective DIR specific targets genes i.e., BGIOSGA034397, BGIOSGA024969 and BGIOSGA036979, respectively. This is the first report in any organism showing the expression of ath-miRf10317-akr and osamiRf10761-akr miRNAs in response to various abiotic stresses. Several miRNA are known to regulate the expression of genes that are involved in abiotic stress tolerance. Expression of miRNA and it targets DIR gene showed inverse expression pattern in response to different abiotic stress treatment. We concluded that dirigent protein, which is involved in lignifications, plays an important role in abiotic stress response in plants.   and paper preparation. All the authors have read and agreed to publish the version of the manuscript. Shourya Mehra: designed the research project and is the corresponding author, performed computational work, Formal analysis, and paper preparation. All the authors have read and agreed to publish the version of the manuscript. Sayan Chatterjee: designed the research project and is the corresponding author, performed computational work, Formal analysis, and paper preparation.

Funding
All the authors have read and agreed to publish the version of the manuscript. Ram Singh Purty: designed the research project and is the corresponding author, performed computational work, Formal analysis, and paper preparation. All the authors have read and agreed to publish the version of the manuscript.

Declaration of competing interest
The authors declare that there are no conflicts of interest.