Major Domestication-Related Phenotypes in Indica Rice are Due to Loss of miRNA-Mediated Laccase Silencing

Domestication of rice ( Oryza sativa ) included conversion of perennial wild species with few seeds to short plants that produced abundant seeds. Most domestication-associated changes were due to variations in transcription factors and other key proteins such as enzymes. Here, we show that multiple yield-related traits associated with indica rice domestication are linked to micro (mi) RNA-mediated regulation. Analysis of small (s) RNA datasets from cultivated indica rice lines, a few landraces, and two wild relatives of rice revealed the presence of abundant 22-nucleotide (nt) reads in wild relatives that mapped to miR397 precursors. miR397 was expressed at very high levels in wild relatives and at negligible levels in high-yielding cultivated lines. In its genera-specific form of 22-nt, miR397 targeted mRNAs encoding laccases that decayed and induced robust secondary cascade silencing in wild species that required RNA-dependent RNA polymerase 6. In wild species of rice, reduced expression of laccases resulted in low lignification. As expected, over-expression of miR397 induced de-domestication phenotypes. At least 26 uncharacterized QTLs previously implicated in rice yield overlapped with laccases and miR397 genes. These results suggest that miRNAs contribute to rice domestication-associated phenotypes. (50 mM Tris-HCl, pH 8, 1.5 M CaCl 2 and 1 M NaCl) for 2 h at 4 o C on a rocker. The samples were centrifuged at 5000 rpm for 30 min at 4 o C and supernatant was collected and subjected to dialysis in 100 mM sodium acetate buffer (pH 4.5) overnight at 4 o C. Laccase purification was performed using concanavalin-A (Sigma-Aldrich) beads at 4 o C. The column was washed for 4 column volumes with wash buffer (1 M NaCl, 5 mM MgCl 2 , 5 mM MnCl 2 , 5 mM CaCl 2 ) and the column was regenerated with regeneration buffers (0.1 M Tris-HCl buffer pH 8.5 with 0.5 M NaCl and 0.1 M sodium acetate buffer pH 4.5 with 1 M NaCl) passed through the column alternatively 4-5 times (started and ended with Tris-HCl buffer). The column was equilibrated for 4 column volumes with equilibration buffer (50 mM Tris-HCl buffer pH 7.5 with 0.5 M NaCl, 1 mM MgCl 2 , 1 mM CaCl 2 and 1 mM MnCl 2 ). The dialyzed sample was passed through the column and following 4 washes in equilibration buffer to wash out the unbound protein, the bound protein was eluted with 50 mM mannose. The presence of laccases was confirmed through LC-MS analysis. In-gel staining of laccases was performed for 2 h, in a 10% (w/v) native gel with separated laccases using 5 mM 4-Hydroxyindole (Sigma, a substrate for laccase).


Introduction
Indica rice was domesticated from Oryza nivara and O. rufipogon in South and East Asia (Sweeney and McCouch, 2007;Wing et al., 2018). Both these wild relatives are still found in many rice-growing areas, often as wild plants next to paddy fields and at the edges of lakes (Sweeney and McCouch, 2007). Some of these wild accessions are perennial. The rice domestication process is reflected in the early history of artificial selection of key agronomically important traits, such as seed shattering, plant architecture, yield, stress tolerance, and grain color (Izawa et al., 2009). In the Oryza complex, four different lines of domestication have been observed. In all these independent domestication events, seed shattering is a common change, mapped to transcription factors now referred to as SHATTERPROOF (SH) genes. Other lines of domestication were linked to other phenotypes such as plant height, panicle structure and grain yield (Sweeney and McCouch, 2007;Wing et al., 2018).
Domestication was accompanied by other changes in genomes in addition to the SH locus, mediated mostly through changes in either sequences or expression levels of transcription factors and their cofactors. Several other genes related mostly to japonica rice domestication, including qSH1, PLOG1, qSW5, GIF1, and Rc, have been cloned and functionally characterized (Kovach et al., 2007;Izawa, 2008). In such instances, simple changes in the promoter sequences of these genes altered their expression (Huo et al., 2017), while in some cases insertion of a transposon was responsible (Li et al., 2017). There were also SNPs that likely alter the key amino acids in important proteins such as enzymes (Kharabian-Masouleh et al., 2012;Ma et al., 2015).
However, it is well established that phenotypic diversity of domesticated crops is not proportional to genomic changes (Meyer and Purugganan, 2013). This has prompted several lines of enquiry into possibilities of other causes that might have contributed to phenotypic diversity and variations associated with domestication. Since the genetic basis for many QTLs that influence the majority of yield-related phenotypes in rice are still unknown, it was speculated that small RNA and epigenetic variations might contribute towards domesticationassociated phenotypes (Campo et al., 2013;Qin et al., 2014). Recently, it has been found that heritable phenotypes of hybrids are likely due to epigenetics and small regulatory RNAs (sRNAs) (Shivaprasad et al., 2012b;Bond and Baulcombe, 2014). sRNAs of the micro(mi)RNA class have been implicated in heritable changes in diverse crops (Li et al., 2012;Shivaprasad et al., 2012b;Groszmann et al., 2011;Houston et al., 2013). Epigenetic regulation was also implicated in anthocyanin biosynthetic pathways in maize (Zea mays) (Hollick and Chandler, 2001). Several epialleles have strong genetic linkage to changes in phenotypes of crops and the DNA methylation status in the regulatory sequences of the genes (Niederhuth and Schmitz, 2014; Manning et al., 2006). sRNAs are key regulators of gene expression across eukaryotes. sRNAs can either be processed from completely complementary dsRNA molecules derived from endogenous or exogenous sources, or from structured RNAs that fold back into dsRNA-like structures (Baulcombe, 2004). In the former case, host RNA-dependent RNA polymerases (RDR) play a crucial role in synthesizing completely complementary dsRNA. The resulting dsRNA molecule is cleaved by the action of multiple Dicer-like proteins (DCL) in plants to generate 20-22-nt miRNAs or small interfering (si)RNAs of 21-24-nt length. The fate of these sRNAs depends on their association with specific Argonaute (AGO) proteins (Baulcombe, 2004).
For example, miRNAs associate with AGO1 to induce targeting of mRNA that shares high complementarity (Baumberger and Baulcombe, 2005). The degree of targeting depends on the nature of the targeted mRNA as well as the extent of complementarity with miRNA. Usually, the target mRNA is cleaved to result in reduced protein accumulation if there is near-complete complementarity. There are ample examples of plant miRNAs participating in feedback loops and subtle spatio-temporal regulations (Willmann and Poethig, 2007). However, most striking among the functions of plant miRNAs is their ability to initiate regulatory cascades to target a large number of mRNAs. This requires conversion of the miRNA-targeted mRNA fragment into a dsRNA molecule by RDR6 (Howell et al., 2007). Most miRNA-mediated targeting does not result in RDR6-dependent biogenesis of secondary siRNAs. Among the reasons known for miRNA-induced RDR6-dependent biogenesis of secondary siRNAs (Manavella et al., 2012;Felippes and Weigel, 2009), length of miRNA, i.e., 22-nt instead of the common 21-nt, is a strong determinant. Once initiated, the cascade silencing generates abundant secondary siRNAs of usually 21-nt in length (Johnson et al., 2009). These siRNAs further associate with specific AGO proteins. The most conserved secondary siRNA loci among plants are derived from the TAS3 locus (Chen et al., 2007) that regulates auxin response factor (ARF) mRNAs, playing crucial roles in development by regulating auxin signaling. The best characterized secondary siRNAs are known as transacting siRNAs (tasi-RNAs) or phased siRNAs (phasi-RNAs), and their targets include a pentatricopeptide and ARF gene mRNAs (Axtell et al., 2006).
Although it is not entirely clear why cascade silencing targets are targeted by multiple small RNA species, it is possible to speculate that the secondary siRNAs play a role in either reinforcing or coordinating the action of primary sRNAs or miRNAs.
In this study, using hundreds of indica rice lines, we explored if domestication-related traits are associated with sRNA-mediated regulation. We profiled sRNAs from two wild species of cultivated indica rice, a few landraces and a few high yielding cultivated lines, selected based on maximum phenotypic variability. We observed an unusual abundance of 22-nt sRNAs in wild species that was mapped to miR397 precursor in chromosome (Chr.) 2. miR397 expression in wild species resulted in selective degradation of laccases, thereby reducing lignification in stem tissues. Complementation of high expressing miR397 variant in cultivated lines resulted in phenotypes that match de-domestication. Surprisingly, miR397 was previously implicated in yield when it was overexpressed in japonica rice through regulation of the brassinosteroid pathway (Zhang et al., 2013). Our results suggest the presence of a previously unknown regulator of indica rice domestication.

Wild relatives of cultivated rice have an unusual abundance of 22-nt sRNAs
To explore if domestication-related traits are associated with small RNA-mediated regulation, we initially performed phenotypic analysis of 290 indica rice landraces grown in the Indian subcontinent along with two wild species of rice that are naturally grown (O. nivara and O. rufipogon) and 12 high-yielding cultivated rice lines that are the results of modern breeding.
These lines exhibited tremendous variation in phenotypes. Among these lines, we selected two high-yielding lines (O. sativa indica Pusa basmati 1, henceforth PB-1, and whiteponni, hereafter WP), six landraces with diverse phenotypes, along with the two wild species of rice for molecular analysis based on maximum phenotypic variations. O. nivara is an annual diploid grass with an AA genome and O. rufipogon is a perennial diploid grass also with an AA genome (Sweeney and McCouch, 2007). Landraces chosen for analysis generally had variations in vegetative growth, development and panicles, phenotypes strongly linked to domesticationassociated phenotypes (Supplemental Fig. 1A).
We profiled sRNAs from three different tissues from each of the six rice lines in biological duplicates (Supplemental Table 1). All vegetative tissues had a strong peak for 21-nt, and this class in plants generally includes miRNAs and siRNAs. There was also a strong peak for 24-nt sRNAs (Supplemental Figs. 1B and 1C). In flagleaf tissues, as expected for a reproductive tissue (Chávez , there was a higher abundance of 24-nt reads (Supplemental Fig. 1B). Most unusually, profiles of sRNAs from O. nivara, and to some extent, O. rufipogon, showed higher accumulation of 22-nt sRNA reads ( Fig. 1A and Supplemental Fig.   1D). This is unusual for multiple reasons. Most plant sRNA profiles lack a strong peak for 22-nt reads since this is not a major sRNA size class (Chávez . Furthermore, 22nt reads are generally produced from a DCL that is not very active (DCL2) during normal plant development. In addition, a higher percentage of 22-nt reads was prominently observed only in flagleaves, indicating tissue-specific expression of these sRNAs. Most importantly, this unusual abundance was restricted to wild species of rice and was not observed among other rice lines.
Among the 22-nt reads, there was a strong bias for 5' U in datasets derived from wild species of rice ( Fig. 1B and Supplemental Fig. 1E), indicating that the majority of these 22-nt sRNAs were likely associated with AGO1 or its close homologs.

Mapping of an abundant 22-nt sRNA-producing locus
We performed differential expression analysis of sRNA loci between rice lines to investigate the unusual abundance of 22-nt sRNAs between the rice lines. This analysis identified 358 differentially expressing sRNA loci (Supplemental Dataset 1). These loci were mapped to all 12 chromosomes of rice. Differentially expressing loci were mapped to many short and long genomic regions; however, one specific locus on chromosome 2, between SSR marker RM12499 and RFLP marker HHUK22, accounted for the most difference between O. nivara and the other rice lines (Figs. 1C and 1D). This region alone accounted for nearly 120 thousand reads per million (RPM) of difference between O. nivara and other rice lines such as PB-1 (Supplemental Fig. 1F).
The mapped region of chromosome 2 encoded a single non-coding RNA expressing a miRNA precursor. This precursor was the single source for the abundant 22-nt reads in this region of chromosome 2. This precursor generates miR397 in both dicots and monocots (Supplemental miR397 is usually 21-nt in length in most plants. The mature miR397 identified here has, in addition to a shift of 3 nucleotides, an extension by 1-nt when compared to other miR397 precursors of other miRBase entries (Supplemental Fig. 2B). The 22-nt miR397 is incorporated into the AGO1 isoform (Wu et al., 2009). We also identified a 21-nt version of miR397 in all our rice datasets; however, accumulation of the 21-nt form is at least 20-fold lower than the 22-nt form. miR397 expressed in O. nivara is 22-nt in length (Supplemental Fig. 2C). To confirm the high accumulation of miR397 in O. nivara, we performed an RNA blot analysis for three different rice lines (Fig. 1F). The tissues derived from O. nivara accumulated very high levels of the 22-nt form of miR397.
Since the presence of 22-nt miR397 has been reported only from rice and maize (Jeong et al., 2011), we performed a miRProf analysis to assess the diversity of miR397 sequences in different plant species (Shivaprasad et al., 2012a). The 22-nt version was abundant only among members of the Oryza genera (Supplemental Fig. 2D). The presence of genera-specific isoforms indicates that this miRNA has an evolutionarily recent origin. We also used miRProf analysis to confirm differential expression of miRNAs across wild species and other lines (Supplemental Fig. 3, Supplemental Dataset 2). Few miRNAs had differential accumulation between rice lines irrespective of the tissues, among which miR397 exhibited maximum variation between lines across three tissues analyzed. Other miRNAs that had variable expression included miR408, miR398 and miR528 across three tissues, and miR1568 and miR168 in flagleaf and panicle tissues. Their targets include genes involved in development and various stress responses. Considering the spectrum of modifications that were introduced through domestication, it is possible that these miRNAs might have contributed to important phenotypes.

miR397 targets rice laccases
In Arabidopsis, maize, poplar (Populus) and Medicago (Medicago truncatula), miR397 targets a family of laccases Lu et al., 2013). Laccases are also the major targets of miR397 in rice, even though the rice miR397 isoform had an extra nucleotide, this made no difference to its targeting ability. Rice has 30 laccases encoded in its genome. To identify the ability of miR397 in targeting rice laccases and other possible previously unidentified targets, we performed psRNA target (Dai and Zhao, 2011) prediction and tapir (Bonnet et al., 2010) tool analysis. We identified 15 laccases as potential targets of miR397 under stringent conditions (tapir score >4), and when the score was relaxed, the number of potentially targeted laccases increased to 28 ( Fig. 2A). These miRNA targets belonged to different clades of laccases (Supplemental Fig. 4), confirming that laccases of diverse sequences were likely targeted by the miRNA. We predicted that miR397 might be able to target multiple clades of laccases due to its ability to target a specific RNA motif encoding a specific domain. Degradome analysis derived from pooled tissues of O. nivara confirmed targeting of at least six laccases by miR397.
However, in a publicly available O. rufipogon degradome dataset (Wang et al., 2012), reads of OsLAC7 but not of other laccases mapped at the expected slicing site. Analysis of target regions of these laccases identified an RNA motif targeted by miR397 (Fig. 2B). This specific RNA motif codes for the INLAAN peptide ( Fig. 2B). Such a peptide sequence is observed only among laccases in plants. Laccases in other plants were also targeted by miR397 family members in the same region.

miR397 induces RDR6-dependant cascade secondary silencing in target loci
Generation of secondary siRNAs from miRNA-targeted RNAs is unusual among plants. These secondary siRNAs, also known as tasi-RNAs or phasi-RNAs generate a cascade silencing effect in miRNA target regions. Since miR397 is 22-nt long, we hypothesized that miRNAtargeted laccase regions might be sources of such secondary siRNAs. We performed tasi-RNA analysis prediction to look for rice regions that have signatures of secondary siRNA generation.
As expected, the well-conserved TAS3 locus indeed produced abundant secondary siRNAs in all lines at similar levels (Table 1). Higher numbers of TAS3-derived siRNAs were observed in seedling tissues. However, there were also eight laccase loci acting as PHAS loci in O. nivara, but not in other lines including those from cultivated rice lines such as PB-1 (Fig. 2C, Table 1).
The presence of secondary siRNAs derived from these laccases, in phase with the miR397 target region, is direct evidence that miR397 targets laccases in O. nivara. This is significant because the presence of phasi-RNAs derived from laccases have not been reported previously, although targeting of laccases by miR397 has been well documented in multiple species where miR397 is expressed in the 21-, but not 22-nt, form. Laccase loci that produced abundant secondary siRNAs in O. nivara included LAC3, LAC7, LAC12 and LAC13 (Figs. 2D and 2E and Supplemental Figs. 5A to 5C). miRNAs of 22-nt initiate secondary silencing in plants that require RDR6 to generate dsRNA substrates from miRNA target regions (Cuperus et al., 2010). To check if laccase-derived secondary siRNAs are indeed dependent on RDR6, we introduced laccase from rice into tobacco plants silenced in RDR6 (Schwach, 2005). Tobacco miR397 abundance was extremely low in publicly available tobacco sRNA datasets as well in our RNA blot analysis (Supplemental Fig. 5D). In WT and RDR6i plants, introduced laccases and miRNAs were expressed at high levels. Secondary siRNAs (D10-) were abundantly expressed in WT tobacco plants, but not in RDR6i plants (Supplemental Figs. 5D and 5E). These results indicate that laccase-derived siRNAs required RDR6 to generate dsRNA intermediates that can be cleaved by DCL4. Multiple secondary siRNAs derived from laccases potentially target parent and other laccases, likely resulting in robust silencing (Supplemental Fig. 5F).
The evolution of the 22-nt form in Oryza genera members might have resulted in robust silencing of laccases that is reinforced due to the formation of a dsRNA and subsequent production of secondary cascade siRNAs. This is confirmed with the requirement of RDR6, an enzyme that is generally associated with silencing signal amplification. It is however unclear why this reinforced regulation is required in Oryza genera members where miR397 is a 22-nt miRNA, but not in other plants where miR397 is expressed in a 21-nt form. One possible explanation for the proposed mechanism could be that the precursor of miR397 in Oryza speciation accumulated genetic changes that altered the precursor structure to accommodate robust silencing of laccases (Supplemental Fig. 2A).

Differential expression of laccases among rice lines
Domestication-associated phenotypes are often complex, involving massive changes in genomes and expression of their contents. To identify transcripts that accumulate differently between rice lines, as well as to identify the outcome of miR397 targeting on laccases and other pathway genes, we performed a transcriptome analysis using two tissues of rice lines with variable expression of miR397. In addition to cultivated rice lines (WP and PB-1) expressing miR397 at lower levels, and wild species that express miR397 at high levels, we also used landraces that have intermediate levels of miR397 (Fig. 3A). Up to 90% of all reads matched to the rice genome, almost to similar levels between wild species, landraces and cultivated high yielding lines (Supplemental Table 2). Many transcription factors, disease/pathogen resistance and starch biosynthesis genes had variable expression between rice lines (Figs. 3B and 3C and Supplemental Datasets 3 and 4).
Among laccases, we observed expression of only 16 laccases out of 30 in rice tissues analyzed in this study. Flagleaves had expression of 12 laccases, while panicle tissues expressed 16 laccases at variable levels ( Fig. 3D, Supplemental Fig. 6 and Supplemental Dataset 5).
Strikingly, almost all laccases were expressed at different levels between rice lines and were among the most differentially expressed genes between wild and cultivated rice lines. miR397targeted laccase such as LAC7 was abundantly expressed in cultivated rice lines such as PB-1, while it accumulated to negligible levels in O. nivara (Fig. 3E). These results, in addition to strongly confirming differential targeting of laccases that depends on expression levels of miR397, also indicate natural variation in the expression of laccases between rice landraces, further suggesting their roles in domestication-associated phenotypes.

Wild rice species accumulate reduced levels of lignin in stem tissues
Since laccase RNAs had lower expression in O. nivara due to miR397 targeting, we hypothesized that the lignification process in this species might be affected. We performed laccase assays for total tissues and observed that cultivated rice lines such as PB-1 and WP had very high laccase activity when compared to wild species (Fig. 3F). We performed phloroglucinol staining in leaf and stem tissues. While there were many morphological differences between leaf samples, such as the vascular bundle density and arrangement, there was little difference in lignin staining (Supplemental Figs.7A and 7B). However, in stem sections we observed that wild species accumulated much lower levels of lignin when compared to landraces and PB-1 (Fig. 3G). This was also well correlated with total lignin content measured by derivatization with thioglycolic acid (Fig. 3H). Similarly, total laccase protein levels, purified using a concanavalin column, also had higher accumulation in cultivated rice line PB-1 when compared to O. nivara ( Fig. 3I and Supplemental Fig. 7C). We separated laccases on a SDS gel and identified several laccases that are expressed only in the cultivated line PB-1 and not in O. nivara. In addition, we performed in-gel staining with 4-hydroxyindole and found that wild species had low accumulation of laccases when compared to PB-1 and WP (Fig. 3J,   Supplemental Fig. 7D).
These results strongly indicate that miRNA-mediated laccase targeting reduced lignification in wild rice lines due to reduced accumulation of laccase RNAs. Together, these results indicate how robust silencing of laccases between rice wild species, landraces and high-yielding cultivated lines is linked to natural variation in expression of a single miRNA that might be the result of selection of phenotypes requiring high lignification.

Over-expression of miR397 in cultivated rice results in morphological abnormalities associated with reduced lignification
To confirm the presence of phenotypes associated with expression of one miRNA in wild species, we generated over-expression (OE) lines in cultivated rice background. We anticipated that OE of miR397 might alter the phenotype of high-yielding lines through reduced lignification and associated changes. We generated several single-copy T-DNA transformants using Agrobacterium-mediated transformation that expressed miR397 under the maize ubiquitin promoter in PB-1 and WP backgrounds. yield-related parameters is not limited to rice but was also observed among other monocots such as maize (Sun et al., 2018) and other plants (Pinçon et al., 2001). Our results strongly suggest that miR397 expression might have been a strong selective force in indica rice domestication (Fig. 4K).
If miR397:laccase interaction was a key feature of indica rice domestication, then it is possible that these genes are part of QTLs identified for rice yield. We overlapped all uncharacterized QTLs previously mapped for rice yield (from Gramene QTL database), and found that more than 26 previously uncharacterized QTLs overlap with laccases and miR397 (Fig. 4J, Supplemental   Fig. 9, Supplemental Dataset 6).

Discussion
The origin and heritability of agronomic traits during crop domestication are complex events in evolution. Domestication of rice, especially that of the japonica subtype, has been extensively studied and the outcome of such findings was duly used in various breeding and genetic modification efforts to improve agronomic traits (Zuo and Li, 2014). In this study, we identified a novel miRNA-mediated regulation of a large family of genes coding for laccase enzymes as a major sRNA variation between wild species of rice, few selected phenotypically diverse landraces and cultivated rice lines. We propose that this regulation is involved in the evolution of various agronomic traits in indica rice.
In our comparative miRNA analysis, few plants showed high abundance of miR397. Expression levels of miR397 might be a direct proxy for laccase expression between plants, since many trees that accumulate high levels of lignins have low levels of miR397 expression (Lu et al., 2013). Surprisingly, among the 60 different species analyzed in this study, wild species of rice had the highest accumulation of miR397. Wild species of rice originated in semi-aquatic ecosystems and strong mechanical support might not be a requirement for these plants. This could explain why wild species of rice had stronger miR397-mediated silencing of laccases when compared to other plants. Also, it is possible that laccases that cannot be targeted by miR397 might have evolved in some plants where mechanical strength was a prerequisite for changing environment. Alternatively, miR397-mediated regulation might have been a result of co-evolution that was altered during the artificial selection of rice. The altered expression of miR397 might have been due to the incorporation of mutations in miR397 promoter sequences since there is an absence of any change in the precursors, and indeed we do see changes in the promoter sequence between wild and cultivated rice species (Supplemental Fig. 2E).
The presence of a 22-nt miRNA form in Oryza genera is extremely interesting since it is not found among other monocots where sequence information is available. Based on the findings in other plants (Chen et al., 2010), it is possible that the evolution of the 22-nt form might have resulted in robust silencing of laccases that is reinforced due to the formation of a dsRNA and subsequent production of secondary cascade siRNAs. This is also confirmed with the requirement of RDR6, an enzyme that is generally associated with silencing signal amplification.
It is, however, hard to speculate why this reinforced regulation is required in Oryza genera members where miR397 is a 22-nt miRNA, but not other plants where miR397 is expressed in a 21-nt form. One possible explanation for the proposed mechanism could be that the precursor of miR397 in Oryza speciation accumulated a genetic change that altered the precursor structure to accommodate robust silencing of laccases (Supplemental Fig. 2A). It is wellestablished that precursor structure is a key requirement for 22-nt miRNA biogenesis (Cuperus et al., 2010).
Surprisingly, little sequence variation in nucleotides and amino acids are seen among individual laccase family members, both within and between rice species and landraces. Interestingly, there are natural variations in miR397-targeted RNA motifs among laccases that specify the extent of miR397 targeting. Due to these variations, closely related laccases might not be targets of miR397. For example, miR397 targets OsLAC13 but not OsLAC5 because of five mismatches in the target region, although OsLAC13 and OsLAC5 belong to the same clade (according to the phylogenetic analysis). OsLAC5 escaped from the control of miR397 likely after the whole genome duplication (Guo et al., 2008). Since the functions of individual laccases of rice are not fully understood, it is not possible to interpret the importance of these interactions; however, variations in the expression of individual laccases between landraces studied here indicate a functional importance for these changes.
Based on our analysis, both through sRNA sequencing, phenotyping and genome modifications, we can define functions for rice laccases that were not well-known previously.
Our observations indicate that laccases regulate lignification and the associated phenotypes  (Ranocha et al., 2002). Lignification also affected seed shattering in rice. For example, OsSH5 enhanced seed shattering in japonica rice by altering lignin deposition at the abscission zone (Yoon et al., 2014).
What are the possibilities that a similar strategy to regulate lignification and mechanical strength is also true for other crops? Among the wild species of crops, few are submerged or prostrate plants with weak stems. In S. pimpinellifolium, a wild relative of cultivated tomato, we observed a higher abundance of miR397 when compared to cultivated tomato lines (21 RPM when compared to 1 in cultivated tomato M82). It is also important to note that S. pimpinellifolium

Plant genetic material
Two wild species of rice (O. nivara and O. rufipogon), two cultivated high-yielding lines of indica rice (whiteponni and Pusa basmati-1) and, six indica landraces from South India, namely kavuni, chomala, kartha, doddabyranel, karigajavali and kagasale, were used in this study.

Plant growth
Rice plants were grown in the greenhouse under natural light conditions and 70% relative humidity

Processing of Illumina sRNA-and RNA-seq reads and expression analysis
The sRNA-Seq reads were processed for adapter removal and filtered for length range of 16-35nt using UEA small RNA Workbench Version-3 (Stocks et al., 2012). miRNA expression analysis and the identification of the loci producing phased siRNAs and sRNA generating loci were performed using miRProf, the TASI-RNA prediction tool (p-value cut-off of 0.01), and SiLoCo (sRNA length 19 to 25 and minimum abundance of 1), respectively, using the UEA small RNA Workbench. RNA-Seq reads were processed to remove adaptors and low-quality bases. The reads obtained after preprocessing were aligned to rice genome (IRGSP-1) and analyzed further as mentioned in Trapnell et al. (2012).

miRNA target prediction, sequence alignments and phased sRNA analysis
The targets of candidate miRNAs were predicted using the psRNAtarget (Dai and Zhao, 2011) tool and tapir (Bonnet et al., 2010). The sequences in the target regions were extracted and the LOGOs of miRNA target regions were generated using WebLogo (Schneider and Stephens, 1990). Phylogenetic analysis was performed using MEGA6 with the maximum likelihood method and 10000 bootstrap replicates.

RT-PCR
cDNA was prepared using Invitrogen SuperScript-III kit with 2 µg of RNA from desired samples.
Hot start RT-PCR was performed using BioRad thermo cycler. The PCR mix for 20 µl reaction

Plasmid clones and constructs
To generate over-expression constructs, cDNA was prepared from O. nivara flagleaf samples using Superscript III (Invitrogen), according to manufacturer's' protocol. Amplified cDNAs were cloned into the pBIN19 vector (Bevan, 1984) Table 3.

RNA blot analysis
The RNA blots were performed as described by Shivaprasad et al. (2012b). Briefly, about 12-µg of total RNA extracted as mentioned earlier was dried and resuspended in 8 µl loading buffer (0.10% w/v bromophenol blue, 0.10% w/v xylene cyanol in 100% de-ionized formamide), heated at 95 o C for 1 min and loaded onto a 15% denaturing polyacrylamide gel (a 19:1 ratio of acrylamide to bis-acrylamide and 8 M urea). The gel was run at 100 V for 3 h and then transferred to a Hybond-N + membrane (GE Healthcare) by electro-blotting (Bio Rad) at 10 V overnight at 4 o C.
The RNA hybridization was performed at 35 o C for 12 h using Ultra Hyb-oligo buffer (Ambion) containing appropriate radiolabeled short DNA oligo probes (Supplemental Table 3), end-labelled with 32 P-ATP (BRIT, India) by polynucleotide kinase (NEB) and purified through MicroSpin G-25 columns (GE Healthcare). The blot was washed twice with 2X SSC, 0.5% w/v SDS for 30 min at 35 o C. The signal was detected after exposure on a phosphor screen using a Molecular Imager (GE Healthcare). For repeated hybridization, the membrane was stripped and re-probed.

Transformation of rice
miR397 and OsLAC3 over-expressing transgenic rice plants were generated using the constructs PUbi:miR397 and PUbi:OsLAC3 using the A. tumefaciens-mediated transformation method described for PB-1 and WP (Sridevi et al., 2008, Sridevi et al., 2005. Transgenic plants were regenerated from hygromycin-resistant calli and maintained in a growth chamber (70% RH at

Laccase assays
Laccase assays were performed using 4-hydroxy indole as substrate. About 500 mg of tissue was ground using liquid nitrogen and the total crude protein was extracted using sodium acetate buffer (pH 4.5). The 300 µl reaction mixture contained 10 µl of crude enzyme and 100 ul 5 mM 4-hydroxy indole. The reaction was set up at 30 o C for 30 min and the absorbance was measured at 620 nm.

Lignin quantification
Cell wall components were enriched by grinding about 300 mg of rice stem tissues in liquid nitrogen. The resulting powder was suspended in 1.5 ml of methanol and stirred vigorously for 1 h. The sample was centrifuged at 12000 rpm for 5 min at 25 o C. The pellet was consecutively treated with 1.5 ml of the following solvents; methanol (twice), 1 M NaCl, 1% w/v SDS, water (twice), ethanol, chloroform:methanol (1:1) and tert-butyl methyl ether, with 15 min stirring followed by 5 min centrifugation as previously mentioned. The pellet thus obtained was dried at 80 o C for 12 h to obtain alcohol insoluble residue (AIR).
The lignin was quantified by derivatization with thioglycolic acid. To 20 mg of AIR, 1 ml of 2 M HCl and 0.2 ml of thioglycolic acid were added and incubated for 4 h at 95 o C with continuous stirring.
The samples were rapidly cooled on ice and centrifuged at 12000 rpm for 10 min at 25 o C.
Supernatant was discarded and the pellet was washed with 1.5 ml of water thrice by disturbing the pellet, followed by centrifugation and discarding the supernatant. The resulting pellet was resuspended in 1 ml of 0.5 M NaOH and was vigorously shaken overnight at 25 o C to extract the lignin thioglycolate. Following centrifugation at 12000 rpm for 10 min at 25 o C, the supernatant was collected and the pellet was resuspended in 0.5 ml of 0.5 M NaOH and centrifuged again. The supernatant obtained was pooled with the supernatant obtained in the previous step. To the combined alkali extract, 0.3 ml of concentrated HCl was added and the sample was incubated at 4 o C for 4 h to facilitate lignin thioglycolic acid (LTGA) precipitation. The samples were centrifuged at 12000 rpm for 10 min and a brown pellet was obtained. The supernatant was discarded and the pellet was dissolved in 1 ml of 0.5 M NaOH. The absorbance was measured at 340 nm and the standard curve was plotted using alkali lignin (Sigma).

Phloroglucinol staining
Stem and leaf sections were stained with phloroglucinol to visualize lignin deposited in the secondary cell walls. The staining solution used was 2% (w/v) phloroglucinol in 95% (v/v) ethanol.
The sections were immersed in the stain for 2 min. The sections were then transferred to 10 N HCl and mounted in 5 N HCl and observed under a bright field microscope (Olympus, BX43).

Laccase extraction and purification
Around 20 g of tissue was crushed in liquid nitrogen and incubated in 200 ml of extraction buffer (50 mM Tris-HCl, pH 8, 1.5 M CaCl 2 and 1 M NaCl) for 2 h at 4 o C on a rocker. The samples were centrifuged at 5000 rpm for 30 min at 4 o C and supernatant was collected and subjected to dialysis in 100 mM sodium acetate buffer (pH 4.5) overnight at 4 o C. Laccase purification was performed using concanavalin-A (Sigma-Aldrich) beads at 4 o C. The column was washed for 4 column volumes with wash buffer (1 M NaCl, 5 mM MgCl 2 , 5 mM MnCl 2 , 5 mM CaCl 2 ) and the column was regenerated with regeneration buffers (0.1 M Tris-HCl buffer pH 8.5 with 0.5 M NaCl and 0.1 M sodium acetate buffer pH 4.5 with 1 M NaCl) passed through the column alternatively 4-5 times (started and ended with Tris-HCl buffer). The column was equilibrated for 4 column volumes with equilibration buffer (50 mM Tris-HCl buffer pH 7.5 with 0.5 M NaCl, 1 mM MgCl 2 , 1 mM CaCl 2 and 1 mM MnCl 2 ). The dialyzed sample was passed through the column and following 4 washes in equilibration buffer to wash out the unbound protein, the bound protein was eluted with 50 mM mannose. The presence of laccases was confirmed through LC-MS analysis. In-gel staining of laccases was performed for 2 h, in a 10% (w/v) native gel with separated laccases using 5 mM 4-Hydroxyindole (Sigma, a substrate for laccase).

Accession Numbers
High throughput sRNA sequencing from seedling, flagleaf and panicle samples, as well as RNA-seq data has been deposited in Gene Expression Omnibus (GEO). The super series containing both small RNA-Seq (GSE111440) and RNA-Seq (GSE111438) can be accessed through GSE111472. Supplemental Figure 9. Overlapping of yield-related QTLs with laccases and miR397. Table 1. Library statistics of sRNA datasets. Table 2. Library statistics of RNA-Seq datasets. Table 3. Primers and oligonucleotide probes used in the study. Table 4. Phenotypic analysis of miR397 OE transgenic plants.