Chao Zhang1,2,3†
Wei Chang1,2†
Xiaodong Li1,2
Bo Yang1,2
Liyuan Zhang1,2
Zhongchun Xiao1,2
Jiana Li1,2,4*
Kun Lu1,2,4*
- 1Chongqing Rapeseed Engineering Research Center, College of Agronomy and Biotechnology, Southwest University, Chongqing, China
- 2Academy of Agricultural Sciences, Southwest University, Chongqing, China
- 3Oil Research Institute of Guizhou Province, Guizhou Academy of Agricultural Sciences, Guiyang, China
- 4Engineering Research Center of South Upland Agriculture, Ministry of Education, Chongqing, China
Harvest index (HI), the ratio of harvested seed weight to total aboveground biomass weight, is an economically critical value reflecting the convergence of complex agronomic traits. HI values in rapeseed (Brassica napus) remain much lower than in other major crops, and the underlying regulatory network is largely unknown. In this study, we performed mRNA and small RNA sequencing to reveal the mechanisms shaping HI in B. napus during the seed-filling stage. A total of 8,410 differentially expressed genes (DEGs) between high-HI and low-HI accessions in four tissues (silique pericarp, seed, leaves, and stem) were identified. Combining with co-expression network, 72 gene modules were identified, and a key gene BnaSTY46 was found to participate in retarded establishment of photosynthetic capacity to influence HI. Further research found that the genes involved in circadian rhythms and response to stimulus may play important roles in HI and that their transcript levels were modulated by differentially expressed microRNAs (DEMs), and we identified 903 microRNAs (miRNAs), including 46 known miRNAs and 857 novel miRNAs. Furthermore, transporter activity-related genes were critical to enhancing HI in good cultivation environments. Of 903 miRNAs, we found that the bna-miR396–Bna.A06SRp34a/Bna.A01EMB3119 pair may control the seed development and the accumulation of storage compounds, thus contributing to higher HI. Our findings uncovered the underlying complex regulatory network behind HI and offer potential approaches to rapeseed improvement.
Introduction
Harvest index (HI), formerly known as coefficient of economics, refers to the ratio of the economic yield to the biological yield of the crop at harvest (Donald, 1962), which can also be understood as the ratio of seeds harvested from a plant to all the above-ground biomass produced. HI thus reflects the distribution of crop assimilates products in economic yield organs and vegetative organs, and also indicates the patency of crop photosynthate transport from “source” organs to “sink” organs (Bennett et al., 2011, 2012), and HI has clear agronomic implications. Since the first Green Revolution, the increase in the yield of major crops such as rice (Oryza sativa), wheat (Triticum aestivum), and barley (Hordeum vulgare) has been mainly due to the increase in the HI (Austin et al., 1980). The allotetraploid crop rapeseed (Brassica napus) is cultivated worldwide for oil, which is extracted from its oil-rich seeds. Luo et al. (2015) examined the B. napus genetic architecture of HI using 35,791 high-throughput single nucleotide polymorphisms (SNPs) genotyped by the Illumina Brassica SNP60 Bead Chip in an association panel with 155 accessions, and a total of nine SNPs on the C genome were identified to be significantly associated with HI; they explained 3.42% of the phenotypic variance in HI. Based on Brassica SNP60 Bead Chip, a natural population (NP) containing 520 materials was used to perform genome-wide association study (GWAS) for traits related to HI. Combined with transcriptomic sequencing (RNA-seq) of materials with high and low HI, candidate genes involved in photosynthesis, inflorescence, and silique development were identified (Lu et al., 2016). In addition, in B. napus, it was also found that the BnaA02.TB1 (BnaA02g14010D) that regulates lateral branch development and the BnaA02.GW2 (BnaA02g18890D) that regulates grain weight may be related to HI (Chao et al., 2019); BnaDwf.C9 (BnaC09g20450D) and BnaC04.BIL1 (BnaC04g41660D) can affect HI by regulating plant height (Wang X. D. et al., 2020; Yang et al., 2021). These previous studies are of great worth for helping us to elucidate the genetic mechanism of HI in B. napus. As seed yield is part of the HI numerator, HI increases with grain yield and has, therefore, received widespread attention in breeding programs (Hay, 1995; Sadras and Lawson, 2011; Beche et al., 2014). Unfortunately, the HI of rapeseed remains much lower than that of other major crops, such as rice, wheat, maize (Zea mays), and soybean (Glycine max) (Hay, 1995; Unkovich et al., 2010). The HI for rapeseed ranges from 0.05 to 0.42 (Lu et al., 2016), which is far below the theoretical biological limit for HI of ∼0.60 in grain crops (Foulkes et al., 2009), indicating the potential to further increase HI. Therefore, improving HI in rapeseed varieties is a major objective for breeders.
Green plant tissues such as leaves and silique pericarps (SP) constitute “source” organs that are photosynthetically active and provide photoassimilates to “sink” organs like grains and seeds via translocation (designated here as “flow”); HI reflects flow and the balance between source tissues and sink organs (Unkovich et al., 2010). In B. napus, flow is thought to be the limiting factor in the accumulation of assimilates in seeds (Shen et al., 2010; Fu and Zhou, 2013; Luo et al., 2015). Unlike many other crops in which leaves provide the main source tissue throughout seed development, in B. napus, SP take over the role of source tissue from senescent leaves, providing nutrients to sustain the growth of the seed after silique formation (Diepenbrock, 2000; Bennett et al., 2011). Several critical sugar transporter gene families have been identified in B. napus, such as SUCROSE TRANSPORTER (BnSUC), SUGARS WILL EVENTUALLY BE EXPORTED TRANSPORTER (BnSWEET), and MONOSACCHARIDE TRANSPORTER (BnMST); transporters encoded by these genes may influence HI positively during the development of siliques by supplying more sugars to the seeds (Jian et al., 2016; Zhang et al., 2020). In addition, the Arabidopsis thaliana AMINO ACID PERMEASE 2 (AAP2) gene encodes an amino acid co-transporter that transports amino acids to the embryo via the seed pod vascular system during the seed-filling phase (Hirner et al., 1998), suggesting that SP may supply not only sugars but also other compounds. However, the underlying molecular regulatory mechanism behind the translocation of photoassimilates and other nutrients during the seed-filling phase is poorly understood in B. napus.
Harvest index reflects the relative allocation of resources to vegetative and reproductive organs and can be influenced by various environmental factors, including water supply (Gajića et al., 2018), temperature stress (Prasad et al., 2006), and nitrogen (N) fertilization (Amanullah and Shah, 2010). For crop species such as maize, shortening the vegetative growth period reduces the accumulation of photosynthetic products in vegetative organs and raises HI (Hütsch and Schubert, 2017). One method employed to modulate the length of the vegetative period relies on the genetic manipulation of photoperiod sensitivity (Li et al., 2020). In addition, enhancing crop adaptability to external stimuli can also result in significant increases in harvestable products and thus contribute to a higher HI. For example, the rice basic leucine zipper transcription factor bZIP58 induces the expression of seed storage protein genes and starch biosynthetic genes, but its transcripts also undergo alternative splicing in response to heat stress. However, more heat-tolerant rice varieties showed limited alternative splicing and higher expression of seed storage protein genes under heat stress, thus contributing to heat tolerance during grain filling (Xu et al., 2020). In rice, a dominant mutation in the DENSE AND ERECT PANICLE 1 (DEP1) gene displayed insensitivity to N during vegetative growth, which increased N uptake and assimilation and thus improved HI and grain yield at moderate levels of N fertilization (Sun et al., 2014). However, which genes regulate HI by promoting reproductive growth and mediating responses to external stimuli remain largely unknown in B. napus.
MicroRNAs (miRNAs) are 20–24-nt single-stranded non-coding RNA molecules that act as key post-transcriptional regulators of gene expression (Voinnet, 2009; Budak and Akpinar, 2015). Extensive studies have revealed that miRNAs play diverse roles during plant development and adjust complex traits in crops (D’Ario et al., 2017; Tang and Chu, 2017). For example, Arabidopsis miR159 lowers the transcript levels of the MYB DOMAIN PROTEIN 33 (MYB33) gene, whose encoded transcription factor normally negatively regulates miR156, thereby modulating vegetative phase change (Guo et al., 2017). Similarly, miR164 induces cleavage of NAM, ATAF1/2, and CUC2 (NAC2) gene transcripts, which act as negative regulators of drought tolerance in rice (Fang et al., 2014). In Arabidopsis, miR397b regulates LACCASE 4 (LAC4) and influences silique numbers and silique length (Wang et al., 2014). Moreover, rice miR396d forms a molecular bridge between BRASSINAZOLE-RESISTANT 1 (BZR1) and GROWTH REGULATING FACTORs (GRFs); OsBZR1 induces the expression of miR396d, which, in turn, represses OsGRF transcript accumulation to modulate plant architecture and grain yield (Tang et al., 2018).
The advent of next-generation sequencing technologies has paved the way to cost-effective and highly efficient methods to identify miRNA–mRNA regulatory networks related to complex traits in crops. For instance, embryo development in peanut (Arachis hypogaea) is highly sensitive to calcium deficiency in the soil, but the cause for the resulting embryo abortion was unknown. However, an integrated miRNA and mRNA profiling (RNA-seq) analysis revealed that a number of miRNA-mediated regulatory networks affecting seed/embryo development, cell division, cell proliferation, and plant hormone signaling all participated in peanut embryo abortion under calcium deficiency (Chen Z. Y. et al., 2019). Joint RNA-Seq and miRNA profiling analyses further established that miRNAs involved in nitrogen-related pathways regulated the thickness of a pod canopy in B. napus (Chen Z. Y. et al., 2019). To date, limited information is available pertaining to miRNA-mediated regulatory networks related to HI.
To better understand the regulatory networks underlying the establishment of HI in B. napus during the seed-filling phase, we determined the mRNA and miRNA transcriptome landscape in SP and seeds during the seed-filling stage in plants grown at two locations, Chongqing and Yunnan. We identified potential gene clusters involved in the regulation of HI, with predicted roles in transporter activity and responses to environmental signals. Furthermore, we discovered several miRNA-mediated regulatory networks in SP and seeds. These results contribute to uncovering the complex regulatory networks behind HI and offering potential solutions for its improvement via genetic engineering or crop breeding.
Materials and Methods
Plant Materials and Field Trials
The B. napus accessions YC24 (SWU47), YC52 (Zhongshuang11), and YC46 (Ningyou12) with stability HI for two consecutive years (2012–13 and 2013–14) were selected from 520 accessions in previous study (Lu et al., 2016) and grown in a randomized block design with three replications at Chongqing Beibei (CQ, 29°45′ N, 106°22′ E, 238.6-m altitude) and Yunnan Lincang (YN, 23°43′ N, 100°02′ E, 1819.5-m altitude) during the 2015–2016 growing season. For simplicity, we denoted the planting location first (CQ or YN), followed by the accession number (24, 52, or 46) to distinguish experimental groups. Planting conditions were as previously described (Lu et al., 2016). Plant materials used for RNA-seq and small RNA-seq of growing status are shown in Supplementary Figure 1. We collected HI-related phenotypic traits, such as HI, biomass yield per plant (BY), seed yield per plant (SY), stem dry weight (ST), canopy biomass yield (CBY), according to the calculation method used previously (Lu et al., 2017), and the HI, BY, and SY of CQ24 and CQ46 at 2016 had been shown in our previous study (Zhang et al., 2020).
RNA Isolation and Library Preparation
For each accession grown at CQ and YN, we harvested seeds (ZS) and SP from the main inflorescence at 30 days after flowering (denoted as 30ZS and 30SP, respectively); there are also leaves (Le, main leaves at 30 days after flowering) and stems (St, main steam at 30 days after flowering) of the same period. For each sample, we collected two biological replicates, each harvested from five independent plants.
Total RNA of each sample was extracted using a CTAB method (Lu et al., 2008). RNA degradation and DNA contamination were monitored by gel electrophoresis on 1% agarose gels. RNA purity was confirmed on a NanoPhotometer spectrophotometer (IMPLEN, CA, United States), and RNA concentration was measured using the Qubit RNA Assay Kit with a Qubit 2.0 Fluorometer (Life Technologies, CA, United States). RNA integrity was assessed using the RNA 6000 Nano Assay Kit on an Agilent Bioanalyzer 2100 system (Technologies, CA, United States). After quality-control, we sent 48 RNA samples for mRNA sequencing and selected 12 samples (including 30SP, 30ZS) for small RNA sequencing to Novogene Corporation (Beijing, China) for library construction and sequencing, as previously described (Chen H. et al., 2019).
Identification of Differentially Expressed Genes
We evaluated the quality of RNA-Seq data with the Trimmomatic software (v.0.36) (Bolger et al., 2014). We removed adapter sequences and low-quality reads from the raw data. All clean reads were then mapped to the B. napus reference genomev.4.11 using the STAR program (v.2.5.3) (Dobin et al., 2013) with default parameters. We then mapped the aligned reads to RNA features using the feature Counts function of Subread (v.1.6.0) (Liao et al., 2014). We identified [differentially expressed gene (DEGs)] with the gene counts generated above with the edgeR package (Robinson et al., 2010) to import, organize, filter, and normalize the data with a false discovery rate (FDR) of 1%. We quantified relative gene expression as fragments per kilobase of exon model per million mapped reads (FPKM) by Cuffdiff within Cufflinks (Trapnell et al., 2012). After calculating their expression fold-change, we identified genes with FDR-adjusted q-value ≤ 0.05 and absolute log2 (fold-change) ≥ 1 as DEGs.
Identification of Differentially Expressed Conserved and Novel MicroRNAs
Clean data were obtained from small RNA sequencing libraries by removing adapter sequences from all reads, as well as reads containing over 10% Ns and low-quality (Q20 < 85%) reads from the raw data. We extracted potential small RNAs 18–30 nt in length and mapped them to the B. napus reference genome using Bowtie2 (v.2.2.9) (Langmead et al., 2009). Mapped reads were further mapped to the released B. napus miRNAs in miRBase22 (v.22.1)2 using basic local alignment sequence tool for nucleotides (BLASTN) with an E-value cutoff of ≤ 1e-5, which identified known miRNAs. We used the Rfam database (v.14.4)3 to remove ribosomal RNAs (rRNAs), transfer RNAs (tRNAs), small nuclear RNAs (snRNAs), and small nucleolar RNAs (snoRNAs). In order to improve the accuracy and efficiency of comparison, we retained only the longest sequence from the remaining reads and merged them with the sequences in miRBase22 to become the final B. napus miRNA dataset. After removing the reads already classified as known miRNAs, we employed miRDeep-P (v.1.1.4) (Yang and Li, 2011) to predict potential novel miRNAs, allowing for 2 bp 3′ overhangs. We named these potential novel miRNAs based on the chromosome they map to and their starting nucleotide position on that chromosome.
To identify differentially expressed miRNAs (DEMs), we normalized miRNA expression across all samples using Salmon (v.0.11.3) (Patro et al., 2017) to obtain the expression of transcript per million (TPM) based on the normalization formula: normalized expression = (actual miRNA counts/total number of mapped reads) × 1,000,000. We calculated the associated fold-change, p-values, and q-values with in-house Perl scripts. miRNAs with FDR-adjusted q-value ≤ 0.05 and absolute log2 (fold-change) ≥ 1 were determined to be DEMs.
MicroRNA Target Prediction
We employed the psRNATarget server (v. 2017)4 to predict the target genes for all known and novel miRNAs that were differentially expressed between comparable groups, with default parameters.
Gene Function Clustering Analysis
Functional annotation was performed by using BLASTX to compare B. napus and Arabidopsis proteins. We used the online tool agriGO (gene ontology analysis toolkit and database for agricultural community, v.2.0)5 (Tian et al., 2017) for GO enrichment analysis. We determined GO classifications by submitting the gene sequences to the BLAST4ID tool of agriGO to obtain the corresponding B. napus locus ID, and then running GO term enrichment analysis.
Weighted Gene Co-expression Network Analysis
Weighted gene co-expression network analysis (WGCNA) package in R was designed in 2008 for helping users create weighted correlation network modules and identify key genes associated with traits in interesting modules (Langfelder and Horvath, 2008). FPKM values of genes were log2 (FPKM + 1) transformed for further calculation of correlation coefficient, determination of gene modules, construction of co-expression network, and correlation analysis of modules and phenotypic traits. In the process of analysis, the soft thresholding power was determined using the pickSoftTreshold function based on the scale-free topology model fit (R2) > 0.9; the automatic blockwiseModules network construction approach was applied to obtain the highly correlated modules, with the following parameters: power = 7; maxBlockSize = 30,000; TOM-type = unsigned; miniModuleSize = 30; reassignThreshold = 0, mergeCutHeight = 0.25. The co-expression networks were displayed using Cytoscape (v.3.5.1) (Maere et al., 2005).
Validation of Transcript Levels by Quantitative Reverse Transcription-PCR
The same RNA samples prepared for mRNA and miRNA sequencing libraries were used for cDNA synthesis and quantitative reverse transcription-PCR (qRT-PCR) detection. qRT-PCR was performed as described previously (Qu et al., 2015). We used gene-specific primers obtained from qPrimerDB6 (Lu et al., 2018). Relative transcript levels were normalized to the B. napus housekeeping genes UBIQUITIN-CONJUGATING ENZYME 21 (Bna.UBC21) and Bna.ACTIN7 (Qu et al., 2016).
For the validation of miRNAs, we added a poly (A) tail and performed reverse transcription from 2-μg RNA in 20-μl reaction volume using the miRcute miRNA First-Strand cDNA Synthesis Kit (Tiangen, Beijing, China). We then diluted the cDNAs 10-fold for RT-qPCR analysis, using the miRcute miRNA qPCR Detection Kit (SYBR) (Tiangen, Beijing, China). We carried out PCR reaction solutions containing 1-μl (∼10 ng) diluted cDNA, 10-μl 2- × -miRcute miRNA premix, a 0.4-μl forward primer (10 μM), a 0.4-μl reverse primer (10 μM), and 8.2-μl ddH2O on a BIO-RAD CFX96 Real-Time System (BIORAD, United States). Cycle conditions were: 95°C for 15 min, followed by five cycles of 94°C for 20 s, 65°C for 30 s, and 72°C for 34 s, and then 40 cycles of 94°C for 20 s and 60°C for 30 s. All reactions were performed in triplicate, with the B. napus U6 snRNA as internal control. Relative gene expression was calculated using the 2–ΔΔCt method (Livak and Schmittgen, 2001). All above-mentioned qRT-PCR assays were performed according to MIQE guidelines (Bustin et al., 2010). Three independent biological replicates, each with three technical replicates, were implemented for each sample.
Results
Physiological Characteristics of Accessions With Different Harvest Index
In the context of this study, we planted the B. napus accessions YC24, YC52, and YC46 at the Beibei, Chongqing (CQ), and Lincang, Yunnan (YN) locations during the 2015–2016 growing season. All phenotypic values are provided in Table 1; in order to distinguish experimental groups, they were named as planting location first (CQ or YN) plus the accession number (24, 52, or 46). We collected HI-related phenotypic traits, such as HI, biomass yield per plant (BY), seed yield per plant (SY), stem dry weight (ST), canopy biomass yield (CBY). We noticed that the HI for all accessions significantly increased at the YN location when compared to that measured at CQ, especially accessions YC24 and YC52. In addition, YC24 displayed a significantly higher HI than the YC52 and YC46 accessions at both locations (Table 1); therefore, YC24 was regarded as a high HI accession, and YC52 and YC46 as low HI accessions. We used these three accessions to elucidate the molecular mechanism of HI.
Table 1. Phenotypic values of the three varieties measured in Chongqing and Yunnan.
mRNA Sequencing Data Analysis Uncovers Harvest Index-Related Differentially Expressed Genes
We collected 48 samples for RNA extraction and subsequent deep sequencing of the transcriptome (RNA-Seq) analysis: ZS, SP, Le, and St in biological duplicates from all three accessions at the two geographic locations. After removal of low-quality reads from the raw data, we retained 40.29∼52.68 million clean reads; 83.97 to 95.23% of which were mapped to the B. napus reference genome (except SP in CQ24 Sample 1, which shows just 61.67%). Sample correlation analysis used FPKM to emphasize the high degree of correlation between biological replicates (Supplementary Table 1). The raw sequencing data were deposited in the BIG Data Center under BioProject accession No. PRJNA602979.
We identified DEGs across SP samples; we detected 2,303 DEGs (1,212 upregulated and 1,091 downregulated) in CQ24 vs. CQ52; 4,680 DEGs (2,472 upregulated and 2,208 downregulated) in CQ24 vs. CQ46; 1,103 DEGs (641 upregulated, 462 downregulated) in YN24 vs. YN52; 1,904 DEGs (1,092 upregulated, 812 downregulated) in YN24 vs. YN46 (Supplementary Table 2). To distinguish HI-related DEGs between high- and low-HI lines, we generated a four-way Venn diagram showing the extent of overlap between the four comparisons listed above. This analysis highlighted 756 DEGs (393 upregulated and 363 downregulated) in SP between high- and low-HI accessions grown at the CQ location, and 343 DEGs (203 upregulated and 140 downregulated) between high- and low-HI accessions grown at the YN location. A subset of 146 DEGs (83 upregulated and 63 downregulated) was identified in SP samples collected at both locations and across all accessions (Figures 1A,B).
Figure 1. Venn diagrams and the number of DEGs of the DEGs from high- and low-HI accessions. (A) A Venn diagram of DEGs between high- and low-HI accessions in ZP at CQ and YN. (B) The number of DEGs between high- and low-HI accessions in ZP. (C) A Venn diagram of DEGs between high- and low-HI accessions in Le at CQ and YN. (D) The number of DEGs between high- and low- HI accessions in Le. (E) A Venn diagram of DEGs between high- and low-HI accessions in ZS at CQ and YN. (F) The number of DEGs between high- and low-HI accessions in ZS. (G) A Venn diagram of DEGs between high- and low-HI accessions in ST at CQ and YN. (H) The number of DEGs between high- and low-HI accessions in ZP in ST.
We repeated the same analysis with the RNA-Seq data obtained from the ZS samples; we identified 2,192 DEGs (1,481 upregulated, 711 downregulated) in CQ24 vs. CQ52; 2,776 DEGs (1,430 upregulated and 1,346 downregulated) in CQ24 vs. CQ46; 850 DEGs (562 upregulated, 288 downregulated) in YN24 vs. YN52; 1,963 DEGs (1,281 upregulated, 682 downregulated) in YN24 vs. YN46 (Supplementary Table 2). The corresponding four-way Venn diagram showed that 438 DEGs (236 upregulated and 202 downregulated) were shared by high- and low-HI accessions grown at the CQ location, while 205 DEGs (149 upregulated and 56 downregulated) were common to high- and low-HI accessions grown at the YN location. Furthermore, 118 genes (83 upregulated and 35 downregulated) were differentially expressed in seed samples collected at both locations and across all accessions (Figures 1C,D).
Leaves are the same “source” organs as the SP; through RNA-seq data, we identified 2,500 DEGs (1,227 upregulated, 1,273 downregulated) in CQ24 vs. CQ52; 6,047 DEGs (3,019 upregulated and 3,026 downregulated) in CQ24 vs. CQ46; 4,660 DEGs (2,185 upregulated, 2,475 downregulated) in YN24 vs. YN52; 6,250 DEGs (2,833 upregulated, 3,417 downregulated) in YN24 vs. YN46 (Supplementary Table 2). Our 4-way Venn diagram showed that 783 DEGs (371 upregulated and 412 downregulated) were shared by high- and low-HI accessions grown at the CQ location, while 1,369 DEGs (588 upregulated and 781 downregulated) were common to high- and low-HI accessions grown at the YN location. And 552 (242 upregulated and 310 downregulated) genes were differentially expressed in Le samples collected at both locations and across all accessions (Figures 1E,F).
Stem as the “flow” organ, which is the limitation for the accumulation of assimilates in B. napus seeds (Shen et al., 2010; Fu and Zhou, 2013; Luo et al., 2015), we also analysis stem RNA-seq data; there are 1,415 DEGs (722 upregulated, 693 downregulated) in CQ24 vs. CQ52; 5,313 DEGs (2,652 upregulated and 2,761 downregulated) in CQ24 vs. CQ46; 1,143 DEGs (623 upregulated, 520 downregulated) in YN24 vs. YN52; 4,866 DEGs (2,559 upregulated, 2,307 downregulated) in YN24 vs. YN46 (Supplementary Table 2). And 418 DEGs (214 upregulated, 204 downregulated) in St between high- and low-HI accessions grown at the CQ location; 311 DEGs (199 upregulated, 112 downregulated) in St between high- and low-HI accessions grown at the YN location. About 227 (132 upregulated, 95 downregulated) genes were differentially expressed in St samples collected from the two locations (Figures 1G,H).
Overall, the transcriptome of the accessions YC24 and YC46 differed by more DEGs than when YC24 was compared to the YC52 accession, both in SP, Le, ZS, and St. In addition, SP and Le samples as “source” organs were characterized by more DEGs than St and ZS samples between high- and low-HI accessions at both locations, indicating that the regulation of HI in “source” organs might be more complicated. Finally, samples harvested at the CQ location exhibited more DEGs than at the YN location in SP, ZS, and St, except Le, suggesting that the seed-filling process might be more complex at the CQ location.
In addition, in order to compare the gene expression differences of the same material in different regions, we compared the DEGs of the same material grown in YN and CQ. The results indicate that these DEGs among the three materials in four different plant tissues are quite different in different environments and the expression levels also varied greatly (Supplementary Table 2), especially in the Le of material YC24; there were 4,599 upregulated DEGs and 4,659 downregulated DEGs (Supplementary Figure 2). At the same time, the DEGs in the same tissue of the three materials also have differences; only 16, 36, and 34 same DEGs were found in the St, Le, and SP, and the absence of the same DEGs was found in the ZS (Supplementary Figure 2). Moreover, combining with Table 1, we found that HI in YN was significantly higher than that in CQ, indicating that HI was easily affected by environmental conditions; this is worthy of further study.
Functional Annotation and Classification of Differentially Expressed Genes
To understand the functions encoded by the DEGs identified between high- and low-HI accessions, we performed a GO enrichment. Studies showed long ago that the photosynthates of the SP are the main sources of seed yield, contributing about 2/3 of the total dry matter of the seed yield, whereas the beak of the seeds is about 8% (Leng et al., 1992), which can seriously affect HI, so we focused on SP and ZS first.
Differentially expressed genes in SP, harvested at the CQ location between high- and low-HI accessions (SP_CQ_High vs. Low), were associated with 93 significantly enriched GO terms (q-value ≤ 0.05), while DEGs for the equivalent samples collected at the YN location (SP_YN_High vs. Low) showed a significant enrichment for 341 GO terms; in addition, we identified 54 significantly enriched GO terms (q-value ≤ 0.05) for SP among DEGs shared by both CQ and YN locations (Supplementary Figures 3A,B). We observed a significant enrichment in several cellular component GO terms related to “chloroplasts” (GO: 0005737, GO: 0009507, GO: 0009536, GO: 0044434, GO: 0044435), and we noticed that GO terms of SP in addition to response to biological stress and a large number of abiotic stresses (GO:0009628, response to abiotic stimulus; GO:0009408, response to heat; GO:0009607, response to biotic stimulus; etc.); there are some interesting pathways enriched in “circadian rhythm” (GO:0007623), “histone H3-K36 demethylation” (GO:0070544), and “response to Karrikin” (GO:0080167), which suggest that circadian rhythm, histone modification affected the plant HI; also, as Karrikin, a new signaling molecule participates in regulating HI (Supplementary Figure 3C).
In ZS samples, a similar analysis revealed 107 significantly enriched GO terms among DEGs between high- and low-HI accessions harvested at the CQ location (ZS_CQ_High vs. Low), and 59 significantly enriched GO terms among DEGs between high- and low-HI accessions harvested at the YN location (ZS_YN_High vs. Low) (Figures 2A,B). We identified 18 common enriched GO terms among DEGs for ZS. Biological processes-related GO terms, such as the “negative regulation of RNA metabolic process” (GO: 1902679), “negative regulation of nucleic acid-templated transcription” (GO: 1903507), “negative regulation of the RNA metabolic process” (GO: 0051253), and the “organonitrogen compound metabolic process” (GO: 1901564), were predominantly enriched in ZS samples harvested at both CQ and YN locations (Figure 2C). In the YN location, we noticed “circadian rhythm” (GO: 0007623) and “alternative mRNA splicing via spliceosome” (GO: 0000380) were identified but not identified in the CQ location. We suspected that, in ZS samples, circadian rhythm and alternative mRNA splicing may participate in the HI regulation but were also influenced by geographical differences and environmental factors.
Figure 2. GO functional classification of the ZS DEGs from high- and low-HI accessions at CQ and YN. (A) ZS_CQ_High vs. Low. (B) ZS_YN_High vs. Low. (C) Common GO terms shared by CQ and YN varieties. The color bars under the figures represent -log10 (p.adjust values).
GO enrichment analysis of St and Le DEGs in the CQ location obviously enriched in “circadian rhythm” (GO: 0007623) as the results in SP and ZS; we speculated that circadian rhythm plays an important role in regulating HI. We also found St and Le DEGs are mainly enriched in hormone-related terms, such as the “ethylene-mediated signaling pathway” (GO: 0009873), “jasmonic acid-mediated signaling pathway” (GO: 0009867) in St at the CQ location, “response to gibberellin stimulus” (GO: 0009739) in St at the YN location. In Le samples, “response to gibberellin stimulus” (GO: 0009739) and the auxin metabolic/biosynthetic process (GO: 0009850 and GO: 0009851) also enriched, suggesting that, although the leaf no longer provides most of the photosynthetic energy during silique ripening (SP provides more), it works with the St to regulate plant growth hormonally (Supplementary Figure 4).
Construction of Co-expression Networks
Through WGCNA, we constructed co-expression networks with all DEGs and 13 phenotypic data (Supplementary Table 3). This analysis yielded 72 gene modules, each represented by different colors in the output (Figure 3). The smallest module is the ME light coral, which included only 55 genes, and the largest module is turquoise, which included 10,008 genes. We focused on the ME light green module that is highly correlated with the HI (r2 = 0.61).
Figure 3. An overview of module and traits corresponding through WGCNA.
GO enrichment analysis showed the ME light green module was mainly enriched in the “ubiquitin-dependent protein catabolic process” (GO: 0006511), “entrainment of the circadian clock” (GO: 0009649), “positive regulation of meiosis” (GO: 0045836), and “production of miRNAs involved in gene silencing by miRNA” (GO: 0035196). This indicated that ubiquitin modification, circadian rhythm, and gene silencing are all strong connection with HI; the results are the same with our previous DEGs analyzed. We used the maximal clique centrality method in the cytoHubba plugin of Cytoscape (v. 3.5.1) to identify hub genes in the ME light green module of interest, and we identified BnaA07g02330D (Bna.A07STY46), which is serine/threonine kinase that phosphorylates transit peptides of chloroplast and mitochondria-targeted pre-proteins and is involved in chloroplast differentiation in Arabidopsis (Giorgia et al., 2011). We speculated that Bna.A07STY46 may participate in retarded establishment of a photosynthetic capacity to influence HI.
Identification of Differentially Expressed MicroRNAs Based on MicroRNA Sequencing
In addition to 48 mRNA libraries, we also sequenced 12 miRNA libraries of samples SP and ZS. We obtained 163,448,662 reads and retained 159,966,665 clean reads after removing low-quality reads and adapters (Supplementary Table 4). We then selected clean reads with a length of 18–30 nt for further analysis. As expected, reads with a length of 21, 22, and 24 nt were more abundant out of all reads. In addition, we observed a higher fraction of 21-nt reads in silique pericarp samples when compared to that in seed samples, while 24-nt reads showed the opposite pattern (Supplementary Figure 5). We identified 903 miRNAs from the 12 miRNA libraries, including 46 known miRNAs and 857 novel miRNAs (Supplementary Table 5).
Based on the selection criteria of an absolute log2 FC > 1 and a q-value < 0.05, we detected DEMs between the same comparison groups as for RNA-Seq analysis. When we compared SP_CQ_High and Low, we identified 11 known DEMs (1 upregulated and 10 downregulated) and 86 novel DEMs (60 upregulated and 26 downregulated), while the comparison of SP_YN_High with Low revealed 15 known DEMs (0 upregulated and 15 downregulated) and 80 novel DEMs (38 upregulated and 42 downregulated). A comparison between ZS_CQ_High and Low yielded 10 known DEMs (8 upregulated and 2 downregulated) and 70 novel DEMs (38 upregulated and 32 downregulated), while we identified 11 known DEMs (11 upregulated and 0 downregulated) and 65 novel DEMs (30 upregulated and 35 downregulated) from a comparison between ZS_YN_High and Low. In addition, SP samples harvested at the CQ and YN locations shared 8 known DEMs and 21 novel DEMs, and ZS samples collected at the two locations saw an overlap consisting of 8 known DEMs and 16 novel DEMs (Supplementary Figure 6).
Integration of Differentially Expressed Genes and Differentially Expressed MicroRNAs
To elucidate the regulatory role of DEMs between high- and low-HI accessions, we first identified the potential target genes of each DEM before combining the expression profiles of DEMs and their target genes for further analysis. We thus obtained miRNA–mRNA interaction pairs (pairs with either negatively or positively correlated expression patterns) through a comparison of high- and low-HI accessions; 130 pairs in SP were harvested at CQ, 68 pairs in SP were collected at YN, 69 pairs in ZS from the CQ location, and 23 pairs in ZS from the YN location. Overall, almost half of the miRNA–mRNA interactions pairs showed a negative correlation, as might be expected from a true mRNA–miRNA pair involving transcript cleavage. We also noticed that several miRNAs had multiple potential target mRNAs, and multiple miRNAs that targeted a single mRNA (Supplementary Figure 7 and Supplementary Table 6). For instance, the upregulated miRNA bna-miRC03_52 controlled the downregulated genes BnaC07g39100D, BnaC08g32460D, and BnaA06g21420D, while the upregulated miRNA bna-miRA02_2282 controlled the downregulated genes BnaA06g21420D and BnaA06g39650D in comparisons between SP_CQ_High and Low (Supplementary Figure 7). These observations indicate that the regulation of miRNA–mRNA pairs is complex.
We performed a GO functional annotation analysis to characterize these differentially expressed target genes. Our results showed enrichment in GO terms related to circadian rhythms, transporter activity response to stress, response to abiotic stimulus, and response to stimulus across all comparisons, with the exception of ZS_CQ_High vs. Low (Supplementary Figure 8 and Supplementary Table 7). We hand-selected miRNA–mRNA pairs associated with the GO terms mentioned above, such as the bna-miRC08_5718–BnaA03g03740D/BnaC03g05240D pair and the bna-miRC01_19092–BnaC09g43920D pair (Figure 4A and Supplementary Figure 8). We did not discover significantly enriched GO terms among miRNA–mRNA pairs from seed samples collected at the CQ and YN locations between high- and low-HI accessions. However, the bna-miR396–BnaA06g21030D and bna-miR396–BnaA01g33410D pairs were shared between the two locations (Figure 4B and Supplementary Figure 8).
Figure 4. Heatmap representation of interesting miRNA-mRNA pairs identified in silique pericarps and seeds. (A) Heatmap representation of miRNA-mRNA pairs-related circadian rhythm, response to stress, response to abiotic stimulus, and response to stimulus in silique pericarps. (B) Heatmap representation of bna-miR396-modulated miRNA-mRNA pairs in seeds. The color bars under the figures represent log2 (FPKM/TPM).
Validation of mRNA and MicroRNA Expression by Quantitative Reverse Transcription-PCR
To validate the quantification of the mRNA and miRNA sequencing data presented here, we analyzed the relative transcript levels of nine randomly selected DEGs (BnaA01g26430D, BnaA08g25340D, BnaA08g26300D, BnaC02g04730D, BnaC02g12960D, BnaC04g50590D, BnaC06g05910D, BnaC09g48250D, and BnaCnng27780D) and 7 DEMs (bna-miR156a, bna-miR164a, bna-miR167a, bna-miR396, bna-miRA01_10807, bna-miRC03_27293, and bna-miRC04_29081) by qRT-PCR on the same RNA used for library construction (Supplementary Table 8 and Figure 5). We observed a high degree of positive correlation between the relative transcript levels of DEGs and DEMs determined by qRT-PCR and their relative expression measured from high-throughput sequencing data. The high-throughput sequencing data used in this study were, therefore, accurate and reliable.
Figure 5. RT-qPCR validation of DEGs and DEMs. Blue bars represent the relative expression level of RT-qPCR, and red points represent FPKM/TPM of sequencing; r represents Pearson correlation coefficient, *: correlation is significant at the 0.05 level, **: correlation is significant at the 0.01 level.
Discussion
Harvest index is a complex agronomic trait of great economic value that depends on interactions between a plant genotype and the environment (Amanullah and Shah, 2010; Unkovich et al., 2010; Gajića et al., 2018; Chao et al., 2019). We determined that the YC24 accession showed a significantly higher HI than the YC52 accession, while the YC52 accession had a significantly higher HI than the YC46 accession at both the CQ and YN locations. These results indicated that HI differences between the YC24, YC52, and YC46 accessions were robust, validating the use of these three accessions to dissect the regulatory mechanism behind HI. YN is a high-yield crop production environment; not surprisingly, HI for the YC24, YC52, and YC46 accessions grown at this location was higher than that from the CQ location, indicating that environmental conditions can have a strong influence on HI. Thus, we carried out transcriptome sequencing to compare the gene expression profile in accessions grown in a standard HI (CQ) and high-HI (YN) environment during the seed-filling stage.
RNA-Seq and Expression Profiles of High- and Low-Harvest Index Materials
Based on gene functional annotation data, we identified significantly enriched GO terms associated with DEGs between high- and low-HI accessions. In B. napus, SP are photosynthetically active, and, as source tissues, they thus provide nutrients to seeds, a sink organ (Diepenbrock, 2000; Bennett et al., 2011). Furthermore, SP mediate maternal control during seed filling (Li et al., 2019). We also noticed that a large fraction of enriched GO terms were specific to the CQ location, as we did not detect them in our analysis of silique pericarp or seed samples collected at the YN location, underscoring the marked environmental sensitivity of HI regulation.
As HI represents the ratio between reproductive organs and vegetative biomass produced, regulating the timing of the phase transition from vegetative to reproductive growth will be beneficial to increasing HI (Hütsch and Schubert, 2017). Based on GO enrichment analysis, we noticed that the circadian rhythm may be associated with HI. In Arabidopsis, the core circadian rhythm comprises the proteins CIRCADIAN CLOCK ASSOCIATED 1 (CCA1), LATE ELONGATED HYPOCOTYL (LHY), TIMING OF CAB EXPRESSION 1 (TOC1), and PSEUDO RESPONSE REGULATORS (PRRs) (Mcclung and Gutiérrez, 2010). In addition, the CONSTANS gene family (CO) controls photoperiodic flowering time together with PRRs (Nakamichi et al., 2007). Transcriptome analysis revealed that genes associated with circadian rhythm were potentially involved in potato (Solanum tuberosum) tuber formation, suggesting that the circadian rhythm may participate in photoassimilate distribution (Shan et al., 2013). In rice, CONSTANS-like (OsCOL9) not only modulated photoperiodic flowering but influenced grain number of the main panicle (Liu et al., 2016). Similarly, TOC1 modulated chickpea (Cicer arietinum) seed yield per plant (Basu et al., 2019), while overexpression of Arabidopsis PRR5 in rice delayed flowering but also significantly increased biomass (Nakamichi et al., 2020), likely due to higher expression of OsPRR37; the overexpression line GCA1OX-5 showed better general combining ability of rice (Liu et al., 2015). In wheat, the Photoperiod-1 (Ppd-1) mutant inactivates a PRR and affects paired spikelet formation (Boden et al., 2015). Moreover, our analysis highlighted additional circadian rhythm-related genes [GLYCINE RICH PROTEIN7 (GRP7), EARLY FLOWERING4 (ELF4), REVEILLE 1 (RVE1)] as being differentially expressed between high- and low-HI lines, strongly suggesting that the circadian rhythm and rhythm-controlled gene regulation might be harnessed to modulate HI in B. napus in the future.
Under the high-HI environment of the YN location, we detected a number of DEGs expressed in SP that were related to transporter activity, such as PLASMA MEMBRANE INTRINSIC PROTEIN (PIP), BILE ACID TRANSPORTER 5 (BAT5), ABC2 HOMOLOG 13 (ATH13), CHLORIDE CHANNEL A (CLC-A), MAJOR FACILITATOR SUPERFAMILY PROTEIN, NITRATE TRANSPORTER 1.7 (NRT1.7), POLYOL/MONOSACCHARIDE TRANSPORTER 5 (PMT5), and VACUOLAR GLUCOSE TRANSPORTER 1 (VGT1). PIP are aquaporins that localize to the plasma membrane and facilitate the flux of water and solutes across the plasma membrane (Kourghi et al., 2017; Wang H. et al., 2020). Multiple studies have shown that PIP genes can affect water balance and solute transport in plant cells (Sommer et al., 2008; Byrt et al., 2017; Macho-Rivero et al., 2018). The BAT5 transporter translocates glucosinolates, which are derived from methionine and sugars, across the chloroplast membranes (Sawada et al., 2009). ATH13 affects the lipid composition of chloroplast membranes and regulates iron distribution within chloroplasts (Manara et al., 2014, 2015). CLC-A regulates the accumulation of nitrate within vacuoles and plays an important role in modifying cytosolic conditions (De Angeli et al., 2006; Manara et al., 2015; Demes et al., 2020). Additional regulators of nitrate balance and transport include members of the NITRATE TRANSPORTER1/PEPTIDE TRANSPORTER (NPF) family, such as AtNPF3.1 (At1g68570) (Yang et al., 2017) and AtNPF6.2 (At2g26690) (Tong et al., 2016); AtNRT1.7/NPF2.13 (At1g69870) also plays important roles in source-to-sink remobilization of nitrate in Arabidopsis (Fan et al., 2009; Liu et al., 2017). PMT5 can transport a wide range of linear polyols (three to six carbon backbones), cyclic polyols (myo-inositol), pyranose, furanose, hexoses, and pentoses across the plasma membrane (Klepek et al., 2005, 2010). Collectively, the upregulation of genes-encoding proteins with transporter activity in SP of high-HI accessions points to their possible involvement in facilitating water and solute fluxes (e.g., sugar, nitrate nutrients, and cations) from the mother plant to the developing seed, and thus in increasing HI. The only downregulated transporter identified in this study was VGT1. However, VGT1 localizes to the vacuolar membrane and mediates the transport of glucose from the cytoplasm to the vacuole (Aluri and Büttner, 2007; Büttner, 2007); a reduction in VGT1 expression may, therefore, promote sugar translocation to sink tissues by a modulating glucose flux to the vacuole.
In seeds, GO terms, such as GO: 0045892, GO: 1902679, GO: 1903507, GO: 0051253, GO: 0045934, GO: 0010629, GO: 0010605, and GO: 1901564, were significantly enriched in samples harvested at both the CQ and YN locations. For the latter GO term (the organonitrogen compound metabolic process, GO: 1901564), the B. napus gene BnaC08g45660D homologous to At1g01090 (PYRUVATE DEHYDROGENASE E1 ALPHA, PDH-E1 ALPHA) was upregulated in high-HI accessions. As PDH-E1 ALPHA has been shown to regulate acyl lipid metabolism (LeClere et al., 2004; Mentzen et al., 2008), we speculated that the upregulation of BnaC08g45660D contributes to the accumulation of storage lipids in seeds. Gene functional annotation of the GO terms GO: 0045892, GO: 1902679, GO: 1903507, GO: 0051253, GO: 0045934, GO: 0010629, and GO: 0010605 determined that the expression pattern of BnaA05g01050D (homologous to At2g46830, CCA1) showed the same trend at both the CQ and YN locations. Studies have shown that CCA1 affects not only the circadian rhythm but also the accumulation of storage lipids (Kim et al., 2019). Thus, we concluded that genes that promote the accumulation of storage materials such as lipids during the seed-filling stage likely contribute to the improvement of HI.
MicroRNA-Mediated Regulatory Networks Related to Harvest Index
MicroRNAs play versatile roles in plant growth and development via miRNA–mRNA interaction networks. In our study, we identified abundant miRNA–mRNA interaction pairs from the comparison of silique pericarp and seed samples from high- and low-HI accessions at CQ and YN locations. Based on the functional analysis of target DEGs, we noticed enrichments for genes in silique pericarp samples harvested at both CQ and YN locations related to responses to stimulus. Thus, we selected genes related to circadian rhythm, response to stress, response to abiotic stimulus, and response to stimulus that also exhibited anti-correlations with their miRNAs for further analysis. MiR164 has been shown to control axillary meristem and floral organogenesis formation in Arabidopsis (Raman et al., 2008; Huang et al., 2012); in our study, we observed a negative correlation between bna-miR164a and BnaA08g26300D. BnaA08g26300D shows homology to the Arabidopsis gene At1g09560 (GERMIN-LIKE PROTEIN 5, GLP5). GLP may control resource allocation between primary and lateral roots by phloem-mediated transport in Arabidopsis (Ham et al., 2012), suggesting that the regulation of GLP5 transcript levels by bna-miR164a may adjust the balance of resources between SP and seeds though the phloem. At5g06530 (ABC TRANSPORTER GENE 22, AtABCG22) also contributes to water transpiration and drought tolerance in Arabidopsis (Kuromori et al., 2011). Hence, we propose that its B. napus homolog BnaC02g01960D, whose transcript levels are regulated by bna-miRA03_11110, might be important for the proper water balance of SP during the seed-filling phase.
The circadian rhythm gene TOC1 contributes to energy metabolism by influencing the phase of the circadian rhythm under environmental fluctuations (Legnaioli et al., 2009; Fung-Uceda et al., 2018) and Arabidopsis; toc1 mutants showed a modified pattern of starch mobilization under light–dark cycles (Flis et al., 2019). We hypothesized that the bna-miRA07_12065–BnaC09g05250D (TOC1) pair might similarly affect the phase of the circadian rhythm to modulate resources allocation in SP. BnaC09g43920D is homologous to At5g13170 (SENESCENCE-ASSOCIATED GENE 29, SAG29, also named SWEET15). In this study, this gene is upregulated in high-HI SP at the YN location via the downregulation of bna-miRC01_19092. SWEET15, together with the transporters SWEET11 and SWEET12, mediates sucrose efflux both intracellularly and intercellularly during seed filling in Arabidopsis (Chen et al., 2015; Eom et al., 2015) and may also regulate senescence under environmental stress (Seo et al., 2011). Moreover, we noted the upregulation of several heat stress-related genes, such as HEAT SHOCK PROTEIN (HSP) Hsp17.6CII (BnaA03g03740D/BnaC03g05240D), HSP26.5 (BnaC06g05390D), CPHSC70-1 (BnaC01g16200D), and HEAT SHOCK FACTOR 4 (HSF4) (BnaCnng56320D), this upregulation being mediated by the downregulated miRNAs bna-miRC08_5718, bna-miRC06_3193, bna-miRC02_6475, and bna-miRC01_11705 in high-HI lines, respectively. Overall, we identified many miRNA–mRNA pairs related to environmental stress in the current study, indicating that increasing plant adaptability to the environment may adjust resources distribution and improve HI.
Similar to our observation with DEGs, numerous miRNA–mRNA pairs identified here differed between seeds and SP, indicating distinct regulatory mechanisms during the seed-filling stage. Among miRNA–mRNA pairs, bna-miR396 and its putative targets BnaA06g21030D and BnaA01g33410D were shared by seed samples harvested from both the CQ and YN locations. MiR396 exerts a strong influence in plant development by regulating complex traits; for example, miR396 regulates the transcript levels of the GRF group to modulate cell proliferation and elongation (Ercoli et al., 2016), grain size and yield (Che et al., 2015; Li et al., 2016; Chen X. L. et al., 2019), and somatic embryogenesis (Szczygieł-Sommer and Gaj, 2019). Moreover, a loss of function in miR396ef resulted in higher grain size and altered plant architecture in rice (Miao et al., 2020), while over-expression of oa-miR396c in rice reduced salt and alkali stress tolerance (Gao et al., 2010). In tomato (Solanum lycopersicum), over-expressing a short tandem repeat mimic for miR396a (STTM396a/396a-88) resulted in earlier flowering and bigger fruits (Cao et al., 2016). In Arabidopsis, the peanut witches’ broom effector PHYLLODY SYMPTOMS 1 (PHYL1) interferes with miR396-mediated regulation of SHORT VEGETATIVE PHASE (SVP) transcript levels to control flower formation (Yang et al., 2015).
In this study, we also identified the bna-miR396–BnaA06g21030D/BnaA01g33410D pair. BnaA06g21030D (SER/ARG-rich protein 34A, SRp34a) encodes a member of the Ser-Arg-rich (SR) protein family, which plays multiple roles in post-transcriptional regulation of gene expression by alternative splicing (Sanford et al., 2004; Richardson et al., 2011). BnaA01g33410D is homologous to EMBRYO DEFECTIVE 3119 (EMB3119), which plays an important role in Arabidopsis growth and development (Meinke, 2020). Although many studies have focused on miR396, the roles of the bna-miR396–Bna.A06SRp34a/Bna.A01EMB3119 pair are largely unknown in the context of HI and should be studied in more detail.
Based on our collective results, we propose a potential model for the regulation of HI in B. napus during the seed-filling stage (Figure 6). Genetic differences and external environmental stimuli affect the expression of circadian rhythm-related genes, stress response genes, and miRNAs. miRNAs, in turn, may modulate circadian rhythm-related genes and stress response genes via a transcript cleavage. The resulting modulation of transcript levels of transporter activity-related genes will affect water and solutes flow to the seeds, driving seed development and the accumulation of storage compounds, thus determining seed yield and HI.
Figure 6. A proposed model for HI determination during the seed-filling stage.
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: BIG Data Center under BioProject accession number PRJNA602979.
Author Contributions
JL and KL designed the experiments. CZ, WC, XL, BY, LZ, and ZX performed the bioinformatics analysis. BY and LZ performed RNA extractions and qRT-PCR experiments. CZ and WC wrote the manuscript. All the authors reviewed the manuscript.
Funding
This work was supported by the National Natural Science Foundation of China (31830067 and 31871653), the Talent Project of Chongqing Natural Science Foundation (cstc2021ycjh-bgzxm0033), the 111 Project (B12006), Chongqing Graduate Student Research Innovation Project (CYB20102), Chongqing Special Postdoctoral Science Foundation (cstc2019jcyj-bsh0102), and Germplasm Creation Special Program of Southwest University.
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/fpls.2022.855486/full#supplementary-material
Footnotes
- ^ http://www.genoscope.cns.fr/brassicanapus/data/
- ^ http://www.mirbase.org/
- ^ http://rfam.xfam.org/
- ^ http://plantgrn.noble.org/psRNATarget/
- ^ http://systemsbiology.cau.edu.cn/agriGOv2/
- ^ https://biodb.swu.edu.cn/qprimerdb/
References
Aluri, S., and Büttner, M. (2007). Identification and functional expression of the Arabidopsis thaliana vacuolar glucose transporter 1 and its role in seed germination and flowering. Proc. Natl. Acad. Sci. U. S. A. 104, 2537–2542. doi: 10.1073/pnas.0610278104
Amanullah Shah, P. (2010). Timing and rate of nitrogen application influence grain quality and yield in maize planted at high and low densities. J. Sci. Food Agric. 90, 21–29. doi: 10.1002/jsfa.3769
Austin, R., Bingham, J., Blackwell, R., Evans, L., Ford, M., Morgan, C., et al. (1980). Genetic improvements in winter wheat yields since 1900 and associated physiological changes. J. Agric. Sci. 94, 675–689. doi: 10.1017/S0021859600028665
Basu, U., Bajaj, D., Sharma, A., Malik, N., Daware, A., Narnoliya, L., et al. (2019). Genetic dissection of photosynthetic efficiency traits for enhancing seed yield in chickpea. Plant Cell Environ. 42, 158–173. doi: 10.1111/pce.13319
Beche, E., Benin, G., Da Silva, C. L., Munaro, L. B., and Marchese, J. A. (2014). Genetic gain in yield and changes associated with physiological traits in Brazilian wheat during the 20th century. Eur. J. Agron. 61, 49–59. doi: 10.1016/j.eja.2014.08.005
Bennett, E. J., Roberts, J. A., and Wagstaff, C. (2011). The role of the pod in seed development: strategies for manipulating yield. New Phytol. 190, 838–853. doi: 10.1111/j.1469-8137.2011.03714.x
Bennett, E. J., Roberts, J. A., and Wagstaff, C. (2012). Manipulating resource allocation in plants. J. Exp. Bot. 63, 3391–3400. doi: 10.1093/jxb/err442
Boden, S. A., Cavanagh, C., Cullis, B. R., Ramm, K., Greenwood, J., Finnegan, E. J., et al. (2015). Ppd-1 is a key regulator of inflorescence architecture and paired spikelet development in wheat. Nat. Plants 1:14016. doi: 10.1038/nplants.2014.16
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Budak, H., and Akpinar, B. A. (2015). Plant miRNAs: biogenesis, organization and origins. Funct. Integr. Genomics 15, 523–531. doi: 10.1007/s10142-015-0451-2
Bustin, S. A., Beaulieu, J. F., Huggett, J., Jaggi, R., Kibenge, F. S., Olsvik, P. A., et al. (2010). MIQE précis: practical implementation of minimum standard guidelines for fluorescence-based quantitative real-time PCR experiments. BMC Mol. Biol. 11:74. doi: 10.1186/1471-2199-11-74
Büttner, M. (2007). The monosaccharide transporter (-like) gene family in Arabidopsis. FEBS Lett. 581, 2318–2324. doi: 10.1016/j.febslet.2007.03.016
Byrt, C. S., Zhao, M. C., Kourghi, M., Bose, J., Henderson, S. W., Qiu, J. E., et al. (2017). Non-selective cation channel activity of aquaporin AtPIP2;1 regulated by Ca2+ and pH. Plant Cell Environ. 40, 802–815. doi: 10.1111/pce.12832
Cao, D. Y., Wang, J., Ju, Z., Liu, Q. Q., Li, S., Tian, H. Q., et al. (2016). Regulations on growth and development in tomato cotyledon, flower and fruit via destruction of miR396 with short tandem target mimic. Plant Sci. 247, 1–12. doi: 10.1016/j.plantsci.2016.02.012
Chao, H. B., Raboanatahiry, N., Wang, X. D., Zhao, W. G., Chen, L., Guo, L. X., et al. (2019). Genetic dissection of harvest index and related traits through genome-wide quantitative trait locus mapping in Brassica napus L. Breed. Sci. 69, 104–116. doi: 10.1270/jsbbs.18115
Che, R. H., Tong, H. N., Shi, B. H., Liu, Y. Q., Fang, S. R., Liu, D. P., et al. (2015). Control of grain size and rice yield by GL2-mediated brassinosteroid responses. Nat. Plants 2:15195. doi: 10.1038/nplants.2015.195
Chen, H., Yang, Q., Chen, K., Zhao, S., Zhang, C., Pan, R., et al. (2019). Integrated microRNA and transcriptome profiling reveals a miRNA-mediated regulatory network of embryo abortion under calcium deficiency in peanut (Arachis hypogaea L.). BMC Genomics 20:392. doi: 10.1186/s12864-019-5770-6
Chen, L. Q., Lin, I. W., Qu, X. Q., Sosso, D., McFarlane, H. E., Londoño, A., et al. (2015). A cascade of sequentially expressed sucrose transporters in the seed coat and endosperm provides nutrition for the Arabidopsis embryo. Plant Cell 27, 607–619. doi: 10.3410/f.725398372.793505270
Chen, X. L., Jiang, L. R., Zheng, J. S., Chen, F. Y., Wang, T. S., Wang, M. L., et al. (2019). A missense mutation in Large Grain Size 1 increases grain size and enhances cold tolerance in rice. J. Exp. Bot. 70, 3851–3866. doi: 10.1093/jxb/erz192
Chen, Z. Y., Huo, Q., Yang, H., Jian, H. J., Qu, C. M., Lu, K., et al. (2019). Joint RNA-seq and miRNA profiling analyses to reveal molecular mechanisms in regulating thickness of pod canopy in Brassica napus. Genes 10:591. doi: 10.3390/genes10080591
D’Ario, M., Griffiths-Jones, S., and Kim, M. S. (2017). Small RNAs: big impact on plant development. Trends Plant Sci. 22, 1056–1068. doi: 10.1016/j.tplants.2017.09.009
De Angeli, A., Monachello, D., Ephritikhine, G., Frachisse, J. M., Thomine, S., Gambale, F., et al. (2006). The nitrate/proton antiporter AtCLCa mediates nitrate accumulation in plant vacuoles. Nature 442, 939–942. doi: 10.3410/f.1040148.489068
Demes, E., Besse, L., Cubero-Font, P., Satiat-Jeunemaitre, B., Thomine, S., and De Angeli, A. (2020). Dynamic measurement of cytosolic pH and [NO3–] uncovers the role of the vacuolar transporter AtCLCa in cytosolic pH homeostasis. Proc. Natl. Acad. Sci. U. S. A. 117, 15343–15353. doi: 10.1101/716050
Diepenbrock, W. (2000). Yield analysis of winter oilseed rape (Brassica napus L.): a review. Field Crops Res. 67, 35–49. doi: 10.1016/s0378-4290(00)00082-4
Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635
Eom, J. S., Chen, L. Q., Sosso, D., Julius, B. T., Lin, I. W., Qu, X. Q., et al. (2015). SWEETs, transporters for intracellular and intercellular sugar translocation. Curr. Opin. Plant Biol. 25, 53–62. doi: 10.1016/j.pbi.2015.04.005
Ercoli, M. F., Rojas, A. M. L., Debernardi, J. M., Palatnik, J. F., and Rodriguez, R. E. (2016). Control of cell proliferation and elongation by miR396. Plant Signal Behav. 11:e1184809. doi: 10.1080/15592324.2016.1184809
Fan, S. C., Lin, C. S., Hsu, P. K., Lin, S. H., and Tsay, Y. F. (2009). The Arabidopsis nitrate transporter NRT1.7, expressed in phloem, is responsible for source-to-sink remobilization of nitrate. Plant Cell 21, 2750–2761. doi: 10.3410/f.1165218.628157
Fang, Y. J., Xie, K. B., and Xiong, L. Z. (2014). Conserved miR164-targeted NAC genes negatively regulate drought resistance in rice. J. Exp. Bot. 65, 2119–2135. doi: 10.1093/jxb/eru072
Flis, A., Mengin, V., Ivakov, A. A., Mugford, S. T., Hubberten, H. M., Encke, B., et al. (2019). Multiple circadian clock outputs regulate diel turnover of carbon and nitrogen reserves. Plant Cell Environ. 42, 549–573. doi: 10.1111/pce.13440
Foulkes, M. J., Reynolds, M. P., and Sylvester-Bradley, R. (2009). Genetic improvement of grain crops-chapter 15: yield potential. Crop Physiol. 68, 355–385. doi: 10.1016/j.fcr.2012.05.014
Fu, T. D., and Zhou, Y. M. (2013). Progress and future development of hybrid rapeseed in China. Eng. Sci. 11, 13–18.
Fung-Uceda, J., Lee, K., Seo, P. J., Polyn, S., De Veylder, L., and Mas, P. (2018). The circadian clock sets the time of DNA replication licensing to regulate growth in Arabidopsis. Dev. Cell 45, 101–113.e4. doi: 10.3410/f.732902266.793569724
Gajića, B., Kresovićb, B., Tapanarovaa, A., Životić, L., and Todorovićc, M. (2018). Effect of irrigation regime on yield, harvest index and water productivity of soybean grown under different precipitation conditions in a temperate environment. Agric. Water Manage. 210, 224–231. doi: 10.1016/j.agwat.2018.08.002
Gao, P., Bai, X., Yang, L., Lv, D. K., Li, Y., Cai, H., et al. (2010). Over-expression of osa-MIR396c decreases salt and alkali stress tolerance. Planta 231, 991–1001. doi: 10.1007/s00425-010-1104-2
Giorgia, L., Irene, L. G., Meurer, J., Jürgen, S., and Serena, S. (2011). The Cytosolic Kinases STY8, STY17, and STY46 Are Involved in Chloroplast Differentiation in Arabidopsis. Plant Physiol. 157, 70–85. doi: 10.1104/pp.111.182774
Guo, C. K., Xu, Y. M., Shi, M., Lai, Y. M., Wu, X., Wang, H. S., et al. (2017). Repression of miR156 by miR159 regulates the timing of the juvenile-to-adult transition in Arabidopsis. Plant Cell 29, 1293–1304. doi: 10.3410/f.727650668.793535158
Ham, B. K., Li, G., Kang, B. H., Zeng, F., and Lucas, W. J. (2012). Overexpression of Arabidopsis plasmodesmata germin-like proteins disrupts root growth and development. Plant Cell 24, 3630–3648. doi: 10.1105/tpc.112.101063
Hay, R. K. M. (1995). Harvest index: a review of its use in plant breeding and crop physiology. Ann. Appl. Biol. 126, 197–216. doi: 10.1111/j.1744-7348.1995.tb05015.x
Hirner, B., Fischer, W. N., Rentsch, D., Kwart, M., and Frommer, W. B. (1998). Developmental control of H+ / amino acid permease gene expression during seed development of Arabidopsis. Plant J. 14, 535–544. doi: 10.1046/j.1365-313x.1998.00151.x
Huang, T. B., López-Giráldez, F., Townsend, J. P., and Irish, V. F. (2012). RBE controls microRNA164 expression to effect floral organogenesis. Development 139, 2161–2169. doi: 10.1242/dev.075069
Hütsch, B. W., and Schubert, S. (2017). Harvest index of maize (Zea mays L.): are there possibilities for improvement. Adv. Agron. 146, 37–82. doi: 10.1016/bs.agron.2017.07.004
Jian, H. J., Lu, K., Yang, B., Wang, T. Y., Zhang, L., Zhang, A. X., et al. (2016). Genome-Wide analysis and expression profiling of the SUC and SWEET gene families of sucrose transporters in oilseed rape (Brassica napus L.). Front. Plant Sci. 7:1464. doi: 10.3389/fpls.2016.01464
Kim, S. C., Nusinow, D. A., Sorkin, M. L., Pruneda-Paz, J., and Wang, X. M. (2019). Interaction and regulation between lipid mediator phosphatidic acid and circadian clock regulators. Plant Cell 31, 399–416. doi: 10.3410/f.734916270.793556708
Klepek, Y. S., Geiger, D., Stadler, R., Klebl, F., Landouar-Arsivaud, L., Lemoine, R., et al. (2005). Arabidopsis POLYOL TRANSPORTER5, a new member of the monosaccharide transporter-like superfamily, mediates H+-symport of numerous substrates, including myo-inositol, glycerol, and ribose. Plant Cell 17, 204–218. doi: 10.1105/tpc.104.026641
Klepek, Y. S., Volke, M., Konrad, K. R., Wippel, K., Hoth, S., Hedrich, R., et al. (2010). Arabidopsis thaliana POLYOL/MONOSACCHARIDE TRANSPORTERS 1 and 2: fructose and xylitol/H+ symporters in pollen and young xylem cells. J. Exp. Bot. 61, 537–550. doi: 10.3410/f.1430961.908059
Kourghi, M., Nourmohammadi, S., Pei, J. V., Qiu, J. E., McGaughey, S., Tyerman, S. D., et al. (2017). Divalent cations regulate the ion conductance properties of diverse classes of aquaporins. Int. J. Mol. Sci. 18:2323. doi: 10.3390/ijms18112323
Kuromori, T., Sugimoto, E., and Shinozaki, K. (2011). Arabidopsis mutants of AtABCG22, an ABC transporter gene, increase water transpiration and drought susceptibility. Plant J. 67, 885–894. doi: 10.1111/j.1365-313x.2011.04641.x
Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. doi: 10.3410/f.1156941.617005
Langmead, B., Trapnell, C., Pop, M., and Salzberg, S. L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10:R25. doi: 10.1186/gb-2009-10-3-r25
LeClere, S., Rampey, R. A., and Bartel, B. (2004). IAR4, a gene required for auxin conjugate sensitivity in Arabidopsis, encodes a pyruvate dehydrogenase E1alpha homolog. Plant Physiol. 135, 989–999. doi: 10.1104/pp.104.040519
Legnaioli, T., Cuevas, J., and Mas, P. (2009). TOC1 functions as a molecular switch connecting the circadian clock with plant responses to drought. EMBO J. 28, 3745–3757. doi: 10.3410/f.1166599.627609
Leng, S. H., Zhu, G. R., Deng, X. L., and Jiang, S. (1992). Studies on the sources of the dry matter in the seed of rapeseed. Acta Agron. Sin. 04, 250–257.
Li, N., Song, D. J., Peng, W., Zhan, J. P., Shi, J. Q., Wang, X. F., et al. (2019). Maternal control of seed weight in rapeseed (Brassica napus L.): the causal link between the size of pod (mother, source) and seed (offspring, sink). Plant Biotechnol. J. 17, 736–749. doi: 10.1111/pbi.13011
Li, Q. Q., Wu, G. X., Zhao, Y. P., Wang, B. B., Zhao, B. B., Kong, D. X., et al. (2020). CRISPRCas9-mediated knockout and overexpression studies reveal a role of maize phytochrome C in regulating flowering time and plant height. Plant Biotechnol. J. 18, 2520–2532. doi: 10.1111/pbi.13429
Li, S. C., Gao, F. Y., Xie, K. L., Zeng, X. H., Cao, Y., Zeng, J., et al. (2016). The OsmiR396c-OsGRF4-OsGIF1 regulatory module determines grain size and yield in rice. Plant Biotechnol. J. 14, 2134–2146. doi: 10.1111/pbi.12569
Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general-purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930. doi: 10.1093/bioinformatics/btt656
Liu, C., Song, G. Y., Zhou, Y. H., Qu, X. F., Guo, Z. B., Liu, Z. W., et al. (2015). OsPRR37 and Ghd7 are the major genes for general combining ability of DTH, PH and SPP in rice. Sci. Rep. 5:12803. doi: 10.1038/srep12803
Liu, H., Gu, F. W., Dong, S. Y., Liu, W., Wang, H., Chen, Z. Q., et al. (2016). CONSTANS-like 9 (COL9) delays the flowering time in Oryza sativa by repressing the Ehd1 pathway. Biochem. Biophys. Res. Commun. 479, 173–178. doi: 10.1016/j.bbrc.2016.09.013
Liu, W. W., Sun, Q., Wang, K., Du, Q. G., and Li, W. X. (2017). Nitrogen limitation adaptation (NLA) is involved in source-to-sink remobilization of nitrate by mediating the degradation of NRT1.7 in Arabidopsis. New Phytol. 214, 734–744. doi: 10.1111/nph.14396
Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2–Δ Δ CT method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Lu, K., Chai, Y. R., Zhang, K., Wang, R., Chen, L., Lei, B., et al. (2008). Cloning and characterization of phosphorus starvation inducible Brassica napus PURPLE ACID PHOSPHATASE 12 gene family, and imprinting of a recently evolved MITE-minisatellite twin structure. Theor. Appl. Genet. 117, 963–975. doi: 10.1007/s00122-008-0836-x
Lu, K., Li, T., He, J., Chang, W., Zhang, R., Liu, M., et al. (2018). qPrimerDB: a thermodynamics-based gene-specific qPCR primer database for 147 organisms. Nucleic Acids Res. 46, D1229–D1236. doi: 10.1093/nar/gkx725
Lu, K., Peng, L., Zhang, C., Lu, J., Yang, B., Xiao, Z., et al. (2017). Genome-wide association and transcriptome analyses reveal candidate genes underlying yield-determining traits in Brassica napus. Front. Plant Sci. 8:206. doi: 10.3389/fpls.2017.00206
Lu, K., Xiao, Z. C., Jian, H. J., Peng, L., Qu, C. M., Fu, M. L., et al. (2016). A combination of genome-wide association and transcriptome analysis reveals candidate genes controlling harvest index-related traits in Brassica napus. Sci. Rep. 6:36452. doi: 10.1038/srep36452
Luo, X., Ma, C. Z., Yue, Y., Hu, K. N., Li, Y. Y., Duan, Z. Q., et al. (2015). Unravelling the complex trait of harvest index in rapeseed (Brassica napus L.) with association mapping. BMC Genomics 16:379. doi: 10.1186/s12864-015-1607-0
Macho-Rivero, M. A., Herrera-Rodríguez, M. B., Brejcha, R., Schäffner, A. R., Tanaka, N., Fujiwara, T., et al. (2018). Boron toxicity reduces water transport from root to shoot in Arabidopsis plants. Evidence for a reduced transpiration rate and expression of major pip aquaporin genes. Plant Cell Physiol. 59, 836–844. doi: 10.1093/pcp/pcy026
Maere, S., Heymans, K., and Kuiper, M. J. B. (2005). BiNGO: a Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics 21, 3448–3449. doi: 10.1093/bioinformatics/bti551
Manara, A., DalCorso, G., Guzzo, F., and Furini, A. (2015). Loss of the atypical kinases ABC1K7 and ABC1K8 changes the lipid composition of the chloroplast membrane. Plant Cell Physiol. 56, 1193–1204. doi: 10.1093/pcp/pcv046
Manara, A., Dalcorso, G., Leister, D., Jahns, P., Baldan, B., and Furini, A. (2014). AtSIA1 and AtOSA1: two ABC1 proteins involved in oxidative stress responses and iron distribution within chloroplasts. New Phytol. 201, 452–465. doi: 10.1111/nph.12533
Mcclung, C. R., and Gutiérrez, R. A. (2010). Network news: prime time for systems biology of the plant circadian clock. Curr. Opin. Genet. Dev. 20, 588–598. doi: 10.1016/j.gde.2010.08.010
Meinke, D. W. (2020). Genome-wide identification of EMBRYO-DEFECTIVE (EMB) genes required for growth and development in Arabidopsis. New Phytol. 226, 306–325. doi: 10.1111/nph.16071
Mentzen, W. I., Peng, J. L., Ransom, N., Nikolau, B. J., and Wurtele, E. S. (2008). Articulation of three core metabolic processes in Arabidopsis: fatty acid biosynthesis, leucine catabolism and starch metabolism. BMC Plant Biol. 8:76. doi: 10.1186/1471-2229-8-76
Miao, C. B., Wang, D., He, R. Q., Liu, S. K., and Zhu, J. K. (2020). Mutations in MIR396e and MIR396f increase grain size and modulate shoot architecture in rice. Plant Biotechnol. J. 18, 491–501. doi: 10.1111/pbi.13214
Nakamichi, N., Kita, M., Niinuma, K., Ito, S., Yamashino, T., Mizoguchi, T., et al. (2007). Arabidopsis clock-associated Pseudo-Response Regulators PRR9, PRR7 and PRR5 coordinately and positively regulate flowering time through the canonical CONSTANS-dependent photoperiodic pathway. Plant Cell Physiol. 48, 822–832. doi: 10.1093/pcp/pcm056
Nakamichi, N., Kudo, T., Makita, N., Kiba, T., Kinoshita, T., and Sakakibara, H. (2020). Flowering time control in rice by introducing Arabidopsis clock-associated PSEUDO-RESPONSE REGULATOR 5. Biosci. Biotechnol. Biochem. 84, 970–979. doi: 10.1080/09168451.2020.1719822
Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., and Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods 14, 417–419. doi: 10.3410/f.727375548.793529669
Prasad, P. V. V., Boote, K. J., Allen, L. H., Sheehy, J. E., and Thomas, J. M. G. (2006). Species, ecotype and cultivar differences in spikelet fertility and harvest index of rice in response to high temperature stress. Field Crops Res. 95, 398–411. doi: 10.1016/j.fcr.2005.04.008
Qu, C. M., Fu, F. Y., Liu, M., Zhao, H. Y., Liu, C., Li, J. N., et al. (2015). Comparative transcriptome analysis of recessive male sterility (RGMS) in sterile and fertile Brassica napus lines. PLoS One 10:e0144118. doi: 10.1371/journal.pone.0144118
Qu, C., Zhao, H., Fu, F., Wang, Z., Zhang, K., Zhou, Y., et al. (2016). Genome-wide survey of flavonoid biosynthesis genes and gene expression analysis between black- and yellow-seeded Brassica napus. Front. Plant Sci. 7:1755. doi: 10.3389/fpls.2016.01755
Raman, S., Greb, T., Peaucelle, A., Blein, T., Laufs, P., and Theres, K. (2008). Interplay of miR164, CUP-SHAPED COTYLEDON genes and LATERAL SUPPRESSOR controls axillary meristem formation in Arabidopsis thaliana. Plant J. 55, 65–76. doi: 10.1111/j.1365-313x.2008.03483.x
Richardson, D. N., Rogers, M. F., Labadorf, A., Ben-Hur, A., Guo, H., Paterson, A. H., et al. (2011). Comparative analysis of serine/arginine-rich proteins across 27 eukaryotes: insights into sub-family classification and extent of alternative splicing. PLoS One 6:e24542. doi: 10.1371/journal.pone.0024542
Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616
Sadras, V. O., and Lawson, C. (2011). Genetic gain in yield and associated changes in phenotype, trait plasticity and competitive ability of south Australian wheat varieties released between 1958 and 2007. Crop Pasture Sci. 62, 533–549. doi: 10.1071/cp11060
Sanford, J. R., Gray, N. K., Beckmann, K., and Cáceres, J. F. (2004). A novel role for shuttling SR proteins in mRNA translation. Genes Dev. 18, 755–768. doi: 10.1101/gad.286404
Sawada, Y., Toyooka, K., Kuwahara, A., Sskata, A., Nagano, M., Satio, K., et al. (2009). Arabidopsis bile acid: sodium symporter family protein 5 is involved in methionine-derived glucosinolate biosynthesis. Plant Cell Physiol. 50, 1579–1586. doi: 10.1093/pcp/pcp110
Seo, P. J., Park, J. M., Kang, S. K., Kim, S. G., and Park, C. M. (2011). An Arabidopsis senescence-associated protein SAG29 regulates cell viability under high salinity. Planta 233, 189–200. doi: 10.1007/s00425-010-1293-8
Shan, J. W., Song, W., Zhou, J., Wang, X. H., Xie, C. H., Gao, X. X., et al. (2013). Transcriptome analysis reveals novel genes potentially involved in photoperiodic tuberization in potato. Genomics 102, 388–396. doi: 10.1016/j.ygeno.2013.07.001
Shen, J. X., Fu, T. D., Yang, G. S., Ma, C. Z., and Tu, J. X. (2010). Genetic analysis of rapeseed self-incompatibility lines reveals significant heterosis of different patterns for yield and oil content traits. Plant Breed. 124, 111–116. doi: 10.1111/j.1439-0523.2004.01069.x
Sommer, A., Geist, B., Da Ines, O., Gehwolf, R., Schäffner, A. R., and Obermeyer, G. (2008). Ectopic expression of Arabidopsis thaliana plasma membrane intrinsic protein 2 aquaporins in lily pollen increases the plasma membrane water permeability of grain but not of tube protoplasts. New Phytol. 180, 787–797. doi: 10.1111/j.1469-8137.2008.02607.x
Sun, H. M., Qian, Q., Wu, K., Luo, J. J., Wang, S. S., Zhang, C. W., et al. (2014). Heterotrimeric G proteins regulate nitrogen-use efficiency in rice. Nat. Genet. 46, 652–656. doi: 10.1038/ng.2958
Szczygieł-Sommer, A., and Gaj, M. D. (2019). The miR396-GRF regulatory module controls the embryogenic response in Arabidopsis via an auxin-related pathway. Int. J. Mol. Sci. 20:5221. doi: 10.3390/ijms20205221
Tang, J. Y., and Chu, C. C. (2017). MicroRNAs in crop improvement: fine-tuners for complex traits. Nat. Plants 3:17077. doi: 10.1038/nplants.2017.77
Tang, Y. Y., Liu, H. H., Guo, S. Y., Wang, B., Li, Z. T., Chong, K., et al. (2018). OsmiR396d miRNA affects gibberellin and brassinosteroid signaling to regulate plant architecture. Plant Physiol. 176, 946–959. doi: 10.1104/pp.17.00964
Tian, T., Liu, Y., Yan, H. Y., You, Q., Yi, X., Du, Z., et al. (2017). agriGO v2.0: a GO analysis toolkit for the agricultural community, 2017 update. Nucleic Acids Res. 45, W122–W129. doi: 10.1093/nar/gkx382
Tong, W., Imai, A., Tabata, R., Shigenobu, S., Yamaguchi, K., Yamada, M., et al. (2016). Polyamine resistance is increased by mutations in a nitrate transporter gene NRT1.3 (AtNPF6.4) in Arabidopsis thaliana. Front. Plant Sci. 7:834. doi: 10.3389/fpls.2016.00834
Trapnell, C., Roberts, A., Goff, L., Pertea, G., Kim, D., Kelley, D. R., et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578. doi: 10.3410/f.14267340.15779565
Unkovich, M., Baldock, J., and Forbes, M. (2010). Variability in harvest index of grain crops and potential significance for carbon accounting: examples from Australian agriculture. Adv. Agron. 105, 173–218. doi: 10.1016/s0065-2113(10)05005-4
Voinnet, O. (2009). Origin, biogenesis, and activity of plant microRNAs. Cell 136, 669–687. doi: 10.1016/j.cell.2009.01.046
Wang, C. Y., Zhang, S. C., Yu, Y., Luo, Y. C., Liu, Q., Ju, C. L., et al. (2014). MiR397b regulates both lignin content and seed number in Arabidopsis via modulating a laccase involved in lignin biosynthesis. Plant Biotechnol. J. 12, 1132–1142. doi: 10.1111/pbi.12222
Wang, H., Zhang, L. Y., Tao, Y., Wang, Z. D., Shen, D., and Dong, H. S. (2020). Transmembrane helices 2 and 3 determine the localization of plasma membrane intrinsic proteins in eukaryotic cells. Front. Plant Sci. 10:1671. doi: 10.3389/fpls.2019.01671
Wang, X. D., Zheng, M., Liu, H. F., Zhang, L., Chen, F., Zhang, W., et al. (2020). Fine-mapping and transcriptome analysis of a candidate gene controlling plant height in Brassica napus L. Biotechnol. Biofuels 13:42.
Xu, H., Li, X. F., Zhang, H., Wang, L. C., Zhu, Z. G., Gao, J. P., et al. (2020). High temperature inhibits the accumulation of storage materials by inducing alternative splicing of OsbZIP58 during filling stage in rice. Plant Cell Environ. 43, 1879–1896. doi: 10.1111/pce.13779
Yang, C. Y., Huang, Y. H., Lin, C. P., Lin, Y. Y., Hsu, H. C., Wang, C. N., et al. (2015). MicroRNA396-targeted SHORT VEGETATIVE PHASE is required to repress flowering and is related to the development of abnormal flower symptoms by the phyllody symptoms1 effector. Plant Physiol. 168, 1702–1716. doi: 10.1104/pp.15.00307
Yang, M., He, J. B., Wan, S. B., Li, W. Y., Chen, W. J., Wang, Y. M., et al. (2021). Fine mapping of the BnaC04.BIL1 gene controlling plant height in Brassica napus L. BMC Plant Biol. 21:359. doi: 10.1186/s12870-021-03137-9
Yang, X. H., Xia, X. Z., Zhang, Z. Q., Nong, B. X., Zeng, Y., Xiong, F. Q., et al. (2017). QTL mapping by whole genome re-sequencing and analysis of candidate genes for nitrogen use efficiency in rice. Front. Plant Sci. 8:1634. doi: 10.3389/fpls.2017.01634
Yang, X. Z., and Li, L. (2011). miRDeep-P: a computational tool for analyzing the microRNA transcriptome in plants. Bioinformatics 27, 2614–2615. doi: 10.1093/bioinformatics/btr430
Keywords: harvest index, transcriptome, miRNA, regulatory network, Brassica napus
Citation: Zhang C, Chang W, Li X, Yang B, Zhang L, Xiao Z, Li J and and Lu K (2022) Transcriptome and Small RNA Sequencing Reveal the Mechanisms Regulating Harvest Index in Brassica napus. Front. Plant Sci. 13:855486. doi: 10.3389/fpls.2022.855486
Received: 15 January 2022; Accepted: 22 February 2022;
Published: 04 April 2022.
Edited by:
Chunyu Zhang, Huazhong Agricultural University, ChinaReviewed by:
Libin Zhang, Huazhong University of Science and Technology, ChinaJinjin Jiang, Yangzhou University, China
Maoteng Li, Huazhong University of Science and Technology, China
Copyright © 2022 Zhang, Chang, Li, Yang, Zhang, Xiao, Li and Lu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jiana Li, ljn1950@swu.edu.cn; Kun Lu, drlukun@swu.edu.cn; orcid.org/0000-0003-1370-8633
†These authors have contributed equally to this work and share first authorship