Overview of sRNA sequencing results
We have already upload data and ensure the deposited data is made public (NCBI) : https://www.ncbi.nlm.nih.gov/sra/PRJNA760803, accession number:PRJNA760803. A total of 22,744,964 (low oleic acid rapeseed materials, A) and 26,060,122 (HOAR materials, B) raw reads were generated from the sequencing machine. After removing the adaptor sequences, filtering out low quality tags, and cleaning up sequences derived from adaptor ligation, 20,912,776 (A) and 23,710,938 (B) clean reads were obtained. Consequently, the bioinformatic analysis of these clean reads were carried out (Table 1). The size distribution patterns of the original and unique reads were displayed in Figure 1A. Small RNAs (24 nt) were the most abundant in all the samples. In addition, the clean reads exhibited 87.45% (A) and 88.26% (B) homology with the reference genome sequences. The sRNA sequencing results were of high quality and reliable and can be used for further functional analysis.
Differentially expressed miRNAs identification and their function analysis
The clustering analysis method was used to investigate the similarity between samples by calculating the differential miRNA distance between low oleic acid rapeseed (A) and HOAR (B) (Figure 1B). As shown, three samples of high or low oleic acid contents were found in a cluster, suggesting that these miRNAs might have similar biological functions.
After applying a P value < 0.05 and an absolute value of log2 (treatment/control) greater than 1.5 to identify differentially expressed miRNAs, 21 differentially expressed miRNAs (8.39% of the total) were detected (Table 2, Supplemental Table 2); their frequencies were calculated using TPM (Table 3). Among them, nine genes (42.86%) were up-regulated and 12 (57.14%) were down-regulated in HOAR materials (Figure 1C). Briefly, 17 miRNAs (including NC_027774.1_24533, NC_027757.1_266, bna-miR156b>bna-miR156c>bna- miR156g(bna-miR156b>c>g), NC_027774.1_24533*, NC_027769.1_17408*, bna-miR156a, bna-miR162a, NC_027769.1_17408, bna-miR172d, bna-miR824, bna-miR166f, bna-miR396a, bna-miR169m, bna-miR160a>bna-miR160b>bna-miR160c>bna-miR160d (bna-miR160a>b>c>d), NC_027760.1_5272, NC_027768.1_15573, and NW_013650328.1_26640*) and 358 target genes were obtained.
The target genes were then subjected to GO functional and KEGG Pathway analyses. In many cases, multiple terms were assigned to the same miRNA. Thus, 133 putative target genes were associated with 21 differentially expressed miRNAs and distributed into the following subcategories: 62 “Biological process”, 32 “Cellular component”, and 39 “Molecular function” (Figure 1D). Under “Biological process”, most of the target genes were related to “transcription”, and “regulation of transcription”. Within the “Cellular component” category, “nucleus” and “cytosol” were observed as much as “cytoplasm”. Among genes in the “Molecular function” category, most potential functions were related to “transcription factor activity”, “DNA binding”, and “Catalytic activity”. The distribution of target genes indicated that rapeseed underwent active metabolization.
The KEGG pathway enrichment analysis indicated that 20 pathways involving 36 target genes were enriched, including “alpha-Linolenic acid metabolism” (7 target genes), “Phagosome” (3 target genes), “Oxidative phosphorylation” (3 target genes), “Oxidative phosphorylation” (3 target genes), and “Protein processing in endoplasmic reticulum” (3 target genes). The fifteen target genes may be related to fatty acid metabolism: “alpha-linolenic acid metabolism” (7 genes), “Oxidative phosphorylation” (3 genes), “Carbon metabolism” (2 genes), “Citrate cycle (TCA cycle)” (1 gene), “Glycerolipid metabolism” (1 gene), and “Glycolysis / Gluconeogenesis” (1 gene). The top 20 KEGG enrichments (Figure 2A) show that α-Linolenic acid metabolism is the most significant, suggesting that bna00592 KEGG pathway may be involved in fatty acid metabolism in rapeseed (Figure 2B).
Expression pattern of bna-miR156b>c>g gene was detected by RT-PCR
To confirm the results of the miRNA sequence analysis, 21 annotated differentially expressed miRNAs were compared to the B. napus genome using BLAST [33] (Figure 2B).
Most expression trends of the RT-PCR analysis results agreed with the miRNA sequencing data (NC_027757.1_266, NC_027760.1_5272, NC_027760.1_5272*, NC_027761.1_6665, NC_027768.1_15573, NC_027769.1_17408, NC_027769.1_17408*, NC_027774.1_24533, NC_027774.1_24533*, NW_013650328.1_26640*, bna-miR162a, bna-miR167a>b, bna-miR169m, bna-miR172d, bna-miR396a, bna-miR824.), In addition, among the 13 miRNAs with significant difference, 9 had target genes (NC_027760.1_5272, NW_013650328.1_26640*, bna-miR156b>c>g, bna-miR166f, bna-miR169m, bna-miR396a, bna-miR824, bna-miR156a, and bna-miR160a>b>c>d), which may be the novel miRNAs related to fatty acids.
Moreover, the expressions of fatty acid metabolism related to differential miRNAs, such as bna- miR396a, bna-miR156b>c>g, and their target genes, were studied in different developmental stages (Figure 2C). The bna-miR396a has opposite expression pattern with its target gene, at first, the bna-miR396 had up-regulated expression, until the bud stage reached the peak, and the expression decreased with the growth stage; bna-miR156b>c>g had an opposite expression pattern with its target gene, it had down-regulated expression in the whole growth stages, and the expression decreased with the growth stages. Differentially expressed miRNAs and their target genes were related to fatty acid metabolism in bna-mi156b>c>g at different developmental stages (Figure 1C). In contrast to the expression pattern of bna-miR156b>c>g and its target gene, the expression level was down-regulated throughout the whole growth stage of the plant, and the expression level gradually decreased with the growth process of rapeseed.
Cloning of OPR genes in rapeseed and bioinformatic analysis
Target gene: bna-miR156b>c>g was cloned by miRNA sequencing, and four copies were detected: GSBRNA2T00012422001, GSBRNA2T00135385001, GSBRNA2T00082938001, and GSBRNA2T00094910001, which were named OPR1, OPR2, OPR3, and OPR4, respectively. Both OPR1 and OPR3 were 1119 bp, OPR2 and OPR4 were 1125 bp and 1122 bp, respectively (Figure 3A).
DNAMAN 7.0 software was used to compare the cloned target sequence with the rapeseed sequence published on the Brassica napus Genome Browser website. Different base position (Figure 3B), homology were more than 99% with the published sequences (Figure 3C). OPR1, OPR3, OPR4 had no base difference with the published sequence, there were 10 base differences between OPR2 and published sequence and the homology was 99.11%. Preliminary identification of OPR2 and OPR3 were located in A genome and OPR4 and OPR1 in C genome was conducted.
The number of four copies of OPR gene amino acids ranged from 372 bp to 374 bp, with a molecular weight of about 41 ku; The encoded amino acids were acidic (< 7), unstable (< 40), exhibited a fat coefficient of about 75, and belonged to fat-soluble proteins. Predicted subcellular localization of proteins were encoded by different copies of OPR genes and we found that these four proteins were located in the cytoplasm. We found that the four copies were all extracellular proteins without a transmembrane structure. Predicted protein secondary structures were summarized in Table 4.
The predicted tertiary structure model of the protein shows that the tertiary structures of OPR1, OPR2, OPR3, and OPR4 were all adapted to the 12-O-plant dienoate reductase model (integrated with the crystal structure of At1g76680 protein in A. thaliana), but their conformations were slightly different (Figure 3D). It has been reported that the protein structure A. thaliana of At1g76680 is similar to that of yeast ScOYE1 [34].
Vector construction and transformation of A. thaliana.
Using the synthesized cDNA as template, four target gene specific fragments were amplified by PCR with high fidelity, and the RNAi fragments with the same length were obtained (Figure 3E). Then, the target genes were recombined with the overexpression vector pCAMBIA1300-35s (OPR-OE), positive strains were screened. After sequencing, the overexpression vectors pCAMBIA1300-35s-OPR1 (OPR1-OE), pCAMBIA1300-35s-OPR2 (OPR2-OE), pCAMBIA1300-35s-OPR3 (OPR3-OE), and pCAMBIA1300-35s-OPR4 (OPR4-OE) were obtained (Figure 4A-B).
The RNAi fragment was recombined with the RNAi vector pCAMBIA1300-RNAi (OPRi) and the positive strains were screened. After sequencing, the RNAi vectors pCAMBIA1300-RNAi-OPR1 (OPR1i), pCAMBIA1300-RNAi-OPR2 (OPR2i), pCAMBIA1300-RNAi-OPR3 (OPR3i), and pCAMBIA1300-RNAi-OPR4 (OPR4i) were obtained (Figure 4C).
The recombinant vector was transformed into A. thaliana and the obtained A. thaliana was detected. The T1 transgenic A. thaliana seeds were screened using hygromycin (Figure 5A) and the results of hygromycin primer identification (Figure 5B) showed that each copy of A. thaliana OPR had been successfully transformed, and the target plasmid T-DNA had been inserted into the genome of A. thaliana.
Analysis of fatty acid composition and fatty acid content
The transformation methods in A. thaliana, reference Clough and Bent (1998) [35]. fatty acid composition was detected by gas chromatography [36]. we obtained fatty acid composition results of the T1 and T2 generations. Fatty acid composition of A. thaliana T-DNA insertion lines in table 5, and transformation materials fatty acid composition in the contrast T1/T2. The contents of oleic acid and stearic acid in OPR1i were significantly increased and the LA content decreased significantly; OPR1-OE will lead to a significant increase in palmitic acid and LA content; OPR2i significantly increased stearic acid content and decreased LA content; and OPR2-OE increased LA content.
The LA content in OPR3i decreased significantly; OPR3-OE significantly increased LA content; OPR4i significantly decreased LA content; and OPR4-OE increased LA content significantly.
Each copy of OPR-OE increased LA content, with an average increase of 12.56% in T1 generation and 7.185% in T2 generation. Subsequently, LA content in OPRi gene was significantly decreased, with an average decrease of 5.98% in T1 generation and 0.86% in T2 generation.
As shown in Table 5, oleic, linolenic, arachidonic, and erucic acids with the same variation trend as that of the fatty acid composition were selected for variance analysis. The results (Table 6) showed that the linolenic acid content in OPR1i significantly increased, while OPR2i significantly increased the linolenic acid content. Both OPR3-OE and OPR4-OE affected the content of arachidonic acid, which decreased significantly. In addition, OPR4i had no significant effect on the arachidonic acid content.