We investigated variations in the genomic and transcriptional profiles of two B. napus lines, Zhongyou821 (ZY821) and Zhongshuang11 (ZS11), which had whole-genome re-sequencing and transcriptome data from two replicates of 11 tissues. These raw data were collected from the BIG Data Center under the BioProject accession codes PRJCA000376 and PRJCA001246 (Lu et al., 2019). To identify comprehensive transcription-level SNPs, we merged all RNA-seq datasets of the same sample for variation calling. Quality was controlled by using the Fastp software with default parameters (Chen et al., 2018).

A total of 6.37 and 7.56 Gb of genomic resequencing clean bases of ZY821 and ZS11 were aligned to the reference genome (B. napus cv. Darmor-bzh, a French winter variety) (Chalhoub et al., 2014) using BWA-MEM v0.7.15-r1140 with default parameters (Li and Durbin 2009). Mapped reads were sorted using Samtools (Li et al., 2009), and duplicate reads of the bam alignment file were marked using Picard tools (http://broadinstitute.github.io/picard). We then used the HaplotypeCaller module of the Genome Analysis Toolkit (GATK) to identify genetic variations in DNA with default parameters (McKenna et al., 2010). Variation calling and genotyping were performed across the Darmorbzh reference genome. We extracted SNPs using the SelectVariants module and filtered them using the VariantFiltration module of GATK with the following parameters: clusterSize, 3; clusterWindowSize, 10; filterExpression, QUAL <30.0 || MQ < 50.0 || QD < 2.

Transcriptional variations were identified by using 111.19 and 107.85 Gb RNA-seq clean bases of ZY821 and ZS11, and variants were called in RNA-seq data using a modified pipeline (https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels). We initially used STAR (Dobin et al., 2013) to construct index files of the reference genome and mapped clean RNA-seq reads to the reference genome with default parameters. Picard tools were then used to add ReadGroup and sort the mapped reads with “AddOrReplaceReadGroups SO = coordinate RGLB = mRNA.” Then, duplicate reads were marked using the MarkDuplicates module with the parameter settings with “CREATE_INDEX = true VALIDATION_STRINGENCY = SILENT.” We eliminated Ns using the SplitNCigarReads module in GATK to maintain grouping information with the parameters of rf as ReassignOneMappingQuality, RMQF 255, RMQT 60, and -U ALLOW_N_CIGAR_READS. We then identified transcriptional variations using the parameters of “-dontUseSoftClippedBases-stand_call_conf 20” in HaplotypeCaller. Finally, transcription-level SNPs were extracted and filtered by GATK using the same module and parameters as the genomic SNPs. All SNP datasets were annotated using SnpEff v4.1 (Cingolani et al., 2014) with upstream and downstream interval lengths set to 0.

To investigate the differences between the two subgenomes of B. napus, we compared SNP transmission between ZY821 and ZS11. Two samples or two levels containing a common homozygous SNP were regarded as having been transmitted. We were initially concerned with the transmission of genomic SNPs at the transcriptional level in ZY821 and ZS11. Because ZY821 is one of parents of ZS11, some ZY821 genetic components were transferred to the ZS11 genome during breeding. We further investigated the subgenome differences in SNPs transmitted from the ancestral parent to the offspring. We analyzed two types of transmission from ZY821 to ZS11, they were genomic SNP of ZY821 transferred to ZS11 genomic level and ZY821 genomic SNP transferred to ZS11 genomic level then transferred to ZS11 transcriptional level, respectively. The transmission ratio/rate was defined as the number of transmitted SNPs divided by the total number of genomic SNPs detected in the gene. In order to study the asymmetric distribution of transmitted SNPs in two subgenomes of B. napus, we analyzed the transmitted SNPs of all syntenic genes between A_n and C_n subgenomes (Chalhoub et al., 2014). Ratios of transmitted SNPs on A_n and C_n in different samples/levels were compared using two-tailed t-tests. A homozygous variant that arises from DNA to RNA transmission is regarded as a novel SNP variation in DNA to RNA.

Pathways and biological functions of genes without SNPs and genes with genomic SNPs that were all transmitted to the transcriptional level were explored via GO enrichment analyses using the PlantRegMap database (http://plantregmap.gao-lab.org/go.php) with the parameter settings of species, B. napus; aspects, biological process, molecular function, cellular component; threshold, and q ≤ 0.01. The GO enrichment results are shown on the ImageGP website (http://www.ehbio.com/ImageGP/).

Filtered high-quality clean RNA-seq data from different tissues and periods of ZY821 and ZS11 were mapped to the “Darmor-bzh” genome (Chalhoub et al., 2014) using STAR with default parameters (Dobin et al., 2013). The mapped reads with mapping quality (MQ) ≤ 30 were filtered using Samtools, and bam files were sorted. As TPM (transcripts per kilobase of exon model per million mapped reads) could correct the inconsistences while comparing the RNA-seq abundance among independent samples (Wagner et al., 2012), StringTie (Pertea et al., 2015) was used to count unique and normalized mapped reads as TPM for each gene with the parameter of “-B -e -G”.

We used linear regression analyses to investigate correlations between gene expression level and SNP numbers in gene regions, new SNP variations in DNA to RNA, and SNP transmission rates, respectively. Genes with different transmission rates were divided into 10 groups at 0.1 intervals. Genes with different numbers of SNP variations were divided into groups of 10 intervals in ZY821 DNA, ZY821 RNA, and ZS11 RNA, whereas those in ZS11 RNA were arranged into groups with 20 intervals since SNP variations occurred frequently. The linear model is y_i = ax_i+b, where y_i is the average TPM of the ith group, and x_i is a different group.

We filtered 6.37 Gb of genome re-sequencing, and 111.19 Gb RNA-seq data were used to identify SNP variation in ZY821. In total, 1,178,526 SNPs were detected in ZY821 genomic DNA, along with 19 chromosomes and 21 scaffolds (Supplementary Table S1). We found that 83,179 SNPs in chromosome chrC03 had the most mutations. Among these SNPs, there were 349,509 SNPs distributed in 40,903 genic regions and 846,691 SNPs located in intergenic regions. 162,124 of the genic SNPs were located in exon regions and 139,244 SNPs in intron regions (Supplementary Table S2). Among the SNPs in exon regions, 69,526 missense and 91,469 synonymous variations were annotated. We also identified 792,189 SNPs at the transcription level in ZY821, and most variations were detected in chromosome chrA03 (63,925 SNPs). We found 536,701 and 144,456 variations in the exon and intron regions, respectively, and 64.05% (349,110) and 34.63% (185,846) of SNPs in exon regions were annotated as synonymous and missense variants, respectively (Supplementary Table S2).

After quality control, we identified 1,823,207 and 686,858 genomic and transcriptional SNPs, respectively, in ZS11 using 7.56 Gb genome re-sequencing and 107.85 Gb RNA-seq data. Among the genomic SNPs, 1,352,662 were located in intergenic regions with the most variations, and 228,122 and 199,500 were located in exons and introns, respectively. Overall, 97,806 (42.87%) and 128,632 (56.39%) SNPs in exons were annotated as missense and synonymous variants, respectively. Among the SNPs at the transcription level, 464,923 and 126,513 variations occurred in the exon and intron regions, respectively. Among the SNPs in exons, 163,198 (35.10%) and 300,178 (64.57%) were missense and synonymous variants, respectively (Supplementary Table S2).

We identified 27,024 shared genes without SNPs at the genomic and transcriptional levels in both ZY821 and ZS11, among which 8,756 and 18,268 were located on A_n and C_n, respectively. Results of GO enrichment analysis of the genes without SNPs showed that 54 terms were enriched in cellular components, biological processes, and molecular functions. Almost all terms were basic biological processes for growth, such as an integral component of membrane and membrane part in cellular component, the regulation of nucleic acid-templated transcription and nucleoside bisphosphate metabolic processes in biological processes, enzyme inhibitor activity, and electron carrier activity in molecular function (Supplementary Table S3), which indicated that these genes were conserved.

The asymmetric characteristics of the A_n and C_n subgenomes of B. napus have been identified (Qian et al., 2014; Lu et al., 2019; Wu et al., 2019). We found that ZY821 and ZS11 had more total SNPs in C_n than in A_n at the genomic level, whereas C_n had fewer SNPs than A_n at the transcriptional level (Supplementary Table S1). The total length of A_n was 315.05 Mb in the reference genome, which was less than that in the C_n (526.93 Mb) subgenome; therefore, we assessed SNP density to reevaluate differences between the two subgenomes. Figure 1A shows the distribution of SNP densities in the whole genome. The density of ZS11 genomic SNPs was the highest at 3.7/kb in scaffold A10_random and that of ZS11 transcriptional SNPs was the lowest (0.18/kb) in the scaffold Unn_random. These might have been due to the frequency of coding genes being the lowest (1/15.61 kb) in the scaffold Unn_random. The mean SNP density in the genome was 1.39 and 0.93/kb at the ZY821 genomic and transcriptional levels, and 2.14, and 0.81/kb at the ZS11 genomic and transcriptional levels (Supplementary Table S1). The genomic and transcriptional SNP density was higher in A_n than in C_n, especially at the transcriptional level, and the SNP density of A_n was almost twice that of C_n (Figure 1B). Taken chromosomes chrA02 and chrC02 as an example due to they had good collinearity, the SNP density of chrA02 was generally higher than that of chrC02, although that of chrA02 was occasionally lower than that of chrC02 in some chromosomal regions (Figure 1C).

We initially detected 61,832 genomic SNPs that were transmitted at the transcriptional level in ZY821 and involved 25,583 coding genes. A total of 37,845 and 23,920 SNPs occurred in A_n and C_n, respectively (Figure 2A), whereas 14,258 synonymous SNPs in A_n while 6,952 synonymous SNPs in C_n. Among syntenic genes in A_n and C_n, the SNP transmission ratio from DNA to RNA was 31.05% in A_n and 30.23% in C_n (Figure 2B). In ZS11, 113,752 SNPs were transmitted from DNA to RNA in A_n, and 63,245 in C_n included 35,774 coding genes (Figure 2C). The ratios of transmitted syntenic genes in A_n and C_n were 55.43 vs. 43.21% (p = 5.45 × 10⁻⁸⁹; Figure 2D). Parental-related transmitted SNPs differed between the two subgenomes. Among 492,916 genomic SNPs from ZY821 transmitted to ZS11, in total 249,185, 241,871, and 1,860 were in A_n, C_n, and the scaffold Unn_random, respectively. These findings indicated that 44.06% of genomic SNP in ZY821was transmitted to the A_n of ZS11, and 40.10% of SNPs in C_n of ZY821 were transmitted to C_n of ZS11 (Figure 2E). Moreover, 76,188 SNPs detected at the ZY821 genomic level were transmitted to the ZS11 genome, followed by the ZS11 transcriptional level; of these, 48,367 (63.48%) of 76,188 occurred in A_n and 27,751 (36.42%) in C_n, respectively (Figure 2F). These results indicated that SNPs in A_n were more prone to transmission than those in C_n, especially transmission SNPs from the genomic level transfer to the transcriptional level.

The results of GO enrichment analyses revealed the functions of genes with a high ratio of SNP transmission. All (100%) genomic SNPs in 3,967 genes were transmitted to the transcriptional level in ZY821. The number of transmitted genes did not substantially differ between A_n and C_n (1,987 vs. 1,971). Only three GO terms were enriched in the biological process, viz. response to endogenous stimulus, response to hormone, and response to organic substances (Supplementary Table S4). Excluding biological processes, no GO terms were enriched in molecular function and cellular components (q < 0.01). Among 9,018 completely transmitted genes detected in ZS11, 5,171 and 3,835 were located in A_n and C_n, respectively. Eight GO terms were enriched in biological processes, four in cellular components, and two in molecular functions (Figure 3; Supplementary Table S5), which are involved in fundamental biological processes and stress responses. All three GO-enriched terms in the biological process of ZY821 were also enriched in ZS11, indicating that these genes play essential roles in the speciation of ZS11.

We analyzed whole gene expression levels of ZY821 and ZS11 to determine correlations between rates of SNP transmission from DNA to RNA and gene expression. We merged RNA-seq data from different tissues and periods and found that 14,717 and 15,822 genes were not expressed in ZY821 and ZS11, respectively, and 10,931 genes of them were shared, which might be because they were tissue- and time-specific. Gene expression levels decreased as the number of genomic SNPs in genes increased (Figures 4A,B), indicating that the expression of a gene with more SNPs would be lower in the B. napus genome. The expression of genes containing more RNA variations was lower in ZY821 and ZS11 during transcription (Figures 4C,D), which was consistent with the number of variations at the genomic and gene expression levels. However, SNP Figures 4E,F density was higher in A_n (Figure 1B), but gene expression was higher in C_n (Figures 4E,F). Further results revealed that transmitted SNPs had significantly associated with gene expression level. We found that 100% of 3,967 and 9,018 genes with SNPs from the genomic level were transmitted to the transcriptional level in ZY821 and ZS11, respectively. We established a linear regression model to fit the details of the SNP transmission ratio from DNA to RNA and gene expression. SNP transmission rates were significantly and positively associated with the gene expression levels. The R ² values of the fitted line were 0.689 and 0.735 in ZY821 and ZS11, respectively (Figures 4G,H), indicating that gene expression increased as the SNP transmission ratio from the DNA to the RNA level increased. These results showed that the ratio of SNPs transmitted from DNA to RNA correlated well with gene expression. The rates and numbers of transmitted SNPs were higher in A_n than in C_n (Figure 2), signifying asymmetric gene expression in polyploid B. napus.