tor sequences, clean sequencing reads from small RNA libraries were summarized for length distribution and common/specific sequences between samples, and mapped to the *L. aurea* transcriptome sequences. Reads with the same sequence were grouped and termed 'unique reads'. The clean sequences were first used for mapping to an entire integrated transcriptome of *L. aurea.* The remaining unique sequences were performed for identifying the known/conserved miRNAs from miRBase 20.0 (<http://www.mirbase.org/index.shtml>) without mismatches and gaps. After that, the sequences matching non-coding RNAs included rRNAs, tRNAs, small nuclear RNAs (snRNAs), and small nucleolar RNAs (snoRNAs) deposited in the Rfam (<http://rfam.xfam.org/>) and NCBI GenBank (<http://www.ncbi.nlm.nih.gov/genbank/>) databases were removed. miRDeep2 software \[[@CR83]\] was then used to predict novel miRNAs from the left unannotated sRNAs and the stem--loop structures of pre-miRNAs were also constructed.

Comparison of miRNA expression profiles between CK and MJ-treated *L. aurea* seedlings {#Sec14}
--------------------------------------------------------------------------------------

The abundance of all miRNAs was normalized to the transcript expression level per million reads (RPM). The log2-transformed RPM value for miRNA expression was used to generate heat map by Heml \[[@CR84]\]. For fold change of MJ to CK, average abundance of three MJ samples and three CK samples was calculated respectively. Then the formula: fold change = log~2~ (MJ100/CK) were performed. The significance (*P*-value) of the miRNA expression difference between the CK and MJ100 libraries was calculated with the formula described previously \[[@CR85]\]. The false discovery rate (FDR) was also estimated, to determine the threshold of *P*-value. Fold change (\|log~2~ (MJ100/CK)\| \>1) and FDR (\<0.01) were combined to identify differentially expressed sequences, which were defined as MJ responsive miRNAs. If the normalized expression was zero, it was adjusted to 0.01.
