Introduction

The root apical meristem (RAM) of the primary root is set up early in embryogenesis. Further root meristematic regions are initiated throughout the life of a plant at the site of lateral root formation and, in the case of legumes such as Medicago truncatula, also at the site of nitrogen-fixing root nodule formation. Unlike Arabidopsis thaliana, M. truncatula has a basic open root meristem (Groot et al. 2004). Therefore, M.truncatula can be used as a model for studies on root stem cells in species with an open meristem organisation. Post-embryonic meristems are formed as a response to spatial and/or environmental cues and share common regulatory mechanisms (Liang et al. 2007; Harris and Dickstein 2010) as well as being dependent on the relative concentrations of several plant hormones, in particular auxin and cytokinin (Su et al. 2011). The interaction of these hormones triggers the establishment of a stem cell niche, which then develops into a root meristem. This is the ultimate source of all new tissue in growing roots or nodules. As cells produced in the RAM migrate beyond its influence, they undergo differentiation into the various cell types of the mature root. This region of differentiating cells is the elongation zone (EZ) and lacks meristematic tissue. It can therefore be expected that genetic factors exhibiting differential expression between EZ and the RAM-containing root tip (RT) may be associated with regulation of the meristematic environment. Transcriptomic profiling between root meristematic and non-meristematic tissue has been used to compare specific gene expression in a variety of species (Kyndt et al. 2012; Haerizadeh et al. 2011).

In addition to studies comparing tissue types, in vitro explant tissue culture systems have enabled the direct comparison of many of the molecular factors involved in stem cell propagation and differentiation between tissues exposed to various exogenous treatments (Nolan et al. 2003; Holmes et al. 2008; Imin et al. 2008). In the M. truncatula wild type cultivar Jemalong, auxin addition to leaf explant tissue induces the formation of callus consisting of undifferentiated cells enriched in root stem cells that subsequently form roots. Conversely, auxin- and cytokinin-treated explants also form callus tissue but do not develop roots (Imin et al. 2007). Previous studies in M. truncatula have shown that a number of genes involved in stem cell niche formation including homologs of WUSCHEL-RELATED HOMEOBOX 5, PLETHORA1 and 2, BABY BOOM 1, LATERAL ORGAN BOUNDARIES-DOMAIN 29 and LATERAL ROOT PRIMORDIA 1 are sharply up-regulated in auxin explant callus tissue in comparison to those treated with auxin and cytokinin (Imin et al. 2007; Holmes et al. 2010). In this manuscript, we make comparisons between the molecular factors active in root forming callus (RFC) and non-root forming callus (NRFC) and between the root meristem (RT) and (EZ) tissues from seedling roots. This strategy represents a robust and simple system for examining the differential expression patterns of genes involved in establishing and maintaining root meristems, particularly those that may respond differentially to auxin or auxin plus cytokinin.

Whole transcriptome studies have allowed the direct comparison of gene expression levels between tissues or between variations in growth conditions and have identified factors involved in processes such as stress response, nutrient uptake and root initiation and growth (Imin et al. 2007; Hao et al. 2011; Lu et al. 2012). Such studies have identified a number of key transcription factors, a major target group of miRNAs (Carrington and Ambros 2003), involved in the formation of root stem cells (Imin et al. 2007; Holmes et al. 2010). More recently, deep sequencing techniques have allowed the rapid elucidation of entire transcriptomes including the smallRNA (smRNA) component containing miRNAs (Zhao et al. 2010; Song et al. 2011). miRNAs have emerged as major regulators of organogenesis through the targeting of many developmental genes including key transcription factors (Sun 2012). They are derived from miR genes which produce single stranded RNA transcribed by RNA polymerase II (Xie et al. 2005). miR gene primary transcripts (pri-miRNA) contain an imperfect hairpin structure, which is recognised by a ribonuclease III-like enzyme of the DICER-LIKE (DCL) family, usually DCL1 (Fahlgren et al. 2007). DCL excises a double stranded 21–24 nt section of the hairpin structure to produce a RNA duplex (the pre-miRNA). The pre-miRNA associates with RNA-INDUCED SILENCING COMPLEX (RISC), which contains an endonuclease subunit belonging to the, ARGONAUTE (AGO) family, commonly AGO1. The strand which is loaded into RISC is dependent on which AGO is present and the respective 5′ base identities of each sequence comprising the pre-miRNA (Mi et al. 2008). The strand which is loaded into AGO (the guide strand) becomes the active miRNA while the opposite strand, known as miRNA* (or passenger strand) is degraded (Eamens et al. 2009). The miRNA guides the RISC complex, via perfect or near perfect base pairing to a target mRNA which is then prevented from being translated, either through cleavage by AGO1 endonucleic activity or translational inhibition (Brodersen et al. 2008). The strong conservation of miRNAs and their targets across species suggests that miRNAs perform fundamental roles in plant physiology. Many conserved miRNAs target transcription factors involved in key processes such as meristem initiation and maintenance, organogenesis, vascular and floral patterning, hormone transport and, in legumes, nodule development (Laufs et al. 2004; Mallory et al. 2005; Williams et al. 2005; Zhou et al. 2007; D’Haeseleer et al. 2011; Kaufmann et al. 2010). A limited number of studies have used deep sequencing to examine miRNA expression in M. truncatula (Szittya et al. 2008; Devers et al. 2011; Wang et al. 2011; Chen et al. 2012b) while others have used high throughput 454 sequencing (Jagadeeswaran et al. 2009; Lelandais-Briere et al. 2009), a technique which produces significantly less data. While in A. thaliana, miRNA profiling through deep sequencing has been performed on various root cell types through cell sorting (Breakfield et al. 2012), to our knowledge no direct comparison of miRNA populations within RT and EZ tissues has been done. Nor have previous studies examined the expression patterns of miRNAs induced in vitro in cells subject to exogenous auxin or auxin and cytokinin treatments.

We deep sequenced the smRNA component of RT, EZ, RFC and NRFC. Small transcriptome expression patterns between RT and EZ as well as RFC and NRFC were analysed. 21 nt sequences with high differential expression patterns were further analysed to identify pre-miRNA hairpin structures. We identified 83 previously reported miRNAs, 24 new to M. truncatula, in 44 families. Thirty-eight additional sequences mapping to genomic locations capable of forming canonical pri-miRNA were also identified, and may represent novel miRNAs. Analysis showed low differential expression in abundantly expressed miRNAs in RT and EZ. This contrasted with many highly differentially expressed miRNAs identified between RFC and NRFC. A number of miRNAs were predicted to target genes previously shown to be differentially expressed between RT and EZ or RFC and NRFC. Conserved miRNA/target relationships were found for miR397 and miR160 while other miRNAs were predicted to target homologues of genes involved in root development.

Materials and methods

Plant materials, growth and explant tissue culture

M. truncatula cv. ‘Jemalong A17’ seeds were obtained from the National Medicago Collection (South Australian Research and Development Institute, Adelaide, South Australia) and were used for all samples. Plants were grown under controlled growth cabinet conditions; 20 °C, 70 % relative humidity, 100 μE m−2 s−1 and a 16 h photoperiod. Five grams of ‘Osmocote Exact’ (Scotts Miracle-Gro Company, Marysville, OH, USA) fertilizer was applied at sowing. Explants were cultured as described in (Imin et al. 2007). Leaves for explants were collected 4–6 weeks after sowing. P4 culture medium, based on Gamborg’s B5 medium as described by (Thomas et al. 1990) contained either 10 μM 1-naphthaleneacetic acid (NAA) (Sigma-Aldrich, St Louis, MO, USA), in the case of RFC cultures, or 10 μM NAA + 4 μM 6-benzylaminopurine (BAP) (Sigma-Aldrich), in the case of NRFC cultures. Explant containing plates were incubated in the dark at 28 °C for 10–14 days. Callus tissue was stored at −80 °C. To collect RT tissue, sterilised A17 seeds were stratified at 4 °C for 24 h and germinated at 20 °C. 12–15 h post germination, seedlings were transferred to 150 mm Petri dishes containing Fahraeus Media (Fahraeus 1957). The media was “slanted” i.e., the medium was poured in Petrie plates held at a 25° angle so only the lower half was covered. In this way, only the roots and lower portion of the hypocotyl were in contact with the medium. Each dish was sealed with Nescofilm (Brando Chem. Kobe, Japan) apart from 40 to 60 mm at the top to allow gas transfer. Dishes were placed in trays at a 15°–25° angle to vertical with black cardboard spacers between to exclude light from the roots. Seedlings were grown for 3 days at 20 °C, 150 μE m−2 s−1, 80 % relative humidity and 16 h photoperiod. The RT (distal 3 mm), which included the root cap, the root meristem and at least part of the transition zone and the EZ (distal 3-13 mm) were removed and stored at −80 °C.

Small RNA extraction

Whole RNA was extracted using a modified “Trizol method” (Sunkar and Zhu 2004). The smRNA isolation proceeded as follows: to 12 μg total RNA (50 μl), 6 μl 50 % PEG (MW 8000) and 7.5 μl 4 M NaCl was added and samples incubated at 4 °C for 30 min. Samples were then centrifuged at 16,100g at 4 °C for 15 min. To the supernatant, 3 × vol 100 % EtOH was added. Samples were precipitated overnight at 4 °C. After spinning at 4 °C at 16,100g for 30 min, the supernatant was removed and pellets were gently washed in 75–80 % EtOH for 5 min. Pellets were air dried and resuspended in 60 μl UltraPure distilled water (Invitrogen, Carlsbad, CA, USA).

Preparation of smRNA and deep sequencing

Small RNA (smRNA) samples were prepared using a modified version of ‘Preparing Samples for Small RNA Sequencing Using Alternative v1.5 Protocol’ (Part # 15002615 Rev. A, February 2009, Illumina, San Diego, CA, USA). Briefly, 200–300 ng low molecular weight RNA was used. 3′ and 5′ adaptors were ligated in separate reactions. The 3′ adaptor ligation was carried out for 2.5 h at 22 °C and the 5′ adaptor ligation was carried out for 1 h at 22 °C and then overnight at 4 °C. Constructs were isolated at each stage using 15 % TBE-Urea Gels (Invitrogen). Reverse transcription was performed using SuperScript® III Reverse Transcriptase (Invitrogen). Deep sequencing was performed with one sample per lane on a Genome Analyser IIx (Illumina). Raw data files are available for viewing at the Gene Expression Omnibus (GEO). Batch accession GSE45726 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE45726).

Data processing

Adaptor trimming of raw data and alignment to the miRBase (Griffiths-Jones 2004) and M. truncatula genome (Mtr 3.5.1) databases was performed using Novoalign v2.07.18 (Hercus and Albertyn 2011). In the case that a sequence did not align alternative databases were used. These were either Genbank (for miRN374, miRN404, miRN404a, miRN477, miRN521a, miRN521 and miRN477a) or MedtrA17_3.5.assemblies (http://jcvi.org/cgi-bin/medicago/download.cgi) (for miR172). Data normalisation was performed using EdgeR as described in (Robinson and Oshlack 2010). Input was combined with 20–25 nt data sets after removing all sequences that did not contain at least five reads in any one tissue. The trimmed mean values used to calculate the scaling factors were; Mg 30 %, Ag 5 %. Scaling factors used for normalisation were; NRFC: 1.0330943, RFC: 1.0744541, RT: 0.8975524, EZ: 1.0037195. The majority of M. truncatula miRNAs listed in miRBase are aligned to the Mtr 3.0 database. For this reason, when identifying new locations for previously published M. truncatula miRNAs, Mtr 3.5.1 was substituted with Mt 3.0. Sorting by read length and collation by read count analysis was performed using in-house Perl scripts. Duplex energy calculations were done using the RNAduplex module of ‘Vienna RNA package’ (Lorenz et al. 2011) using the default parameter settings. Duplex energies were calculated as kilocalories per mol. Individual miRNA hairpin structures were produced using Mfold-quikfold (RNA 3.0, linear) (Zuker 2003). When multiple hairpin diagrams were generated, a modified version of UNAFold.pl was used (Markham and Zuker 2005). * strand prediction was based on quikfold structures and allowed for a miRNA 2 nt 3′ overhang. Target prediction was performed using the online version of TAPIR (Bonnet et al. 2010). The minimum score threshold was set to nine and the free energy ratio 0.7. Database for target prediction was the Affymetrix GeneChip Medicago Genome Array MTGI_V11 (Affymetrix, Santa Clara, CA, USA). Degradome mapping was performed using Cleaveland v3 (Addo-Quaye et al. 2009). Degradome datasets were root (GSE26218) and flower (GSM769293) available from the NCBI gene expression omnibus (http://www.ncbi.nlm.nih.gov/geo/).

Mature miRNA validation by quantitative RT-PCR

Mature miRNAs were validated via qRT-PCR using Taqman® MicroRNA Assays [Applied Biosystem (ABI)] with 100 % probe sequence complementarity. A Taqman Assay® probe specific to M. truncatula U6 was used as a housekeeping control. qRT-PCR was performed for 40 cycles on an ABI-7900. SuperScript® III Reverse Transcriptase (Invitrogen) was used for cDNA synthesis. Assays were multiplexed using both miRNA-specific and random primers.

Validation of miRNA-guided cleavage of target mRNAs by 5′ RACE

Target validation employed a modified version of 5′ RACE (Llave et al. 2002) using the GeneRacer Kit (Invitrogen). Prior to sequencing, cleavage products were detected via PCR using the GeneRacer 5′ forward primer and a gene-specific reverse primer; APETALA2 (Mtr.41294.1.S1_at) GGC-TAT-GGT–GGT-GCC-TGT-GAG-GAC-TT, HB1 (TC175771) TTC-AAA-ACT-GCT-TAT-GAG-ATG-TCC-GTA-ACT-TTC-TAC. SuperScript® III Reverse Transcriptase was used for cDNA synthesis.

Results

Analysis of total small RNA population

Four smRNA transcriptome libraries for M. truncatula were constructed for deep sequencing. Two were isolated from 10-day-old RFC or NRFC explants when root-forming stem cell niches were being generated in the RFC only. A third library was constructed from the tissue containing the meristematic region (distal 3mm of the root) and a fourth from the EZ (distal 3-13 mm of the root). Both were collected from 3-day-old seedlings. Each library produced between 24.5 million and 29.2 million reads (Table 1). Adaptor trimming and sorting by length revealed, for each tissue set, 24 and 21 nt sequences as the most common (Fig. 1). This distribution, with a predominant clustering from 21 to 24 nt is consistent with other smRNA deep sequencing studies across a range of species (Li et al. 2011; Peng et al. 2011; Yu et al. 2011; Chen et al. 2012a). At the time of writing, miRBase contained 7,703 plant sequences, of which 1,039 (58.6 %) were 21 nt. In order to maximise the likelihood of identifying novel miRNAs, we restricted further analysis to this set. The total number of individual sequences for any given sample ranged from 318 to 440 K (Table 1). Exact alignment to the M. truncatula 3.5.1 genome occurred for 62.1 % (RT), 63.5 % (EZ), 62.7 % (RFC) and 60.8 % (NRFC) of sequences. In order to analyse sequence frequency, identical sequences were collated (Table 2). A comparison of the four samples showed that between 84 and 87 % of sequences, occurred only once, whereas only ~1.0 % or less had ≥50 counts.

Table 1 Read counts per tissue
Full size table
Fig. 1
figure 1

Sequence length distribution. Total read counts of each tissue sample distributed, after adaptor trimming, by sequence length. The distribution, with peaks at 24 and 21 nt, was typical of smRNA deep sequencing profiles in plants

Full size image
Table 2 Frequency of 21nt read counts
Full size table

Identification of conserved miRNA, novel to M. truncatula

The full dataset of 21 nt sequences was aligned against miRBase (release 19). We identified 83 conserved sequences in 44 families. Fifty-nine of these had been previously reported in M. truncatula and 24, in 8 families, were novel to M. truncatula (Suppl. Table S1, Suppl. Fig. S1). Members of conserved families miR165, miR181 and miR397 were identified for the first time in M. truncatula. miRNAs with the highest read counts were from highly conserved families. Isoforms of miR159, miR172 had extremely high counts (≥100 K reads) and isoforms of miR166, miR319 and miR396 (two forms) were also highly abundant (≥10 K reads). Conserved isoforms of miR172, miR166, miR319 and miR396 were identified that have not previously been reported in M. truncatula. Counter to previous genome wide small RNA studies (Chen et al. 2012b; Zhou et al. 2012), we found members of the miR156/157 family under-represented while counts matching isoforms of miR395 were negligible. A novel isoform of mtr-miR166 was the most prevalent within its family and had the third highest count overall. Interestingly, this sequence is located on the opposite arm to mtr-miR166g within the reported precursor and, in other species, has extensively been reported as the passenger strand (Zhang et al. 2009; Song et al. 2009; Ma et al. 2010). Two other passenger strand sequences; miR1507* and miR1509*  also had counts markedly higher than sequences previously reported as guide strands (Suppl. Table S2). Additionally, we found 11 previously reported M. truncatula miRNA sequences at novel genomic locations (Suppl. Table S2, Suppl. Fig. S2).

Identification of non-conserved miRNAs, novel to M. truncatula

In order to identify novel miRNAs, 1,197 non-conserved sequences from our data showing high expression levels were examined. The genomic region of each was searched for nearby, highly complementary sequences (duplex sequences). These sequences were analogous to, but not necessarily the same as, the miRNA* strand and were treated as indicators of the presence of a characteristic miRNA hairpin structure. In general, duplexes producing ΔE ≤ −28 kcal/mol were associated with hairpin structures conforming to canonical pre-miRNA (Suppl. Table S3, Suppl. Fig. S3). Of the 38 novel miRNAs predicted, only four had duplex ΔE greater than −28 kcal/mol and none greater than −24 kcal/mol. To validate this procedure, an identical analysis of all plant miRNAs in miRBase was performed (data not shown) revealing 50 % of the miRNAs had duplex energy less than ΔE −28 kcal/mol. The miRNA* is degraded in vivo during miRNA biogenesis. Nevertheless, it is often detected in deep sequencing data, albeit at much reduced levels compared to the miRNA (Mi et al. 2008; Breakfield et al. 2012). We calculated the * sequences for each predicted structure and found 18 for 35 predicted miRNAs in the deep sequencing data (Suppl. Table S3), suggesting many of the stemloop structures were targeted for miRNA processing by DCL1.

Validation of mature miRNAs by qRT-PCR

Validation of mature miRNAs was performed by qRT-PCR analysis (Fig. 2) using Taqman® MicroRNA Assays. Conserved miRNAs miR172, miR166-5p and miR166-3p were detected. Congruent with the deep sequencing data miR166-5p had stronger expression than miR166-3p. Also validated were a predicted novel miRNA, miRN304 (identifiers prefixed 'miRN' are used throughout the text to denote sequences without an exact match in miRBase and do not signify pre-existing family names). This is contained within the same precursor as mtr-miR2087, a previously validated M. truncatula miRNA (Szittya et al. 2008). Both sequences were mapped to singular genomic locations. As would be expected for miRNAs processed from the same precursor the expression pattern between tissues was very similar.

Fig. 2
figure 2

qRT-PCR validation of mature miRNA sequences detected by deep sequencing. Normalised Taqman assay data. Conserved and novel 21 nt sequences were detected in the small RNA libraries of each sample. A Probe specific to U6 mRNA (not shown) was included as a mean centred control for normalisation. Error bars indicate 1 standard error of three biological replicates

Full size image

Validation of miRNAs-guided cleavage of target mRNAs by 5′ RACE

To further validate our data, we detected two predicted cleavage products using 5′ RACE. Our data contained four isoforms of miR172. A conserved sequence (ath-miR172d) represented 96 % of miR172 isoforms and was amongst the most prolific in our data but had not previously been reported in M. truncatula. Using TAPIR, each isoform was predicted to target the APETELA2 (AP2) domain transcription factor Mtr.41294.1.S1_at. In A. thaliana, miR172 targets multiple AP2 domain genes regulating flowering time through both cleavage and translational attenuation (Aukerman and Sakai 2003). We detected cleavage at an identical 21 nt sequence to the target region of A. thaliana APETALA2 (At4g36920) (Fig. 3a). miRN304 is contained on the same precursor as miR2087 and its sequence is unique to M. truncatula. We predicted miRN304 to target HB1 (TC175771), a LOB-like HD-Zip I transcription factor involved in lateral root formation (Ariel et al. 2010). Cleavage of HB1 by miRN304 was confirmed (Fig. 3b).

Fig. 3
figure 3

5′ RACE confirmation of predicted miRNA target sites. a confirmation of cleavage of Mt- APETALA2 by miR172. Targeting was predicted for four isoforms of miR172 (miR172d shown) b confirmation of MtHB1 cleavage by novel miRNA miRN304. In both cases cleavage was between positions 10–11 (arrows) relative to miRNA

Full size image

Differential expression analysis

Differential expression for RT/EZ (Table 3) and RFC/NRFC (Table 4) was calculated for all conserved and predicted novel miRNAs (including * strands) listed in Suppl. Tables S1 and S2. Forty-one had ≥1.5-fold change between RT and EZ. Only miRN482 and miRN1112 had differential expression greater than threefold. The majority (68 %) of miRNAs with ≥1.5-fold change were up regulated in the RT. Differential expression of miRNAs between RFC and NRFC was much higher. Thirty-seven miRNAs were expressed greater the fivefold. Of these, 15 had differential expression fold changes >50.

Table 3 Differential expression between root tip (RT) compared to elongation zone (EZ)
Full size table
Table 4 Differential expression between root forming callus (RFC) compared to non-root forming callus (NRFC)
Full size table

Prediction of miRNA targets: differentially expressed genes

Genes previously identified as differentially expressed between RT and EZ (Holmes et al. 2008) plus RFC and NRFC (Holmes et al. 2010) were used as a database for target prediction. Targets were predicted for all conserved and novel miRNAs identified listed in Suppl. Tables S1 and S3. Thirty-six miRNA/target interactions were identified. Fourteen miRNAs were from conserved families and five were novel. Only four miRNAs were predicted to target genes differentially expressed between RFC and NRFC. Twelve potential miRNA/target pairs showed an inverse miRNA/target relationship (an increase in miRNA expression correlates with a down regulation of the target). Three highly similar isoforms of miR396 had multiple targets including transcription factor homologues (Mtr.31411.1.S1_at and Mtr.27107.1.S1_at). Other predicted targets have potential roles in root formation. These include: ENDO-1 (Nicol et al. 1998), Mtr.13919.1.S1_s_at; ARF1 (Inukai et al. 2005), (Mtr.39233.1.S1_at) and; a LACCASE-LIKE protein (Mtr.41285.1.S1_at), (Zhang et al. 2012).

Prediction of miRNA targets: whole transcriptome

We predicted targets for conserved miRNAs not previously identified in M. truncatula (Suppl. Table S4) and predicted novel M. truncatula miRNAs (Suppl. Table S5). Target predictions were disregarded if they had insignificant expression levels or indistinct expression patterns or if they were considered not be involved in the regulation of root development or growth. We identified 32 targets for 16 conserved miRNAs (Suppl. Table S4) and 28 targets for 15 novel miRNAs (Suppl. Table S5). Transcription factors involved in stress responses and development featured amongst the predicted targets of these miRNAs. Additionally, the sequences from both tables were also mapped to the available M. truncatula degradome databases contained at GEO datasets. Of the novel sequences, we were able to match only miRN404 to degradome data. The current M. truncatula degradome resources are limited to one flower (Zhai et al. 2011) and one root (Devers et al. 2011) database. As more become available, we feel it is likely that data matching more of these predictions will emerge. This is illustrated by the fact that no match was identified for miRN304 despite confirmation of HB1 cleavage by 5′ RACE. We also confirmed cleavage of mtr-APETALA (Mtr.49221.1.S1_at) but were unable to locate a coinciding tag within the root degradome database. The non-conserved nature, however, of many these miRNAs suggests that until further validation is undertaken, these miRNA/target predictions should be treated with caution.

Discussion

The development of global gene expression tools such as deep sequencing has enabled the study of differential expression within entire transcriptomes. In M. truncatula transcriptional profiling has previously been performed using Affymetrix GeneChip® technology to identify genes involved in root formation. Holmes et al. (2008, 2010) compared root meristematic (RT) and non-meristematic (EZ) tissue as well as in vitro RFC and NRFC tissue identifying a number of differentially expressed genes involved in the regulation of root meristem and root stem cell niche formation. miRNAs perform a diverse range of functions including regulation of development and cell fate (Carlsbecker et al. 2010; Boualem et al. 2008; Guo et al. 2005) by targeting key developmental genes including transcription factors (Sun 2012). We looked for known or putative miRNAs that potentially target genes involved in root development. Of particular interest were the AP2-domain PLETHORA1 and 2, BABYBOOM1, GRAS transcription factors SHORTROOT and SCARECROW, homeodomain containing WOX5 and members of the AUXIN RESPONSE FACTORS (ARFs) family, all of which are conserved and perform crucial developmental functions. As miRNAs have been previously shown to be involved in plant development we believed it likely that miRNA profiling in the same tissues would show a similar differential expression pattern.

Root tip comprises multiple tissue types. Primarily, it consists of the RAM which contains the functional initials: pluripotent stem cells also known as the quiescent centre (QC). These serve to maintain the surrounding cells, the structural initials, also stem cells, in a pluripotent state and to slow their division. The RT also contains some of the daughter cells of the structural initials, which have moved beyond the influence of the QC. Additionally, RT tissue contains cells of the root cap. In contrast, the EZ consists of elongating daughter cells of the RAM undergoing differentiation (Tsukagoshi et al. 2010). When leaf explants are cultured on media containing auxin, cells at the wounded leaf margins are induced to function as pluripotent stem cells able to form callus and ultimately roots (Rose et al. 2006). If the media also contains cytokinin, a repressor of root growth, calli formation occurs but root formation does not.

The growth of in vitro roots requires the formation of stem-cell niches and ultimately meristems at the site of initiation and it is likely that auxin plays a key role by altering the redox environment to that required by the QC (Jiang et al. 2003). The RFC tissue was collected at the early stages of auxin-induced calli formation, a time point at which they are very likely enriched in cells that have committed to root formation, especially when compared to NRFC tissue.

Previous comparative studies showed 608 genes differentially expressed (over twofold) between RT and EZ and 274 between RFC and NRFC (Holmes et al. 2008, 2010). We deep sequenced the small RNA component of these tissues to identify miRNAs that target genes previously identified as involved in root development. After initial data analysis we retained only 21 nt sequences with at least 50 reads in any one tissue sample (Tables 3, 4). We found the greatest expression change was between RFC and NRFC. Between RT and EZ, no miRNA had differential expression more than fivefold compared to 37 miRNAs between RFC and NRFC. We found 17 sequences with ≥50-fold differential expression between RFC and NRFC. Of the 24 miRNAs expressed with ≥10-fold change between RTC and NRFC, 12 were higher in RFC and another 12 in NRFC. Of those up-regulated in NRFC, nine are conserved all of which, apart from miR160, have been shown to be up-regulated in response to drought or salinity (Sun 2012; Trindade et al. 2010). Overall, these nine miRNAs also had much higher read counts than those up-regulated in RFC suggesting strong induction. It is possible therefore that many conserved miRNAs are stress responsive and that the repression of root growth by cytokinin shares genetic elements with abiotic stress responses. Some evidence for this has been seen previously in studies suggesting that increased cytokinin levels may instigate drought stress responses in A. thaliana (Nishiyama et al. 2011) and rice (Peleg et al. 2011).

We compared predicted targets for all conserved and novel miRNAs with the data produced by the studies of Holmes et al. (2008, 2010). Initially, we looked for miRNA/target interactions with inverse expression patterns as this is the most common, but not exclusive, mode of mRNA action. When we applied the criterion of at least a two-fold expression change in the target transcript between tissues, we identified 12 potential miRNA/target pairs with inverse expression patterns (Table 5). Notable amongst these is miR397 and its predicted target Mtr.9478.1.S1_at, a Cu+-binding/laccase homologue. Although miR397 is highly conserved and has been widely studied, no family member has previously been reported in M. truncatula. A BLAST search also showed very strong conservation of the predicted target region amongst laccase-like genes across multiple species. Comparison of the M. truncatula predicted target site and those validated in A. thaliana (Abdel-Ghany and Pilon 2008) show very high conservation. We note with interest that the miR397, miR408 and miR398 family members represented in our data share highly similar expression patterns. Notably, each was strongly up-regulated in the cytokinin treated NRFC tissue. All have been shown to be involved in maintenance of copper homeostasis (Abdel-Ghany and Pilon 2008; Yamasaki et al. 2009).

Table 5 miRNAs predicted to target differentially expressed genes
Full size table

Other miRNAs and their targets did not show the same inverse expression pattern but may also be part of a functional miRNA/target interaction. miR160 is conserved and targets ARFs in multiple species (Sun 2012). Eight ARFs were predicted to be targeted by mtr-miR160 (isotypes a, b, d, e). Three were ARF10 homologues (Mtr.39233.1.S1_at, Mtr.25836.1.S1_at, Mtr.10054.1.S1_at) and five ARF17 homologues (Mtr.33918.1.S1_at, Mtr.33918.1.S1_at, Medtr5g061220, Mtr.21419.1.S1_at, Medtr5g060630). Mtr.39233.1.S1_at, a homologue of A. thaliana ARF10, is up-regulated in meristematic tissue (Holmes et al. 2008) and controls stem cell niche size in the root cap. In the root tip, miR160 is up-regulated in conjunction with ARF10 presumably to restrict the size of the niche (Wang et al. 2005). In A. thaliana, miR160 has also been shown to be involved in adventitious root development where it is induced by ath-ARF6 and targets ath-ARF17 in a self regulating feed-back loop (Gutierrez et al. 2009).

miRNAs were predicted to target important developmental genes

Alignment of all 21 nt read sequences against the miRBase database revealed 24 sequences not previously reported in M. truncatula (Suppl. Table S1, Suppl. Fig. S2) but which are conserved in other plant species. We used duplex energy calculation and RNA precursor folding to predict 16 conserved miRNAs, novel to M. truncatula in 8 families (Suppl. Table S3 and Suppl. Fig. S2). In these cases, the duplex energy threshold was relaxed to −23 kcal/mol as sequence conservation is a strong miRNA predictor, therefore less folding stringency was required than for novel miRNA prediction. Target prediction showed the likelihood that the miRNA/target relationship conserved across plant species is largely maintained in the M. truncatula miRNAs predicted here. We predicted targets which we considered to be biologically significant for every conserved miRNA listed in Suppl. Table S3 except for miR166i. This result reinforces conservation as a strong predictor of miRNAs. We identified three members of the miR171 family; mtr-miR171b, mtr-miR171f and a novel sequence (provisionally designated miR171i) conserved across a high number of plant species. Each of these sequences target homologues (Mtr.37696.1.S1_at and Mtr.41007.1.S1_at) of SCARECROW-LIKE, a GRAS family transcription factor. The region of greatest conservation is small (~30 nt) but contains the target region for each of the miR171 isoforms. In A. thaliana, SCARECROW is expressed in the cortical/endodermal initials, endodermis and QC (Hashimura and Ueguchi 2011) where it contributes to stem cell maintenance (Sabatini et al. 2003). In A. thaliana, miR171 has been shown to target SCARECROW-LIKE-II (At2g45160) and SCARECROW-LIKE-III (At3g60630) (Llave et al. 2002).

Amongst the predicted targets for all novel (Suppl. Table S5) and conserved (Suppl. Table S4) miRNAs, we identified 11 targets that encode transcription factors involved in various developmental processes. We also predicted miRNA targeting of nine F-box domain annotated proteins. These are known to be involved in the ubiquitination pathway and are often activated during stress response (Sepulveda-Garcia and Rocha-Sosa 2012). Included were miRN395, which shares a 19/21 nucleotide identity with miR1510 and is conserved in Glycine max and G. soja. miRN395 is predicted to target (Mtr.46578.1.S1_at), a homologue of A. thaliana F-box protein SLOWMO, required for auxin homeostasis. SLOWMO mutants in A. thaliana show increased lateral root density (Lohmann et al. 2010). miRN246 was predicted to target an EXPANSIN homologue (Mtr.22752.1.S1_s_at), a class of protein involved in cell elongation and asymmetric cell division, a process seen during functional initial division in root meristems. Multiple EXPANSIN proteins have been shown to be active in roots of developing seedlings (Shin et al. 2005).

miR156, miR159 and miR172 were strongly expressed in all tissues

The majority of well-documented conserved miRNAs identified have been reported to have multiple roles covering both development and stress response. While the function of many of these have been well characterised we note a number had very strong expression in our samples. We found strong expression of mtr-miR156 and miR164a, b, c in RT compared to the other tissue samples. Members of the miR156 and miR172 families operate in the same pathway and are conserved across a broad range of plant species. miR156 targets the SBP/SPL family of transcription factors (Xing et al. 2010) and miR172 targets AP2 genes (Aukerman and Sakai 2003). We confirmed cleavage of mtr-APETALA2 (Mtr.41294.1.S1) through 5′ RACE (Fig. 3) and degradome analysis (Suppl. Table S3). mtr-miR159 had the highest reads counts of all 21 nt sequences. In A. thaliana, miR159 is known to target MYB transcription factors during seed and anther development (Achard et al. 2004) and seed germination (Reyes and Chua 2007), although it has also been shown to have strong expression in RTs (Allen et al. 2007).