Newly identified miRNAs may contribute to aerenchyma formation in sugarcane roots

Abstract Small RNAs comprise three families of noncoding regulatory RNAs that control gene expression by blocking mRNA translation or leading to mRNA cleavage. Such post‐transcriptional negative regulation is relevant for both plant development and environmental adaptations. An important biotechnological application of miRNA identification is the discovery of regulators and effectors of cell wall degradation, which can improve/facilitate hydrolysis of cell wall polymers for second‐generation bioethanol production. The recent characterization of plant innate cell wall modifications occurring during root aerenchyma development triggered by ethylene led to the possibility of prospection for mechanisms of cell wall disassembly in sugarcane. By using next‐generation sequencing, 39 miRNAs were identified in root segments along the process of aerenchyma development. Among them, 31 miRNAs were unknown to the sugarcane miRBase repository but previously identified as produced by its relative Sorghum bicolor. Key putative targets related to signal transduction, carbohydrate metabolic process, and cell wall organization or biogenesis were among the most representative gene categories targeted by miRNA. They belong to the subclasses of genes associated with the four modules of cell wall modification in sugarcane roots: cell expansion, cell separation, hemicellulose, and cellulose hydrolysis. Thirteen miRNAs possibly related to ethylene perception and signaling were also identified. Our findings suggest that miRNAs may be involved in the regulation of cell wall degradation during aerenchyma formation. This work also points out to potential molecular tools for sugarcane improvement in the context of second‐generation biofuels.


| INTRODUC TI ON
miRNAs are a class of relevant plant development regulators (Borges & Martienssen, 2015;Chen, 2004;Guleria & Yadav, 2011). Posttranscriptional regulation is achieved through sequence-specific base-pairing of 20 to 22 nucleotides of a small RNA and a mRNA.
Dicer-like proteins are plant ribonucleases responsible for processing primary miRNAs leading to its mature functional conformation.
Within the cytoplasm, the active guide strand (i.e., the one that effectively binds to the target mRNA) is then separated from its complementary strand (named miRNA* or star miRNA). miRNA* is degraded and mature miRNA is loaded onto the argonaute protein, giving rise to the RNA-induced silencing complex (RISC). This protein complex leads the miRNA guide strand to target mRNAs due to base complementarity within its 3' UTR sequence (Li & Zhang, 2016). The mRNA can be then targeted for degradation, or its translational process can be blocked (Borges & Martienssen, 2015).
The wide range of pathways regulated by miRNA within cellular metabolism has led to the assumption that plant miRNA can serve as important tools for plant improvement (Zhang, 2015). miRNA identification involved in plant tolerance to environmental stresses and carbohydrate metabolism can be especially relevant for crops, since increased biomass, higher digestibility, and stress tolerance are desirable traits, especially for biofuel industries. Some of these crops are Panicum virgatum, Sorghum bicolor, Zea mays, and Saccharum sp.
For instance, P. virgatum plants overexpressing Z. mays miR156 displayed distinct biomass composition and presented higher saccharification yields due to increased starch content (Chuck et al., 2011;Fu et al., 2012). Data mining for miRNA regulation among cell wall-related transcripts in Sorghum bicolor has led to the identification of 10 putative targets, including laccases and cellulose synthases (Rai et al., 2016). However, the underlying miRNA regulation of cell wall metabolism in sugarcane lacks experimental evidence.
Sugarcane is an essential crop in Brazil due to its widespread use in sugar and ethanol production. Also, further expansion of production is possible without need to use land presently occupied by preserved biomes or food production as projected to 2045 (Jaiswal et al., 2017).
Sugar is obtained through culm milling followed by sucrose extraction, whereas ethanol results from sugar fermentation. Fermentation can occur using two different substrates, which defines first-generation (1G) and second-generation (2G) ethanol. The fermentation of the sucrose stored within sugarcane culm has supported Brazil as the second larger 1G producer in the world. Although available commercially, 2G ethanol production still needs further improvements, since it requires more effective depolymerization of structural polysaccharides of the cell wall (Buckeridge & De Souza, 2017).
Recently, the characterization of aerenchyma formation within sugarcane roots uncovered the potential use of innate pathways for cell wall deconstruction, facilitating the access to cell wall polymers and paving the way toward successful 2G ethanol production (Leite et al., 2017;Tavares, Souza, & Buckeridge, 2015). Lysigenous aerenchyma is a developmental process activated by the hormone ethylene. Its action leads to the opening of gas spaces within parenchymatic tissues due to programmed cell death and cell wall modifications (Takahashi, Yamauchi, Colmer, & Nakazono, 2014;Tavares et al., 2018;Yamauchi, Colmer, Pedersen, & Nakazono, 2018;Yamauchi, Shimamura, Nakazono, & Mochizukic, 2013). Immunolocalization showed arabinoxylan debranching, homogalacturonan hydrolysis from the middle lamella, and β-glucan mobilization during the formation of aerenchyma in sugarcane roots (Leite et al., 2017). The negative regulation played upon homogalacturonan hydrolysis by an ethylene response factor, RAV1, represents an important trigger to cell wall attack (Rahji et al., 2011;Tavares et al., 2019). Thus, aerenchyma formation seems to rely on cell targeting induced by ethylene and auxin balance. This is followed by cell expansion and separation, programmed cell death, and hemicellulose and cellulose hydrolysis. Each step configures a conserved set of pathways-named "modules"-shared between other endogenous cell wall degradation events (Grandis, Souza, Tavares, & Buckeridge, 2014;Tavares et al., 2015). Although cell wall changes are subtle and sugarcane pectin content is rather low (De Souza, Leite, Pattathil, Hahn, & Buckeridge, 2013), the effects of partial pectin degradation within plant tissues, especially in the middle lamella, may be quite relevant for saccharification and bioenergy production (Latarullo, Tavares, Maldonado, Leite, & Buckeridge, 2016). During aerenchyma formation in sugarcane roots, pectin degradation is thought to be a result of the attack by acetyl esterases, endopolygalacturonases, β-galactosidases, and α-arabinofuranosidases, followed by the action of β-glucan-/callose-hydrolyzing enzymes Tavares et al., 2018Tavares et al., , 2019. Concomitantly, there are modifications in arabinoxylan (by α-arabinofuranosidases), xyloglucan (by xyloglucan endotransglucosylase/hydrolase), xyloglucan-cellulose interactions (by expansins), and partial hydrolysis of cellulose . Thus, the precise control underlying this process might be a key factor to understand modulation of cell wall changes in sugarcane. One question is to what extent epigenetics might be involved in the control of aerenchyma formation in sugarcane and other grasses.
In this study, we identified 39 expressed miRNAs within sugarcane roots undergoing aerenchyma formation. Transcripts related to carbohydrate metabolic process and cell wall organization or biogenesis were predicted to be targeted by the sequenced miRNA. Among those, miRNAs expressed during aerenchyma formation and presumably target transcripts related to several cell wall polysaccharides previously detected in sugarcane (Leite et al., 2017). A more significant fraction of these transcripts is related to pectin degradation. The diversity of miRNA targeting pectin-related and ethylene regulation-related transcripts corroborates previous data showing decreasing cellular adhesion and hormone signaling as fundamental mechanisms during the onset of aerenchyma formation Leite et al., 2017;Tavares et al., 2018Tavares et al., , 2019. Our results represent an important step toward the identification of miRNAs involved in sugarcane cell wall degradation, as well as its hormonal control. It also paves the way for the development of biotechnology in the biofuel field.
Precursor sequence displaying A/U percentage outside 30%-70% range (Zhang, Pan, Cox, Cobb, & Anderson, 2006) and those with minimum fold energy index (MFEI) values below 0.7 (Zhang et al., 2006) or below the minimum free energy of folding randomization (computed by Randfold; Bonnet, Wuyts, Rouze, & Peer, 2004) were discarded. Precursor miRNA not present in all biological replicates from the same root segment was not used in the data analysis. miRNA abundance was expressed in terms of reads per million (RPM).
All primer sequences used in this study are listed in Table S1.

| Root miRNA sequencing and precursor prediction
The development of the lysigenous aerenchyma was observed along the root segments. The root cortex was intact within S1 and S2, and the aerenchyma area was noticeable in S3 by the opening of gas spaces, which increased in S4 (Figure 1). This was in accordance with previous results Leite et al., 2017;Tavares et al., 2018).
Aiming at characterizing sugarcane miRNA transcriptome and identifying candidates related to aerenchyma formation, 4 miRNA libraries were obtained from sequenced RNA extracts from sugarcane roots with developing aerenchyma. A total of 12 libraries (3 biological replicates from 4 root segments) produced 180435212 reads after Ion Proton sequencing, ranging from 73 (S1) to 31.5 (S4) million reads (Table   S2). Non-miRNA filtering (taking out rRNA, tRNA, snRNA, snoRNA, ln-cRNA, and tasiRNA) led to 119961617 total reads in all libraries.
The remaining sequences after the filtering and trimming steps showed lengths ranging from 20 to 22 nt (Table S3), which is in accordance with the expected range for miRNA (Li & Zhang, 2016).
The dominance of 21-nt length sequences was similar to the tendency observed elsewhere for sugarcane miRNA sequences from four different cultivars  but distinct from the 22-and 24-nt bias shown by other reports with the same plant (Carnavale-Bottino et al., 2013;Lin, Chen, Qin, & Lin, 2014;Thiebaut et al., 2014). Furthermore, the majority of mature miRNA presented U as the first nucleotide (Table S3), which is in accordance with previously reported data on sugarcane miRNA .

| Reference genome read mapping and miRNA prediction
The pipeline described above allowed the identification of 39 miR-NAs, from 19 different families. Five out of 39 have already been deposited as sugarcane miRNAs in miRBase (Table 1). The remaining 34 miRNAs that have been identified in S. bicolor were named on the basis of the observed homology between predicted sugarcane miRNA and those from S. bicolor available on miRBase.
F I G U R E 2 Pipeline for miRNA prediction using mirDeep2 and target identification through psRNATarget. Ion Proton raw reads were uploaded on FastQC for quality analysis, and further, rRNA, tRNA, snRNA, snoRNA, lncRNA, and tasiRNA were removed. Filtered reads were uploaded on miRDeep2 mapper module together with Sorghum bicolor as the reference genome. Collapsed and mapped reads output were used as input on miRDeep2 core together with mature and precursor miRNA sequences from S. bicolor and Saccharum sp retrieved from miRBase v. 21. Precursor, mature, and star miRNA sequences were predicted. The mature miRNA expression level was measured as reads per million if the corresponding precursor presented a minimum fold energy index (MFEI) higher than 0.7. Target prediction was performed by uploading mature miRNA sequences from previous steps and using the sugarcane EST database (SUCEST) and Arabidopsis thaliana as queries By mapping miRNAs onto S. bicolor genome, it was observed that the majority of the miRNAs are predicted to be encoded by an independent transcriptional unit (intergenic genomic location), whereas only 8 are located into gene sequences (exon, intron, or utr according to S. bicolor genomic location) ( miRNAs not yet deposited into miRBase. miR159b, comprised within EST SCAGFL3025B10.g sequence, was one of the few precursors predicted to be encoded by an exon (according to S. bicolor genomic coordinates).

| miRNA expression profiles and target prediction by using A. thaliana genome and sugarcane EST as queries
The majority of the identified mature sequences were expressed in all root segments (Figure 3), which represents more than 86% of the total miRNA identified in the present study. This indicates that miR-NAs are present during aerenchyma development. The expression level in terms of reads per million (Figure 4 and Table S4) shows four general expression patterns regarding the range of expression (blue boxes). However, there is an overall tendency of lower miRNA expression levels on S1, except for the most basal clade. Furthermore, miRNAs from the same family do not present similar expression levels.
The set of identified miRNA was used to search for predicted targets by using as queries the SUCEST database and the A. thaliana genome. Putative targets were predicted for all the 39 miRNAs identified in the present study [747 (Table S5) and 4,728 (Table S6) targets] for SUCEST and A. thaliana genome, respectively. In most cases, the predicted targets are inhibited by cleavage instead of translation blockage. Five miRNAs and 8 putative targets were selected for qRT-PCR validation (Tables S7 and S8), but there were no significant differences among root segments using one-way ANOVA.
Gene ontology term enrichment for the complete set of putative targets identified in A. thaliana recovered several enriched functional categories and processes, including biosynthetic process (14.5%), cellular nitrogen compound metabolic process (13.7%), anatomical structure development (8.8%), response to stress (7.6%), and transport (5.8%) ( Figure 5 and Table S9). Cell wall organization or biogenesis and carbohydrate metabolic process were the 15th and 12th most represented categories, with 2.3% and 3.5% of assigned targets respectively.
All of the hydrolysis modules (cell expansion, cell separation, hemicellulose, and cellulose hydrolyzes) directly related to the cell wall (Grandis et al., 2014) are predicted to be targeted by miRNAs F I G U R E 3 Distribution of the number of shared and exclusive expressed miRNAs during aerenchyma formation in sugarcane root segments (S1 to S4) F I G U R E 4 Hierarchical heatmap clustering miRNA expression level along with the sugarcane root segments during aerenchyma formation (S1 to S4). Blue boxes enclose the four major expression clusters (Table 2). Cell separation, which comprises pectin hydrolysis, presented the highest number and diversity of predicted targets (50%), separating S1 from the other root segments and PC2 (35.2% of variation) splitting root segments consistently with the order of aerenchyma formation, that is, from S1 to S4. The miRNA species that most contributed to the separation of S1 samples in PC1 were miR6222-5p, miR395f, and miR2118-5p, whereas miR396, miR166b and miR166d, miR393b, miR397-5p, miR167b, miR164b, and some miR171 (e, h, and j) were important for segregating samples from the other segments ( Figure 6b and Table S10). In PC2, miR172e and miR399j influenced the separation of S1 and S2, while miR528, miR6223-5p, and most of the members of the miR395 family were important for segregating S3 and S4 (Figure 6b and Table S10). Although it is difficult to find clear patterns of miRNA expression among the identified miRNA members/families and segments (Figure 4), the PCA confirms that their global expression differs along the aerenchyma formation.
Indeed, S1 is the meristem, S2 is the beginning of the transformation of the parenchymatic cells into aerenchyma, and S3 and S4 are tissues already displaying aerenchyma regions and therefore the action of the cell wall degrading enzymes.  (Table 1).
The knowledge about ethylene signaling makes aerenchyma formation an interesting subject for miRNA transcriptomic studies.
Although ethylene has been pointed out as the major hormonal trigger of aerenchyma formation, its levels decrease as the aerenchyma develops (Tavares et al., 2018). Key players of the signal cascade have been described (He, Morgan, & Drew, 1996) and one scERF1 named scRAV1 acting as negative regulator upon pectin hydrolysis at the onset of aerenchyma formation was identified . miR156 and miR172 are known negative regulators of AP2/ ERF transcription factors implicated in flower organogenesis (Yant et al., 2010), developmental timing (Wu et al., 2009), nodulation (Nova-Franco et al., 2015;Wang et al., 2014), ripening (Bi, Meng, Ma, & Yi, 2015), and apical dominance (Schwab et al., 2005). Thirteen miRNAs were predicted to target transcripts related to ethylene production, perception, signal transduction, and ethylene-mediated transcriptional activity (Table 3). Among those, miR172, which was only detected when the aerenchyma is already visible (S3 and S4), targets transcripts from the last two of the mentioned categories.
Here, we observed that cell wall organization or biogenesis and carbohydrate metabolic process were well-represented ontology categories among putative miRNA target genes ( Figure 5 and Table S9).
Those include transcripts related to cell expansion and separation, and cellulose and hemicellulose hydrolysis modules. Thirty-one out of 39 miRNAs are predicted to target at least one cell wall-modification-related transcript (Table 2). miRNA6222-3p (Tables S5 and S6) stands out for targeting transcripts coding for 5 different enzyme activities.
Except for cellulose hydrolysis, all other cell wall degradation modules are predicted to be targeted by this miRNA. In addition, all degradation modules are targeted by members from miR171 and miR395 families (Table 2).
Hemicellulose hydrolysis-related targets are around 22% of cell wall-related transcripts targeted by miRNA during aerenchyma formation. Arabinoxylan debranching and β-glucan hydrolysis occur during this process Leite et al., 2017).
miR6223-5p showed decreasing levels along the root segments, which suggests that the targeted feruloyl esterase mRNA would be de-repressed as the aerenchyma develops. This esterase breaks ferulic acid bridges between arabinose residues that decorate arabinoxylan backbones. Besides being responsible for the cross-linking among arabinoxylan chains, ferulic bridges may also anchor lignin to arabinoxylan in grasses (Buanafina, 2009). As a result, when ferulic acid is absent, arabinosyl residues can be further hydrolyzed, allowing the enzyme attack on xylan. Consistent with this, the levels of arabinose decrease from S1 to S4, whereas the lignin content remains stable along the root (Leite et al., 2017). Thus, the control of expression of feruloyl esterases by miRNAs might be important because it is likely to promote loosening of interchain connections.
Such modification would allow the action of hydrolases, producing the composite that confers impermeability to gasses in the aerenchyma (Leite et al., 2017).
Cell separation encompasses half of the targeted cell wall-related transcripts (Table 2) and mostly relies on pectin degradation Leite et al., 2017;Tavares et al., 2019). This module includes carbohydrate-active enzymes such as glycosyl hydrolases, esterases, and lyases (Grandis et al., 2014. The most abundant pectin in plant cell walls is homogalacturonan and its de-acetylation and de-methylation by pectin acetyl and methylesterases, respectively. This occurs prior to the attack on α-1,4 linkages between galacturonic acids by endopolygalacturonases or by eliminative cleavage by pectate lyase (Kohli & Gupta, 2015;Latarullo et al., 2016;Yadav, Yadav, Yadav, & Yadav, 2009).
β-galactosidase is also thought to increase cell wall porosity through the hydrolysis of galactan side chains from heavily decorated pectin Leite et al., 2017). Transcripts coding for all 5 enzyme activities are among putative targets for miRNA expressed during aerenchyma formation ( Table 2).
Also, overexpression of miR156 in alfalfa led to disturbed pectin content (Aung et al., 2015). These findings, together with the observation by Leite et al. (2017) that galactose consistently decreases and β-galactosidase increases  during aerenchyma formation are consistent with the hypothesis that miR156 might be a negative regulator of pectin degradation through β-galactosidase within sugarcane roots ongoing aerenchyma formation.
The cell separation module comprises not only pectinases but also fasciclin-like arabinogalactan proteins. They are a class of hydroxyproline-rich proteins highly abundant within cell walls and plasma membrane (Harpaz-Saad et al., 2011;Shi, Kim, Guo, Stevenson, & Zhu, 2003). Ethylene precursor 1-aminocyclopropane-1-carboxylic acid is a putative key mediator of the role played by fasciclin-like proteins upon cell wall regulation within roots (Xu, Rahman, Baskin, & Kieber, 2008). These proteins have also been implicated in programmed cell death (Gao & Showalter, 1999), which is a feature of aerenchyma formation in sugarcane roots (Leite et al., 2017). Fasciclin-like arabinogalactan-related mRNAs alone were predicted to be targeted by 12 miRNAs. It has been suggested that the role played by β-galactosidases on pectin porosity could be mediated by the interaction with arabinogalactan-like proteins due to the removal of galactosyl side chains (Dean et al., 2007;Lamport, Kieliszewski, & Showalter, 2006;Seifert & Roberts, 2007). The transcriptional pattern of miR171i, miR171k, and miR166d displayed a tendency to increase along the sugarcane root segments (Figure 4). Endopolygalacturonase is also a key enzyme for aerenchyma development in sugarcane roots .
Endopolygalacturonase mRNA levels increase upon transcriptional de-repression by scRAV1 at the same time as ethylene declines, and major changes regarding pectins become visible (Leite et al., 2017;Tavares et al., 2018). Curiously, the blockage of ethylene perception using the drug 1-methylcyclopropene (1-MCP) delays but is unable to prevent aerenchyma development (Tavares et al., 2018).
miR6223-5p, miR528, and miR2118 decreasing levels from S1 to S4 might indicate that translational activity upon target transcripts (pectin acetyl esterase and endopolygalacturonase) is de-repressed ( Figure 4). These results suggest that, besides the hormone and transcriptional repression played by ethylene through scRAV1, endopolygalacturonase hydrolysis could also be subjected to miRNA regulation. Furthermore, the connection between ethylene perception and action, and pectin degradation through scRAV repression might also be influenced by miRNA regulation. The present work showed that miR156 is predicted to target scRAV1  transcription factor. This result is consistent with the observation that RAV1 expression was affected under miR156 overexpression in rice (Xie et al., 2012). Moreover, alteration in this miRNA led to a reduction in lignin content, higher cellulose accumulation, increased sugar yield, and improved digestibility in P. virgatum, maize, alfalfa, and poplar (Chuck et al., 2011;Fu et al., 2012;Rubinelli, Chuck, Li, & Meilan, 2013 The global analysis of miRNA expression levels revealed a clear distinction among all root segments, and S1 was completely separated from S2, S3, and S4 in PC1 ( Figure 6). Several miRNAs influenced this separation (Table S10), but, probably, they are more related to the different developmental root zones rather than with the formation of aerenchyma itself. PC2 discriminates S2 completely from the other segments, suggesting that this segment is where most of the epigenetic control via miRNA of aerenchyma formation takes place in sugarcane roots. S3 and S4 were closely grouped, with miR172e and miR399j being relevant in this separation. In these segments, the ethylene signaling, programmed cell death, and pectin degradation modules are already advanced, and the formation of air spaces is visible (Leite et al., 2017). The miRNA expression levels alone (Figure 4) make it difficult to define clear patterns able to explain the aerenchyma formation. However, the global analysis reinforces the fact that the segments differ among each other and miRNA expression possibly contributes to this process.

| CON CLUS IONS
Our results show that miRNA epigenetic control in sugarcane is predicted to target several genes within a range of GO categories, including biosynthetic process, cellular N metabolic process, anatomical structure development, to name but a few. Our analyses, focused on the cell wall-and hormone-related processes in the roots of sugarcane, suggest the existence of a layer of epigenetic control via miRNA. Considering the modularity of aerenchyma formation (ethylene action, cell expansion, cell separation, hemicellulose, and cellulose hydrolysis), the sugarcane miRNAs found in this work seem to be involved in hormone perception, cell separation, and hemicellulose hydrolysis.
Since sugarcane is an important bioenergy crop, our findings may have some impact on the production of bioenergy in the future, as the aerenchyma formation process offers opportunities to gain control of the whole-plant cell wall hydrolysis and biomass saccharification.

ACK N OWLED G M ENTS
The authors thank the Tropical Medicine Institute facility from the

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest associated with the work described in this manuscript.