Global Analysis of the Genetic Variations in miRNA-Targeted Sites and Their Correlations With Agronomic Traits in Rapeseed

MicroRNAs (miRNAs) and their target genes play vital roles in crops. However, the genetic variations in miRNA-targeted sites that affect miRNA cleavage efficiency and their correlations with agronomic traits in crops remain unexplored. On the basis of a genome-wide DNA re-sequencing of 210 elite rapeseed (Brassica napus) accessions, we identified the single nucleotide polymorphisms (SNPs) and insertions/deletions (INDELs) in miRNA-targeted sites complementary to miRNAs. Variant calling revealed 7.14 million SNPs and 2.89 million INDELs throughout the genomes of 210 rapeseed accessions. Furthermore, we detected 330 SNPs and 79 INDELs in 357 miRNA target sites, of which 33.50% were rare variants. We also analyzed the correlation between the genetic variations in miRNA target sites and 12 rapeseed agronomic traits. Eleven SNPs in miRNA target sites were significantly correlated with phenotypes in three consecutive years. More specifically, three correlated SNPs within the miRNA-binding regions of BnSPL9-3, BnSPL13-2, and BnCUC1-2 were in the loci associated with the branch angle, seed weight, and silique number, respectively; expression profiling suggested that the variation at these 3 miRNA target sites significantly affected the expression level of the corresponding target genes. Taken together, the results of this study provide researchers and breeders with a global view of the genetic variations in miRNA-targeted sites in rapeseed and reveal the potential effects of these genetic variations on elite agronomic traits.

Introduction

Traditional crop breeding requires phenotypic evaluations and its success is mainly dependent on the experience of the breeder. The recent development of molecular breeding by design has led to increased breeding efficiency (Zeng et al., 2017). On the basis of advances in functional genomic research and genetic studies on families or natural populations, a large number of single nucleotide polymorphisms (SNPs) contributing to agriculturally important traits have been identified. Re-sequencing followed by SNP genotyping is one of the key strategies used for developing molecular markers for specific agronomic traits (Morrell et al., 2012). Including SNPs that change protein sequences as well as those that affect transcriptional and post-transcriptional regulation involving cis-elements is also important. For example, SNPs in gene promoters targeted by transcription factors regulate transcription, whereas functional SNPs in exon–intron junctions, protein–RNA binding sites, and microRNA (miRNA)-targeted sites contribute to post-transcriptional regulation (Zhang W. et al., 2018). Therefore, these functional SNPs are key genetic resources for crop breeding.

MicroRNAs, which are small RNAs comprising 20–24 nucleotides, function as important post-transcriptional regulators of gene expression in plants and animals (Bartel, 2004; Dugas and Bartel, 2004; Kidner and Martienssen, 2005; Chuck et al., 2009). Plant miRNAs are generated from stem-loop precursors by the DCL1 complex. Mature miRNAs are loaded into the RNA-induced silencing complex to enable AGO1 to cleave mRNAs complementary to specific miRNAs, and their miRNA^∗ strand is subsequently degraded (Shukla et al., 2008; Chen et al., 2016). Inhibited gene expression via the translational arrest mediated by the AGO1 complex is also a common phenomenon among plants (Brodersen et al., 2008). Additionally, miRNAs negatively regulate the expression of specific genes by targeting the untranslated region (UTR) or the coding region of mRNAs. Several studies revealed that miRNAs and their target genes play important roles in the vegetative phase and the transition to reproductive development (Wang et al., 2009, 2014; Tripathi et al., 2018), signal transduction pathways (Navarro et al., 2006; Liu et al., 2009; Inui et al., 2010), plant morphogenesis (Liu et al., 2010, 2019; Mao et al., 2014; Ren et al., 2018), and biotic and abiotic stress responses (Lu et al., 2013; Li et al., 2019). For example, miR156 is conserved among species and affects shoot branching and plant architecture by regulating the expression of SPL genes (Wang et al., 2009; Ferreira e Silva et al., 2014; Liu et al., 2017). The overexpression of miR156 affects several traits, including plant height, flowering time, leaf number, root length, grain yield, and branch number (Wang et al., 2009; Wang and Zhang, 2017). Thus, miR156–SPL modules coordinately regulate diverse traits in crops.

MicroRNAs control many biological processes by silencing target genes. Moreover, they have been broadly applied for manipulating eukaryotic gene expression. The variability of the DNA sequences at miRNA target sites alters the complementation between miRNAs and their target genes, which influences the target gene expression levels and the traits relevant to these genes (Mallory et al., 2004). However, the extent of the genetic variations at miRNA target sites and the agronomic traits that are altered by the variations remains unclear.

Rapeseed (Brassica napus) is one of the most important oilseed crops in the genus Brassica. It is mainly grown in Europe, Asia, North America, and Oceania, where it is an important oilseed crop. Through natural and artificial selection, rapeseed has become a multipurpose crop (i.e., vegetable crop, forage crop, and oilseed crop) (Allender and King, 2010). It contains two divergent sub-genomes, A and C, which are shaped by a whole-genome triplication followed by extensive diploidization (Parkin et al., 2005; Mandáková and Lysak, 2008). Thus, its genome structure is complex, which has hindered genomic research and the high-throughput identification of high-quality molecular markers. Researchers and breeders have bred many elite rapeseed germplasms through the laborious selection of target traits. These germplasms have been exploited to breed new cultivars and have shown great potential for increasing crop yield and quality. However, the molecular mechanism underlying trait selection remains unclear. Studying the genetic variants generated during the formation of elite rapeseed germplasms may produce valuable information regarding the genetic changes, thereby facilitating the efficient selection of elite germplasms for breeding.

Current research on important agronomic traits includes the detection of global genetic variations and genome-wide association studies following by the fine-mapping of the causal genetic variants affecting a particular trait. However, in this study, we focused on elucidating the genetic variations in the miRNA-complementary region of target genes in 210 elite rapeseed accessions. We detected three SNPs in miRNA target sites that are, respectively, correlated with the branch angle, silique number, and seed weight. The results of this study provide an overview of the genetic variations in miRNA target sites as well as additional functional SNPs relevant for rapeseed molecular breeding.

Materials and Methods

Plant Materials

For this study, 650 rapeseed accessions were obtained from the Center for Shaanxi Hybrid Rapeseed and Shanghai Academy of Agricultural Sciences, China. On the basis of the information provided by the rapeseed accession donor and our own observations, 210 elite accessions were selected from the 650 accessions as the core collection of genetically diverse rapeseed accessions with varying phenotypic traits for the analysis of genomic variations (Supplementary Table 1). To obtain homozygous genotypes, each accession was self-fertilized for more than five generations and one representative plant of each accession was selected for DNA re-sequencing. These accessions were assigned to different germplasm categories according to their growing seasons. More specifically, the accessions were designated as winter (7), semi-winter (198), and spring (5) oilseed rape (Supplementary Table 1). The accessions originated from Asia, Europe, and North America (Supplementary Table 1). Their seeds were sown in experimental field plots at the Farm Station in Yangling, Shaanxi province, China in 2016, 2017, and 2018. The accessions were grown using a randomized complete block design with two replicates. Each accession was grown in plots with four rows, which were separated by 25 cm. The 22–25 plants in each row were separated by 8–12 cm. Between 5 and 10 representative plants in the middle of each plot were selected at harvest to analyze the plant height, silique number, seed number, seed weight, oil content, branch number, and branch angle. The flowering time, ripening time, cold resistance, Sclerotinia sclerotiorum resistance, and stem lodging resistance were analyzed in the field at a representative stage. The average trait values of the inbred line replicates in each year were used for phenotypic and correlation analyses in a single environment. The correlations among traits were evaluated using the cor.test function in R (version 3.6.1).

DNA Re-sequencing and Alignment With the Reference Genome

Developing leaves were collected from individual plants for a genomic DNA extraction using the Plant Genomic DNA Kit (Tiangen, Beijing, China). The concentration and quality of the extracted genomic DNA were determined using the NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific). The DNA samples with an OD₂₆₀/OD₂₈₀ ratio of 1.8–2.0 and a total content of more than 1.5 μg were suitable for constructing sequencing libraries. An Illumina sequencing library was constructed for each accession using the IlluminaTruSeq Library Construction Kit. Paired-end sequencing was performed using the Illumina HiSeq platform, with a read length of 150 bp at each end (Novogene, Tianjin, China). Raw reads were filtered according to the following criteria: paired-end reads with >10% “N” bases and Q ≤ 10. The clean reads for each accession were mapped to the B. napus “Darmor” reference genome (version 4.1), which was downloaded from BRAD¹ (Chalhoub et al., 2014). The PCR-duplicated reads were removed by Picard MarkDuplicates. The SNPs and INDELs were called for all 210 B. napus accessions using GATK.² The identified variants were annotated using ANNOVAR.³ Additionally, AGAT-master was used to convert the gff3 file to a gtf file⁴.

Predication of miRNA-Complementary Target Regions

The miRNA target sites were predicted using psRNATarget,⁵ with default parameters and a maximum expectation value of 3.5. The targets were all protein-coding transcripts in the B. napus transcript library in BRAD (see text footnote 1). To ensure target sites were accurately identified, the putative target sites were further evaluated using the available data generated by degradome sequencing or the 5′ RNA ligase-mediated rapid amplification of cDNA ends. Incorrectly predicted target sites were removed.

Extraction of SNPs and INDELs From miRNA Target Sites and Correlation Analyses

The physical positions of target genes were obtained from the gff3 file (see text footnote 1). The BLASTN (blastn-short) algorithm was used to determine the locations of miRNA target sites within the reference genome. Variants were extracted from the genotype data using perl packages. Genotypes were filtered using PLINK2. Of the target site regions, 330 raw SNPs were detected and then a pseudo-VCF was created. The TASSEL 5 program⁶ was used to convert the VCF format to the Hapmap format. Finally, 199 high-quality SNPs (MAF > 0.05, integrity ≥ 0.7) were obtained for the single-SNP correlation analysis, which was conducted using the general linear model in the TASSEL 5 program.

Quantitative Real-Time PCR

Total RNA was extracted from young leaves using the TRIzol reagent, after which DNase I (Takara) was used to remove DNA contaminants. The total RNA (2 μg) was reverse transcribed into cDNA using PrimeScript Reverse Transcriptase (TaKaRa) and an oligo (dT) primer. The qRT-PCR was performed using the SYBR Green Real-time PCR Master Mix (TOYOBO, Japan) and specific primer pairs (Supplementary Table 2) to quantify relative transcript levels. The samples analyzed by qRT-PCR are listed in Supplementary Table 3. The ACTIN7 gene was used as an internal control for analyzing gene expression. The qRT-PCR was completed with three biological replicates and three technical replicates.

Genome-Wide Association Study

GWAS was performed for each trait using a mixed linear model (MLM) in TASSEL 5 program. The SNPs detected in more than 80% of accessions and MAF > 0.05 were used. The Linkage Disequilibrium (LD) between each pair of SNPs and the LD decay values were calculated using PopLDdecay.⁷ More details about GWAS were described in previous reports (Sun et al., 2016).

Bioinformatics Tools

The information of main Bioinformatics tools used in this study is listed in Supplementary Table 4.

Results

Overview of the Genetic Variations in 210 Rapeseed Accessions

To investigate the genetic variations in rapeseed, we collected 210 elite rapeseed accessions exhibiting 12 elite agronomic traits from worldwide collections, including 164 inbred lines, 20 doubled haploid lines, 14 near isogenic lines, 4 ethyl methanesulfonate-mutated lines, and 8 self-selected lines (Supplementary Table 1 and Figure 1A). All elite accessions had one or more desirable agronomic trait and were derived from interspecific crosses, mutagenesis, and natural variations (Supplementary Table 1). The analysis of re-sequencing data confirmed that almost all accessions were homozygous and suitable for further study.

FIGURE 1

Figure 1. Summary of the 210 rapeseed accessions exhibiting 12 important agronomic traits and the pipeline for analyzing genetic variations. (A) Percentages of different rapeseed categories. (B) Pipeline for the DNA re-sequencing and correlation analysis. Genome REF refers to the B. napus “Darmor” reference genome, whereas transcript REF refers to the B. napus transcript reference sequences; both were downloaded from the BRAD website. The number in each frame represents the sum of the values. Cutadapter, FastQC, BWA, SortSam, Picard, GATK, and psRNATarget were the main bioinformatic tools used in the pipeline, whereas TAIR, BRAD, and miRBase were the key databases used in this study. (C) Schematic representation of the 12 important agronomic traits. The arrows indicate a life cycle.

The DNA re-sequencing analysis of 210 rapeseed accessions generated 2.94 Tb of raw data. After a quality control step, 14 billion high-quality clean reads were obtained. The number of clean reads per accession ranged from 56 to 86 million, with a mean of 65 million clean reads. All clean reads from each accession were aligned to the B. napus “Darmor” variety reference genome (version 4.1). The mapping rate ranged from 83 to 99%, with a mean of 98%, and the effective coverage depth for each line ranged from 9- to 25-fold, with a mean of 12-fold (Supplementary Table 1). Duplicated reads caused by the PCR amplification were filtered and the remaining non-redundant reads were used for genotyping. The number of non-redundant reads for each accession ranged from 40 to 74 million, with a mean of 57 million reads (Supplementary Table 5 and Figure 1B). On the basis of the alignment with the reference genome, 7.15 million SNPs and 2.89 million insertions/deletions (INDELs) were identified in the 210 rapeseed accessions. These variants were distributed on 19 chromosomes and scaffolds. Additionally, 2.01 million SNPs and 1.63 million INDELs were rare variants in these rapeseed accessions [minor allele frequency (MAF) ≤ 0.05, integrity ≥ 0.7]. Approximately 40.22% of the variants were rare events.

We subsequently annotated the genomic variants in rapeseed using ANNOVAR (Supplementary Table 6). Most of the variants were in intergenic regions (52.94% of the SNPs and 43.98% of the INDELs). Regarding the genic regions, 10.83% of the SNPs and 20.37% of the INDELs were in introns, 15.86% of the SNPs and 3.84% of the INDELs were in exons, 1.36% of the SNPs and 2.60% of the INDELs were in 3′-UTRs, and 0.85% of the SNPs and 1.94% of the INDELs were in 5′-UTRs. Moreover, 8.87% of the SNPs and 13.57% of the INDELs were in the 1 kb region upstream of the transcription start site, whereas 7.51% of the SNPs and 10.20% of the INDELs were in the 1 kb region downstream of the transcription termination site. Accordingly, these variants may affect transcription.

Prediction of miRNA Target Sites in Rapeseed Transcripts

In general, miRNAs negatively regulate gene expression by targeting the 3′-UTR or the coding region of mRNAs. Therefore, variants in exons and UTRs might significantly affect gene expression via miRNA binding. To reveal the genetic variations in miRNA target sites in rapeseed, we first identified the gene regions bound by miRNAs using psRNATarget with cruciferous miRNAs, with a maximum expectation value of 3.5 (Jones-Rhoades and Bartel, 2004; Alves et al., 2009; Yu et al., 2011; Dai et al., 2018; He J. et al., 2018; Zhang X. D. et al., 2018). A total of 163 crucifer miRNAs from four species (B. napus, Brassica rapa, Brassica oleracea, and Arabidopsis thaliana) were retrieved from the miRbase database (version 22)⁸ (Supplementary Table 7). The protein-coding transcripts were downloaded from BRAD (see text footnote 1). Together with previously reported degradome data and the results of the 5′ RNA ligase-mediated rapid amplification of cDNA ends (Lu et al., 2008; Vaucheret, 2009; Chen et al., 2011; Ng et al., 2011; Kinoshita et al., 2012; Bari et al., 2013; Song et al., 2013; Zhang et al., 2013, 2015; Liu H. et al., 2014; Park et al., 2014; Windels et al., 2014; Xian et al., 2014; Ma et al., 2015, 2017; Ai et al., 2016; Ouyang et al., 2016; Leng et al., 2017; He F. et al., 2018; Kuo et al., 2018; Tang et al., 2018; Zheng et al., 2019), we identified 1,644 high-confidence target site sequences (Supplementary Table 8). After being generated from stem-loop precursors by the DCL1 complex, the mature miRNAs are loaded into an AGO1-dominant RNA-induced silencing complex (RISC) and repress the target transcripts mainly in two ways: cleavage and translational repression (Brodersen et al., 2008). For the miRNA target sites, 1284 were predicted to be cleaved and 360 could repress translation by the AGO1. In addition, 97% of the target site sequences were 20–22 nucleotides long, with the most (77%) was 21 nucleotides long (Figure 2A).

FIGURE 2

Figure 2. Features of the variants in miRNA target sites. (A) Length distribution of putative miRNA target sites. (B) Percentages of the A, G, C, and T variations. (C) Percentages of the different types of nucleotide substitutions among the SNPs. C↔T: C to T or T to A substitution; G↔A: G to A or A to G substitution; C↔T: G to T or T to G substitution; A↔C: A to C or C to A substitution; A↔T: A to T or T to A substitution; C↔G: C to G or G to A substitution. (D) Percentage of SNPs in each miRNA target site. (E) Annotation of the variants identified in miRNA target sites. (F) Number of rare and common variants in miRNA target sites. (G) Altered target binding caused by genetic variants within miRNA target sites. UPE: target accessibility energy to unpair the target site (kcal/mol); decreases in energy increase the possibility that small RNAs will be able to interact with (and cleave) target mRNAs. Allelic SNPs are presented in blue and are marked by red arrows. UPEs are calculated using psRNAtarget.

To functionally characterize the miRNA target genes in rapeseed, we annotated the protein-coding genes in B. napus via a BLASTN search of the Araport11 genome release in The Arabidopsis Information Resource (TAIR) database.⁹ The annotation of rapeseed genes is relatively limited, but the functions of these genes may be determined by analyzing the similarity between the mRNA sequences of rapeseed and A. thaliana (85% homology) (Cavell et al., 1998). In total, 101,040 protein-coding genes have been reported in rapeseed (Chalhoub et al., 2014). Here, we found 67,288 genes were homologous to A. thaliana genes and 33,752 genes were specific to rapeseed (Supplementary Figure 1A and Supplementary Table 9). Among the 1,644 identified miRNA target genes, 1,394 were homologous to A. thaliana genes, whereas 250 were specific to rapeseed (Supplementary Figure 1B).

Next, we performed a gene ontology (GO) enrichment analysis of the miRNA target genes homologous to A. thaliana genes. The overrepresented GO terms in the biological process category were response to stimulus, metabolic process, and biological regulation (Figure 3A). The main GO terms in the molecular function category were adenyl nucleotide and ribonucleotide binding and ATP binding.

FIGURE 3

Figure 3. Gene ontology (GO) analyze the miRNA target genes. (A) The top 14 GO terms (biological process category) among the miRNA target genes. The enriched processes are listed on the right. (B) The top 10 GO terms (biological process category) among the 328 variant miRNA target genes in B. napus. The enriched processes are listed on the left. (C) The top 10 GO terms (biological process category) among the 124 rare variant target genes in rapeseed. The enriched processes are listed on the right. Rich factor: number of candidate genes/total number of genes. BP: biological process.

Genetic Variations in the miRNA Target Sites in Rapeseed Accessions

To identify the genetic variations in the miRNA target sites in the examined rapeseed accessions, we extracted the variants from the genotype data. We detected 310 raw SNPs and 79 INDELs in the miRNA target sites (Supplementary Table 10). A global survey of SNPs in the miRNA target sites (integrity ≥ 0.7) revealed that 24.68, 26.60, 21.79, and 26.92% of the single-base variations resulted in an A, G, C, and T, respectively (Figure 2B). The C↔T (U) and A↔G substitutions were the most common variations, which is consistent with the fact that U/G at an mRNA site can still partially pair with G/U of a mature miRNA (Figure 2C). The SNPs were more frequent at the first, tenth, and fifteenth nucleotides of the target site than at the other nucleotides. The SNPs were the least frequently detected at the twelfth and twenty-third nucleotides (Figure 2D). To clarify the potential effects of these variants, the detected SNPs and INDELs (integrity ≥ 0.7) were functionally annotated by ANNOVAR according to the “Darmor” gene model. Most of the SNPs (80.12%) were in exons and 11.93% were in 3′-UTRs (Supplementary Table 6 and Figure 2E). This result was consistent with miRNAs targeting the 3′-UTR or the coding region of mRNAs. These SNPs were linked to 257 miRNA target genes, of which 214 were homologs of A. thaliana genes and 43 were specific to rapeseed (Supplementary Figure 1C). Most of the INDELs (94.05%) were in exons and UTRs. More specifically, 38.10, 28.57, and 27.38% of the INDELs were in the 3′-UTRs, exons, and 5′-UTRs, respectively. The most common GO term among these miRNA target genes with variants (integrity ≥ 0.7) was single-organism process and cellular process (Figure 3B).

We next analyzed the rare variants (MAF ≤ 0.05, integrity ≥ 0.7) in the miRNA target sites. A total of 106 rare SNPs and 30 rare INDELs were detected (Figure 2F). Approximately 33.50% of the variants were rare. Comparison of the rare variant ratio between miRNA target site and whole genome, we found miRNA target site had less rare variants. The 106 rare SNPs and 30 rare INDELs were linked to 124 miRNA target genes. The most commonly assigned GO term among these rare variant target genes was single-organism process (P < 0.05) (Figure 3C).

We subsequently focused on the SNPs affecting the complementation between miRNAs and their target genes. For example, the SNP A10C (A to C at the tenth nucleotide position) in the miR156–SPL5 target site increased the target accessibility energy to unpair the target site (UPE) (Figure 2G); decreases in energy increase the likelihood that small RNAs will be able to interact with (and cleave) target mRNAs (Dai et al., 2018). Therefore, the SNP A10C might hinder the ability of miR156 to cleave SPL5 in plants. The presence of two or more variants in one miRNA target site is another type of variation. Both G7A and T13C were within the miR5719-NLR target site (Figure 2G), which decreased the UPE from 14.456 to 12.499, thereby increasing the ability of miR5719 to bind and cleave NLR genes. Both SPL5 and NLR genes are reportedly important for plant development and disease resistance (Ferreira e Silva et al., 2014; Shao et al., 2019; Jiang et al., 2020). Overall, specific variants that impair the binding of miRNAs, including A10C, G7A, and T13C, affect the expression and biological function of target genes.

Phenotypes of 12 Important Agronomic Traits in Rapeseed

To reveal the possible biological effects of the variants in the miRNA target sites, we analyzed the correlation between the variants and rapeseed agronomic traits. We examined the following 12 important agronomic traits in 210 rapeseed accessions: plant height, oil content, branch angle, branch number, silique number, seed number, seed weight, ripening time, flowering time, cold resistance, stem lodging resistance, and S. sclerotiorum resistance (Figure 1C). The 210 rapeseed accessions were grown in the field in 2016, 2017, and 2018. The phenotypic variations among the accessions are summarized in Table 1. We then used one-way ANOVA to analyze the environmental effect variance to phenotypic variance. We found that seven traits over the 3 years showed no significant difference and the other five traits (branch angle, branch number, seed number, seed weight and S. sclerotiorum resistance) had at least 1 year showed significant difference (Supplementary Table 11). These results suggested that environmental factors in different years affected the above five traits. All phenotypic data collected in different environments had a normal or skewed distribution (Supplementary Figures 2A–C), suggesting the traits were quantitatively inherited.

TABLE 1

Table 1. Phenotypic variations analysis among 210 B. napus accessions in the correlation panel.

Next, we determined the pairwise correlations among the 12 traits (Supplementary Table 12). Silique number and branch number were the most positively correlated traits, whereas seed number and seed weight had the weakest negative correlation. The results of the correlation analysis imply that most of the traits were not independent of each other.

Correlations Between Agronomic Traits and SNPs in miRNA Target Sites

To examine the correlations between the genetic variations in target sites and the phenotypic variance, we selected 199 high-quality SNPs (MAF > 0.05, integrity ≥ 0.7) in miRNA target sites for further analyses. We used the general linear model to clarify the correlations between the SNPs in miRNA target sites and the 12 important agronomic traits. A P value ≤ 0.005 (1/199 SNPs) was set as the threshold for the correlation study.

The single-SNP correlation approach revealed 223 (86 in 2016, 79 in 2017, and 58 in 2018) significant correlations between miRNA target sites and the examined phenotypes (P ≤ 0.005) (Supplementary Table 13). These significant SNPs partially explained the observed phenotypic variance (average R² = 6.68, 6.69, and 7.34% in 2016, 2017, and 2018, respectively). The 223 correlated pairs involved 103 SNPs in 91 miRNA target genes. The relevant miRNA target genes included multiple transcription factor genes from the MYB, NAC, SPL, NF-YA, and NBS-LRR families. These transcription factors are related to various biological processes (e.g., growth and development and responses to abiotic or biotic stresses) in cells during the life cycle (He and Hannon, 2004; Ferreira e Silva et al., 2014; Si et al., 2016; Liu et al., 2020). These results indicate that the predicted candidate genes as well as the corresponding miRNAs might function together to regulate rapeseed growth and development.

Eleven significant SNP–trait correlations were consistently detected in the three consecutive years included in this study. According to the (contribution rate) R² values, three SNPs (chrC04__46917368, chrA03__16353138, and chrC03__8380169) explained the phenotypic variance more than the other SNPs (Table 2). More specifically, chrC04__46917368, which was within the miR156–SPL9 target site, explained 12.57–17.51% of the phenotypic variance in the three consecutive years. The SNP chrA03__16353138, which was correlated with the silique number and was located in the miR164–CUC1 target site, explained 10.28–13.06% of the phenotypic variance. The SNP chrC03__8380169, which was correlated with the seed weight and was detected in the miR156–SPL13 target site, explained 12.01–23.57% of the phenotypic variance. Considering they were detected in multiple years with high R² values, the miRNA target genes SPL9, CUC1, and SPL13 as well as the related miRNAs may contribute to branch, silique, and seed developmental processes, respectively.

TABLE 2

Table 2. Co-detected SNPs significantly correlated with specific agronomic traits in three consecutive years.

Also, we investigated the linkage disequilibrium (LD) of each pair of SNPs, and performed GWAS for branch angle, silique number and seed weight using the mixed linear model (MLM), considering the genetic background, and found that the SNP of chrC04__46917368 (variation in BnSPL9-3 targeted by miR156), chrA03__16353138 (variation in BnCUC1-2 targeted by miR164) and chrC03__8380169 (variation in BnSPL13-2 targeted by miR156) were in the loci associated with branch angle, silique number and seed weight, respectively (Supplementary Figure 3).

Functional Analysis of the SNPs in miRNA Target Sites

A previous study proved that mutations in miRNA target sites in A. thaliana change the RNA abundance of the corresponding target genes (Mallory et al., 2004). To verify the effects of SNPs in the miRNA target sites on rapeseed agronomic traits, we first designed SNP-specific primers to validate the variants in miRNA target sites. The sequencing of the PCR products confirmed the authenticity of 98% of the target site variants. The three SNPs mentioned above were selected for further analyses. Consistent with the genotyping data, chrC04__46917368 was a T to G substitution (T9G) at the ninth target site nucleotide of BnSPL9-3 (BnaC04g48150D) (Figure 4D). In contrast, chrA03__16353138 was an A to C substitution (A1C) at the first target site nucleotide of BnCUC1-2 (BnaA03g33770D) (Figure 4H). The SNP chrC03__8380169 was an A to G substitution (A20G) at the twentieth target site nucleotide of BnSPL13-2 (BnaC03g16490D) (Figure 4L).

FIGURE 4

Figure 4. Haplotypes of candidate SNPs for the branch angle, silique number, and seed weight in rapeseed. Boxplots of the branch angle, silique number, and seed weight according to the SNP genotyping of BnSPL9-3 (A–C), BnCUC1-2 (E–G), and BnSPL13-2 (I–K). (D,H,L) Validation of variants in the candidate genes by Sanger sequencing. The variants are indicated by a red arrow. REF represents the reference sequence. Differences between the genotypes were analyzed using the Wilcoxon rank-sum test R package; “n” represents the number of accessions, whereas “a” and “b” represent significant differences at the 5% level.

We analyzed the effects of these three SNPs on agronomic traits in the 210 rapeseed accessions. The accessions with T9G in BnSPL9-3 had smaller branch angles than the accessions with T9 (Figures 4A–C). The BnCUC1-2 gene is a homolog of AtCUC1. An earlier study proved that knocking down AtCUC expression in A. thaliana affects silique development and results in abnormal pods (Liu et al., 2020). The accessions with A1C in BnCUC1-2 produced fewer siliques than the accessions with A1 (Figures 4E–G). The BnSPL13-2 gene is homologous to OsSPL13, which regulates seed size in rice (Si et al., 2016). The accessions with A20G in BnSPL13-2 had a lower seed weight than the accessions with A20 (Figures 4I–K). According to the correlation analysis, these three SNPs were correlated with the observed phenotypic differences.

Finally, to determine whether the SNPs in miRNA target sites affect gene silencing, we analyzed the expression levels of the miRNA-targeted genes in two groups of accessions with extremely different phenotypes by quantitative real-time PCR (qRT-PCR) (Figures 5A–C). Each group contained four accessions (Supplementary Table 3). The qRT-PCR data indicated the BnSPL9-3 expression level was significantly higher in the accessions with T9G and the smallest branch angle than in the accessions with T9 and the largest branch angle (Figure 5D). The BnCUC1-2 expression level was significantly lower in the accessions with A1C and the fewest siliques than in the accessions with A1 and the most siliques (Figure 5E). The BnSPL13-2 expression level was significantly lower in the accessions with A20G and the lowest seed weight than in the accessions with A20 and the highest seed weight (Figure 5F). These results demonstrated that the three SNPs were functional and that a single SNP within the miRNA target sites was able to change target gene expression levels.

FIGURE 5

Figure 5. Effects of SNPs in the miRNA target site on target gene expression levels. (A–C) Effects of miRNA target site SNPs on the branch angle (A), silique number (B), and seed weight (C) determined by comparing groups with extremely high values (four accessions) and low values (four accessions) for these traits. (D–F) Comparison of the BnSPL9-3 (D), BnCUC1-2 (E), and BnSPL13-2 (F) expression levels between the two groups with extremely different phenotypes. ACTIN7 was used as an internal control for analyzing gene expression. Data are presented as the mean ± standard deviation (n = 3). Significant differences were determined by Student’s t-test: *P < 0.05, **P < 0.01.

Discussion

Mutations are a source of genetic variations. In many cases, genetic variations are related to phenotypic diversity. There are substantial genetic variations among crop populations. Additionally, previous studies revealed that sequence variations in mature and precursor miRNA sequences affect miRNA biogenesis (Ryan et al., 2010), but very little is known about the effects of genetic variations at target sites, especially in rapeseed. Our rapeseed population included 210 elite accessions collected from the major rapeseed-producing regions worldwide. All accessions were homozygous and suitable for genotyping and phenotyping. We re-sequenced the rapeseed population and detected millions of variants, of which 327 SNPs and 79 INDELs were in 357 miRNA-targeted sites. Hence, the genetic variation in the target sites is apparently extensive. Among the identified target site variants, most of the SNPs (80.12%) were in exons and 16.83% were located within UTRs (Supplementary Table 6), reflecting the fact the miRNAs targeted the 3′-UTR or the coding region of mRNAs.

Association analyses generally use common variants (MAF > 0.05) to detect significantly associated variants. The rare variants are eliminated from the genotypic data. However, in this study, we revealed that approximately 40.22% of the variants (SNPs and INDELs) in the whole genome were rare and 33.50% of the variants in the miRNA target sites were rare (MAF ≤ 0.05). The rare variants in a population have important but unknown roles influencing plant growth. The GO term enrichment analysis of the genes with these rare variants indicated that most were related to single-organism processes (Figure 3C). How the rare variants modulate phenotypes remains to be investigated.

Genetic variations at target sites can alter miRNA–mRNA interactions and affect the expression of miRNA target genes (Smith et al., 2015). The complementarity and thermodynamics of binding are critical factors related to the interaction between miRNAs and their target genes (Gong et al., 2012; Smith et al., 2015). Therefore, sequence variants, such as SNPs in the seed region or even in the cleavage site, may alter miRNA–mRNA interactions and affect the expression of miRNA target genes. Compared with animal miRNAs and their targets, the complementarity between plant miRNAs and their targets is greater and there are fewer mismatches in the binding region. It is possible that a SNP in the miRNA-targeted regions will disrupt gene regulation (Liu N. et al., 2014). The SNP A10C within the miR156–SPL5 target site and the SNPs G7A and T13C within the miR5719–NLR target site were predicted to change the accessibility of the miRNA targets (Figure 2G). Accordingly, these SNPs might affect the expression and biological function of the target genes. Compared with prior studies that often focused on the change of protein coding genes, the large-scale detection of variants in miRNA target site could provide new clues for elucidating the miRNA-dependent gene expression regulation.

In this study, a SNP-based correlation analysis revealed 103 independent genetic variants, 223 significantly correlated pairs, and 91 miRNA target genes. Moreover, 11 SNP–trait correlations were detected in three consecutive years. Based on genome-wide association study using mixed linear model, we found the variation of miRNA target sites at BnSPL9-3, BnCUC1-2, and BnSPL13-2 were in the loci associated with their correlated traits, respectively. However, due to the large LD block, these variations in miRNA target sites were not the only candidates in the trait loci, fine mapping were needed for further investigation. Notably, BnSPL9-3, BnCUC1-2, and BnSPL13-2 expression levels varied between the accessions exhibiting extremely different phenotypes. A previous study demonstrated that the overexpression of zma-miR164e in transgenic A. thaliana downregulates CUC1 expression (Liu et al., 2009). The consequences of this altered CUC1 expression included the production of thin and short pods and the inability to produce normal pods and seeds, which ultimately seriously affected the number of pods per plant. The BnSPL13-2 gene is a homolog of OsSPL13, which regulates seed size in rice (Si et al., 2016). These earlier studies further supported that BnSPL13-2 may play the conserved role in controlling seed size. Therefore, these genes are worth considering to be used in crop improvement.

The miRNA target genes can also be miRNA biogenesis genes. For example, the AGO1 protein binds to miRNAs that post-transcriptionally repress target genes via the cleavage of mRNA and/or translational repression. The AGO1 mRNA is targeted by miR168. Thus, variants in target sites might have important effects on miRNA biogenesis. A recent study in B. rapa indicated that a single-base variation (T to C) can alter the leaf hair phenotype (Zhang X. D. et al., 2018). Such SNPs are functional and can be used for designing molecular markers relevant for breeding (Zhang W. et al., 2018). In the current study, we identified a set of SNPs in miRNA target sites. Additionally, three of the SNPs (T9G, A1C, and A20G) in the target sites led to significant changes to the correlated traits in the rapeseed population (Figure 5). These SNPs may be candidate functional SNPs. In addition to the B. napus population, we also analyzed a B. rapa population (300 accessions) and a B. oleracea population (101 accessions) (NCBI accession number: PRJNA312457). Some SNPs in rapeseed miRNA target sites also exist in B. oleracea or B. rapa homologs, implying some of these genetic variations were derived from B. oleracea or B. rapa. Therefore, these SNPs appear to be conserved among Brassica crops and may contribute to the elite traits of B. oleracea or B. rapa.

Mutations in miRNA-targeted genes influence the domestication and selection of high-yielding and high-quality cultivars. A large branch angle is conducive to robust photosynthesis because of the increased exposure to sunlight. Increases in silique number and seed weight are crucial factors for high yields. Crop breeders have focused on the silique number and seed weight to improve their varieties. On the basis of the identification of SNPs in miRNA target sites, future investigations should explore the possibility of introducing specific mutations in miRNA target genes to enhance crop traits. This study generated a substantial amount of genotyping and phenotyping data for 210 rapeseed accessions, revealed the considerable genetic variations in miRNA target sites in rapeseed, and produced information regarding the correlations between SNPs in miRNA target sites and elite agronomic traits.

Conclusion

In this study, through re-sequenced 210 elite rapeseed accessions, we identified large amount of variants in genome. Furthermore, we detected 330 SNPs and 79 INDELs in 357 miRNA target sites, some of which were rare variants. These variants possibly affect phenotypic variation. Three SNPs within the miRNA-binding regions of BnSPL9-3, BnSPL13-2, and BnCUC1-2 were highly related to the branch angle, seed weight, and silique number, respectively. This outcome indicates that variation in the target sites is apparently extensive and distinctive. The genome-wide genetic variations in miRNA-targeted sites provide researchers and breeders a choice to detecting the allelic variation in miRNA target site related to elite traits.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: PRJNA748869).

Author Contributions

PX conceived the project. PX and YH designed the experiments. PX and YtZ performed the experiments. YfZ, JJ, LY, and JM analyzed and interpreted the data. PX, XY, and YH wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the National Programs for Science and Technology Development of China (Grant No. 2016YFD0101900) and the Natural Science Foundation of China (Grant Nos. 31471883 and 31571261).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.741858/full#supplementary-material

Supplementary Figure 1 | Comparison of protein-coding gene and miRNA target site sequences in B. napus and A. thaliana. (A) Comparison of protein-coding genes between B. napus and A. thaliana. (B) Putative target sites. (C) Target genes with sequence variations in B. napus.

Supplementary Figure 2 | Frequency distribution of the 12 examined agronomic traits in rapeseed in three consecutive years. (A–C) Frequency distribution of 12 traits in 2016 (A), 2017 (B), and 2018 (C). PH: plant height; OC: oil content; BA: branch angle; BN: branch number; SNPP: silique number; SNPS: seed number; SW: seed weight; RT: ripening time; FT: flowering time; RC: cold resistance; RSL: stem lodging resistance; RS: Sclerotinia sclerotiorum resistance.

Supplementary Figure 3 | Genome-wide association studies of branch angle, silique number and seed weight. Local Manhattan plot (top) and LD heatmap (bottom) surrounding the SNP. Arrow indicates the position of chrC04__46917368 (A), chrA03__16353138 (B), and chrC03__8380169 (C). The local Manhattan was shown of the SNPs in flanking region (−1 Mb to +1 Mb) of the targeted sites. LD heatmap of SNPs in the qSPL9 (D), qCUC1 (E), and qSPL13 (F) loci, respectively. The local LD was calculated of the SNPs in flanking region (−1 kb to +1kb) of the targeted sites.

Supplementary Table 1 | Summary of all accessions sequenced in this study.

Supplementary Table 2 | Primers used in this study.

Supplementary Table 3 | Accessions analyzed by qRT-PCR.

Supplementary Table 4 | Software and websites used in this study.

Supplementary Table 5 | Number of non-redundant reads for each accession.

Supplementary Table 6 | Annotation of the SNPs and INDELs identified in B. napus.

Supplementary Table 7 | Crucifer family miRNA sequences.

Supplementary Table 8 | Putative miRNA target site sequences in rapeseed transcripts.

Supplementary Table 9 | Global annotation of protein-coding genes in B. napus. The BLAST algorithm was used to compare the Arabidopsis thaliana transcripts with B. napus transcripts (threshold: e-50).

Supplementary Table 10 | SNPs and INDELs in miRNA target sites in rapeseed accessions.

Supplementary Table 11 | Analysis of variance of the 12 observed agronomic traits. Significant differences were determined by F-test (ANOVA) using the SAS software. The letter “A” and “B” represent significant differences at the 1% level.

Supplementary Table 12 | Correlations among the 12 observed agronomic traits. The correlation analysis of the different traits was performed using the cor.test function in R (version 3.6.1). Significance levels: ^∗P < 0.05; ^∗∗P < 0.01; ^∗∗∗P < 0.001.

Supplementary Table 13 | Summary of significantly correlated SNPs.

Footnotes

1. ^ http://brassicadb.cn/
2. ^ https://github.com/broadinstitute/gatk
3. ^ https://hpc.nih.gov/apps/ANNOVAR.html
4. ^ https://github.com/NBISweden/AGAT
5. ^ https://www.zhaolab.org/psRNATarget/
6. ^ https://www.maizegenetics.net/tassel
7. ^ https://github.com/BGI-shenzhen/PopLDdecay
8. ^ http://www.mirbase.org/
9. ^ https://www.arabidopsis.org

References

Ai, Q., Liang, G., Zhang, H., and Yu, D. (2016). Control of sulfate concentration by miR395-targeted APS genes in Arabidopsis thaliana. Plant Divers. 38, 92–100. doi: 10.1016/j.pld.2015.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Allender, C. J., and King, G. J. (2010). Origins of the amphiploid species Brassica napus L. investigated by chloroplast and nuclear molecular markers. BMC Plant Biol. 10:54. doi: 10.1186/1471-2229-10-54

PubMed Abstract | CrossRef Full Text | Google Scholar

Alves, L. J., Niemeier, S., Hauenschild, A., Rehmsmeier, M., and Merkle, T. (2009). Comprehensive prediction of novel microRNA targets in Arabidopsis thaliana. Nucleic Acids Res. 37, 4010–4021. doi: 10.1093/nar/gkp272

PubMed Abstract | CrossRef Full Text | Google Scholar

Bari, A., Orazova, S., and Ivashchenko, A. (2013). miR156- and miR171-Binding sites in the protein-coding sequences of several plant genes. Biomed Res Int. 2013:307145. doi: 10.1155/2013/307145

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartel, D. P. (2004). MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116, 281–297. doi: 10.1016/S0092-8674(04)00045-5

CrossRef Full Text | Google Scholar

Brodersen, P., Sakvarelidze-Achard, L., Bruun-Rasmussen, M., Dunoyer, P., Yamamoto, Y. Y., Sieburth, L., et al. (2008). Widespread translational inhibition by plant miRNAs and siRNAs. Science 320, 1185–1190. doi: 10.1126/science.1159151

PubMed Abstract | CrossRef Full Text | Google Scholar

Cavell, A. C., Lydiate, D. J., Parkin, I. A. P., Dean, C., and Trick, M. (1998). Collinearity between a 30-centimorgan segment of Arabidopsis thaliana chromosome 4 and duplicated regions within the Brassica Napus genome. Genome 41, 62–69. doi: 10.1139/g97-097

PubMed Abstract | CrossRef Full Text | Google Scholar

Chalhoub, B., Denoeud, F., Liu, S., Parkin, I. A. P., Tang, H., Wang, X., et al. (2014). Early allopolyploid evolution in the post-Neolithic Brassica Napus oilseed genome. Science 345, 950–953. doi: 10.1126/science.1253435

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Xie, J., Chen, B., Quan, M., Li, Y., Li, B., et al. (2016). Genetic variations and miRNA-target interactions contribute to natural phenotypic variations in Populus. New Phytol. 212, 150–160. doi: 10.1111/nph.14040

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Z., Bao, M., Sun, Y., Yang, Y., Xu, X., Wang, J., et al. (2011). Regulation of auxin response by miR393-targeted transport inhibitor response protein1 is involved in normal development in Arabidopsis. Plant Mol. Biol. 77, 619–629. doi: 10.1007/s11103-011-9838-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Chuck, G., Candela, H., and Hake, S. (2009). Big impacts by small RNAs in plant development. Curr. Opin. Plant Biol. 12, 81–86. doi: 10.1016/j.pbi.2008.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Dai, X., Zhuang, Z., and Zhao, P. X. (2018). psRNATarget: a plant small RNA target analysis server (2017 release). Nucleic Acids Res. 46, W49–W54. doi: 10.1093/nar/gky316

PubMed Abstract | CrossRef Full Text | Google Scholar

Dugas, D. V., and Bartel, B. (2004). MicroRNA regulation of gene expression in plants. Curr. Opin. Plant Biol. 7, 512–520. doi: 10.1016/j.pbi.2004.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Ferreira e Silva, G. F., Silva, E. M., Da Silva Azevedo, M., Guivin, M. A. C., Ramiro, D. A., Figueiredo, C. R., et al. (2014). microRNA156-targeted SPL/SBP box transcription factors regulate tomato ovary and fruit development. Plant J. 78, 604–618. doi: 10.1111/tpj.12493

PubMed Abstract | CrossRef Full Text | Google Scholar

Gong, J., Tong, Y., Zhang, H., Wang, K., Hu, T., Shan, G., et al. (2012). Genome-wide identification of SNPs in microRNA genes and the SNP effects on microRNA target binding and biogenesis. Hum. Mutat. 33, 254–263. doi: 10.1002/humu.21641

PubMed Abstract | CrossRef Full Text | Google Scholar

He, F., Xu, C., Fu, X., Shen, Y., Guo, L., Leng, M., et al. (2018). The MicroRNA390 TRANS-ACTING SHORT INTERFERING RNA3 Module mediates lateral root growth under salt stress via the auxin pathway. Plant Physiol. 177, 775–791. doi: 10.1104/pp.17.01559

PubMed Abstract | CrossRef Full Text | Google Scholar

He, J., Xu, M., Willmann, M. R., McCormick, K., Hu, T., Yang, L., et al. (2018). Threshold-dependent repression of SPL gene expression by miR156/miR157 controls vegetative phase change in Arabidopsis thaliana. PLoS Genet. 14:e1007337. doi: 10.1371/journal.pgen.1007337

PubMed Abstract | CrossRef Full Text | Google Scholar

He, L., and Hannon, G. J. (2004). MicroRNAs: small RNAs with a big role in gene regulation. Nat. Rev. Genet. 5, 522–531. doi: 10.1038/nrg1379

PubMed Abstract | CrossRef Full Text | Google Scholar

Inui, M., Martello, G., and Piccolo, S. (2010). MicroRNA control of signal transduction. Nat. Rev. Mol. Cell Biol. 11, 252–263. doi: 10.1038/nrm2868

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, Y., Peng, J., Wang, M., Su, W., Gan, X., Jing, Y., et al. (2020). The Role of EjSPL3. EjSPL4, EjSPL5, and EjSPL9 in regulating flowering in loquat (Eriobotrya Japonica Lindl.). Int. J. Mol. Sci. 21:248. doi: 10.3390/ijms21010248

PubMed Abstract | CrossRef Full Text | Google Scholar

Jones-Rhoades, M. W., and Bartel, D. P. (2004). Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol. Cell 14, 787–799. doi: 10.1016/j.molcel.2004.05.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Kidner, C. A., and Martienssen, R. A. (2005). The developmental role of microRNA in plants. Curr. Opin. Plant Biol. 8, 38–44. doi: 10.1016/j.pbi.2004.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Kinoshita, N., Wang, H., Kasahara, H., Liu, J., MacPherson, C., Machida, Y., et al. (2012). IAA-Ala Resistant3, an evolutionarily conserved target of miR167, mediates arabidopsis root architecture changes during high osmotic stress. Plant Cell 24, 3590–3602. doi: 10.1105/tpc.112.097006

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuo, Y., Lin, J., Li, Y., Jhu, M., King, Y., and Jeng, S. (2018). MicroR408 regulates defense response upon wounding in sweet potato. J. Exp. Bot. 70, 469–483. doi: 10.1093/jxb/ery381

PubMed Abstract | CrossRef Full Text | Google Scholar

Leng, X., Wang, P., Zhu, X., Li, X., Zheng, T., Shangguan, L., et al. (2017). Ectopic expression of CSD1 and CSD2 targeting genes of miR398 in grapevine is associated with oxidative stress tolerance. Funct. Integr. Genomics. 17, 697–710. doi: 10.1007/s10142-017-0565-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Cao, X., Zhu, Y., Yang, X., Zhang, K., Xiao, Z., et al. (2019). Osa-miR398b boosts H2 O2 production and rice blast disease-resistance via multiple superoxide dismutases. New Phytol. 222, 1507–1522. doi: 10.1111/nph.15678

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H., Guo, S., Xu, Y., Li, C., Zhang, Z., Zhang, D., et al. (2014). OsmiR396d-Regulated OsGRFs function in floral organogenesis in rice through binding to their targets OsJMJ706and OsCR4. Plant Physiol. 165, 160–174. doi: 10.1104/pp.114.235564

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Cheng, X., Liu, P., and Sun, J. (2017). miR156-targeted SBP-box transcription factors interact with DWARF53 to regulate TEOSINTE BRANCHED1 and BARREN STALK1 expression in bread wheat. Plant Physiol. 174, 1931–1948. doi: 10.1104/pp.17.00445

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, M., Tan, X., Yang, Y., Liu, P., Zhang, X., Zhang, Y., et al. (2020). Analysis of the genetic architecture of maize kernel size traits by combined linkage and association mapping. Plant Biotechnol. J. 18, 207–221. doi: 10.1111/pbi.13188

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, N., Tu, L., Tang, W., Gao, W., Lindsey, K., and Zhang, X. (2014). Small RNA and degradome profiling reveals a role for miRNAs and their targets in the developing fibers of Gossypium barbadense. Plant J. 80, 331–344. doi: 10.1111/tpj.12636

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Q., Zhang, Y., Wang, C., Luo, Y., Huang, Q., Chen, S., et al. (2009). Expression analysis of phytohormone-regulated microRNAs in rice, implying their regulation roles in plant hormone signaling. FEBS Lett. 583, 723–728. doi: 10.1016/j.febslet.2009.01.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, T., Tang, J., Chen, L., Zeng, J., Wen, J., Yi, B., et al. (2019). Differential expression of miRNAs and their targets in wax-deficient rapeseed. Sci. Rep. 9:12201. doi: 10.1038/s41598-019-48439-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Jia, L., Mao, Y., and He, Y. (2010). Classification and quantification of leaf curvature. J. Exp. Bot. 61, 2757–2767. doi: 10.1093/jxb/erq111

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, X., Guan, Q., and Zhu, J. (2013). Downregulation of CSD2 by a heat-inducible miR398 is required for thermotolerance in Arabidopsis. Plant Signal. Behav. 8:e24952. doi: 10.4161/psb.24952

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y., Gan, Q., Chi, X., and Qin, S. (2008). Identification and characterization of MicroRNAs and Their Targets in grapevine (Vitis vinifera). Agric. Sci. China 7, 929–943. doi: 10.1016/S1671-2927(08)60132-7

CrossRef Full Text | Google Scholar

Ma, C., Burd, S., and Lers, A. (2015). miR408 is involved in abiotic stress responses in Arabidopsis. Plant J. 84, 169–187. doi: 10.1111/tpj.12999

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, Z., Jiang, J., Hu, Z., Lyu, T., Yang, Y., Jiang, J., et al. (2017). Over-expression of miR158 causes pollen abortion in Brassica campestris ssp. chinensis. Plant Mol. Biol. 93, 313–326. doi: 10.1007/s11103-016-0563-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Mallory, A. C., Reinhart, B. J., Jones-Rhoades, M. W., Tang, G., Zamore, P. D., Barton, M. K., et al. (2004). MicroRNA control of PHABULOSA in leaf development: importance of pairing to the microRNA 5′ region. EMBO J. 23, 3356–3364. doi: 10.1038/sj.emboj.7600340

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandáková, T., and Lysak, M. A. (2008). Chromosomal phylogeny and karyotype evolution in x=7 crucifer species (Brassicaceae). Plant Cell 20, 2559–2570. doi: 10.1105/tpc.108.062166

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, Y., Wu, F., Yu, X., Bai, J., Zhong, W., and He, Y. (2014). MicroRNA319a-targeted Brassica rapa ssp. pekinensis TCP genes modulate head shape in Chinese cabbage by differential cell division arrest in leaf regions. Plant Physiol. 164, 710–720. doi: 10.1104/pp.113.228007

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrell, P. L., Buckler, E. S., and Ross-Ibarra, J. (2012). Crop genomics: advances and applications. Nat. Rev. Genet. 13, 85–96. doi: 10.1038/nrg3097

PubMed Abstract | CrossRef Full Text | Google Scholar

Navarro, L., Dunoyer, P., Jay, F., Arnold, B., Dharmasiri, N., Estelle, M., et al. (2006). A plant miRNA contributes to antibacterial resistance by repressing auxin signaling. Science 312, 436–439. doi: 10.1126/science.1126088

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, D. W.-K., Zhang, C., Miller, M., Palmer, G., Whiteley, M., Tholl, D., et al. (2011). cis- and trans-Regulation of miR163 and target genes confers natural variation of secondary metabolites in two arabidopsis species and their allopolyploids. Plant Cell 23, 1729–1740. doi: 10.1105/tpc.111.083915

PubMed Abstract | CrossRef Full Text | Google Scholar

Ouyang, X., Hong, X., Zhao, X., Zhang, W., He, X., Ma, W., et al. (2016). Knock out of the PHOSPHATE 2 Gene TaPHO2-A1 improves phosphorus uptake and grain yield under low phosphorus conditions in common wheat. Sci. Rep. 6:29850. doi: 10.1038/srep29850

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, Y. J., Lee, H. J., Kwak, K. J., Lee, K., Hong, S. W., and Kang, H. (2014). MicroRNA400-Guided cleavage of pentatricopeptide repeat protein mrnas renders arabidopsis thaliana more susceptible to pathogenic bacteria and fungi. Plant Cell Physiol. 55, 1660–1668. doi: 10.1093/pcp/pcu096

PubMed Abstract | CrossRef Full Text | Google Scholar

Parkin, I. A. P., Gulden, S. M., Sharpe, A. G., Lukens, L., Trick, M., Osborn, T. C., et al. (2005). Segmental structure of the Brassica napus genome based on comparative analysis with Arabidopsis thaliana. Genetics 171, 765–781. doi: 10.1534/genetics.105.042093

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren, W., Wang, H., Bai, J., Wu, F., and He, Y. (2018). Association of microRNAs with types of leaf curvature in Brassica rapa. Front. Plant Sci. 9:73. doi: 10.3389/fpls.2018.00073

PubMed Abstract | CrossRef Full Text | Google Scholar

Ryan, B. M., Robles, A. I., and Harris, C. C. (2010). Genetic variation in microRNA networks: the implications for cancer research. Nat. Rev. Cancer 10, 389–402. doi: 10.1038/nrc2867

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, Z.-Q., Xue, J.-Y., Wang, Q., Wang, B., and Chen, J.-Q. (2019). Revisiting the origin of plant NBS-LRR genes. Trends Plant Sci. 24, 9–12. doi: 10.1016/j.tplants.2018.10.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Shukla, L. I., Chinnusamy, V., and Sunkar, R. (2008). The role of microRNAs and other endogenous small RNAs in plant stress responses. Biochim. Biophys. Acta 1779, 743–748. doi: 10.1016/j.bbagrm.2008.04.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Si, L., Chen, J., Huang, X., Gong, H., Luo, J., Hou, Q., et al. (2016). OsSPL13 controls grain size in cultivated rice. Nat Genet. 48, 447–456. doi: 10.1038/ng.3518

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, L. M., Burbano, H. A., Wang, X., Fitz, J., Wang, G., Ural-Blimke, Y., et al. (2015). Rapid divergence and high diversity of miRNAs and miRNA targets in the Camelineae. Plant J. 81, 597–610. doi: 10.1111/tpj.12754

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, J., Gao, S., Sun, D., Li, H., Shu, X., and Yang, Z. (2013). miR394 and LCR are involved in Arabidopsis salt and drought stress responses in an abscisic acid-dependent manner. BMC Plant Biol. 13:210. doi: 10.1186/1471-2229-13-210

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, C., Wang, B., Yan, L., Hu, K., Liu, S., Zhou, Y., et al. (2016). Genome-Wide association study provides insight into the genetic control of plant height in rapeseed (Brassica napus L.). Front. Plant Sci. 7:1102. doi: 10.3389/fpls.2016.01102

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, M., Bai, X., Niu, L.-J., Chai, X., Chen, M.-S., and Xu, Z.-F. (2018). miR172 Regulates both vegetative and reproductive development in the perennial woody plant Jatropha curcas. Plant Cell Physiol. 59, 2549–2563. doi: 10.1093/pcp/pcy175

PubMed Abstract | CrossRef Full Text | Google Scholar

Tripathi, R. K., Bregitzer, P., and Singh, J. (2018). Genome-wide analysis of the SPL/miR156 module and its interaction with the AP2/miR172 unit in barley. Sci. Rep. 8:7085. doi: 10.1038/s41598-018-25349-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaucheret, H. (2009). AGO1 homeostasis involves differential production of 21-nt and 22-nt miR168 species by MIR168a and MIR168b. PLoS One 4:e6442. doi: 10.1371/journal.pone.0006442

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Czech, B., and Weigel, D. (2009). miR156-regulated SPL transcription factors define an endogenous flowering pathway in Arabidopsis thaliana. Cell 138, 738–749. doi: 10.1016/j.cell.2009.06.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., and Zhang, Q. (2017). Boosting rice yield by fine-tuning SPL gene expression. Trends Plant Sci. 22, 643–646. doi: 10.1016/j.tplants.2017.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Wu, F., Bai, J., and He, Y. (2014). BrpSPL9 (Brassica rapa ssp. pekinensis SPL9) controls the earliness of heading time in Chinese cabbage. Plant Biotechnol. J. 12, 312–321. doi: 10.1111/pbi.12138

PubMed Abstract | CrossRef Full Text | Google Scholar

Windels, D., Bielewicz, D., Ebneter, M., Jarmolowski, A., Szweykowska-Kulinska, Z., and Vazquez, F. (2014). miR393 is required for production of proper auxin signalling outputs. PloS One 9:e95972. doi: 10.1371/journal.pone.0095972

PubMed Abstract | CrossRef Full Text | Google Scholar

Xian, Z., Huang, W., Yang, Y., Tang, N., Zhang, C., Ren, M., et al. (2014). miR168 influences phase transition, leaf epinasty, and fruit development via SlAGO1s in tomato. J. Exp. Bot. 65, 6655–6666. doi: 10.1093/jxb/eru387

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, X., Wang, H., Lu, Y., De Ruiter, M., Cariaso, M., Prins, M., et al. (2011). Identification of conserved and novel microRNAs that are responsive to heat stress in Brassica rapa. J. Exp. Bot. 63, 1025–1038. doi: 10.1093/jxb/err337

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeng, D., Tian, Z., Rao, Y., Dong, G., Yang, Y., Huang, L., et al. (2017). Rational design of high-yield and superior-quality rice. Nat. Plants 3:17031. doi: 10.1038/nplants.2017.31

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, C., Xian, Z., Huang, W., and Li, Z. (2015). Evidence for the biological function of miR403 in tomato development. Sci. Hortic. 197, 619–626. doi: 10.1016/j.scienta.2015.10.027

CrossRef Full Text | Google Scholar

Zhang, W., Mirlohi, S., Li, X., and He, Y. (2018). Identification of functional single-nucleotide polymorphisms affecting leaf hair number in Brassica rapa. Plant Physiol. 177, 490–503. doi: 10.1104/pp.18.00025

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, X. D., Sun, J. Y., You, Y. Y., Song, J. B., and Yang, Z. M. (2018). Identification of Cd-responsive RNA helicase genes and expression of a putative BnRH 24 mediated by miR158 in canola (Brassica napus). Ecotoxicol. Environ. Saf. 157, 159–168. doi: 10.1016/j.ecoenv.2018.03.081

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y.-C., Yu, Y., Wang, C.-Y., Li, Z.-Y., Liu, Q., Xu, J., et al. (2013). Overexpression of microRNA OsmiR397 improves rice yield by increasing grain size and promoting panicle branching. Nat. Biotechnol. 31, 848–852. doi: 10.1038/nbt.2646

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, L., Zhang, X., Zhang, H., Gu, Y., Huang, X., Huang, H., et al. (2019). The miR164-dependent regulatory pathway in developing maize seed. Mol. Genet. Genomics. 294, 501–517. doi: 10.1007/s00438-018-1524-4

PubMed Abstract | CrossRef Full Text | Google Scholar