Mining for genes related to pistil abortion in Prunus sibirica L.

Plant materials

Seven-year-old P. sibirica clones in the National Forest Germplasm Resource Preservation Repository for Prunus species of Shenyang Agricultural University (Kazuo, Liaoning, China) were selected as the experimental materials. They were grown with a row spacing of 3m and an in-row distance of 2m. Clone No. 595 (experimental group) was selected as the material with abortive pistils (APt), and clone No. 28 (control group) was selected as the material with normal pistils (NPt). Three plants with good and uniform growth conditions were selected from each clone, and flower bud samples of the white stage (about 30 buds per sample) were collected for transcriptome sequencing on March 27, 2021. There were no extreme climatic conditions such as low temperature, high temperature or flood in the test site. For the convenience of description, the sequencing results from three repeats each of P. sibirica with abortive pistils and normal pistils were denoted as APt1, APt2, APt3, and NPt1, NPt2, NPt3, respectively.

Total RNA extraction and quality control

Total RNA was extracted from mixed flower buds samples from 3 plants using the EASYspin Plant microRNA Kit (Aidlab, China). RNA concentration and purity were measured using a NanoDrop 2000 (Thermo Fisher Scientific, Wilmington, DE, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit with the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA).

Library preparation and transcriptome sequencing

After the samples were tested for their quality, library construction for sequencing was performed. The detailed process was as follows: First, mRNA was isolated from total RNA by Oligo(dT)-attached magnetic beads, and was randomly fragmented in the fragmentation buffer. Then, the fragmented mRNA and random hexamer primer were used to synthesize the first-strand cDNA. The second-strand cDNA was synthesized by adding buffer, dNTPs, RNase H, and DNA polymerase I. The cDNA was purified using AMPure XP beads. Next, the double-stranded cDNA was end-repaired and A-tailed for adapter-ligated. AMPure XP beads were used to select 300∼400 bp fragments. Finally, cDNA libraries were obtained by PCR enrichment.

Qubit 2.0 and Agilent 2100 systems were used to examine cDNA concentration and insert size to ensure the library quality. The high-quality libraries with concentration greater than 2 nM were obtained by qPCR. Three biological replicates were used for each group for sequencing. The constructed libraries were sequenced on an Illumina Novaseq 6000 platform, 150 bp paired-end reads were generated, and the amount of data was about 6G clean data for each sample. The statistical power analysis of this experimental design was calculated using the R package RNASeqPower (Table S1).

Raw data quality control

Raw data (raw reads) of fastq format were first processed using in-house Perl scripts. In this step, clean data (clean reads) were obtained from the raw data by removing reads containing adapters, reads containing poly-N, and low-quality reads. At the same time, Q20, Q30, GC content and sequence duplication level of the clean data were calculated. All downstream analyses were based on high-quality clean data.

Sequence alignment

Sequence alignment and follow-up analysis were conducted using the P. sibirica genome (https://www.rosaceae.org/Analysis/10254124) as a reference genome. Hisat2 tools soft (Kim, Langmead & Salzberg, 2015) were used to map the reference genome.

Correlation assessment of replicates

Pearson’s correlation coefficient was used to evaluate the reproducibility of biological replicates (Liu et al., 2018). A closer R² value of 1 indicates better reproducibility between the two samples.

Gene functional annotation

In order to obtain the annotation information of unigenes, the unigene sequences were compared against the NCBI non-redundant protein sequences (Nr), NCBI non-redundant nucleotide sequences (Nt), protein family (Pfam), clusters of orthologous groups of proteins (KOG/COG), Swiss-Prot protein sequence database, KEGG ortholog (KO) and Gene Ontology (GO).

Quantification of gene expression levels

Quantification of gene expression levels was estimated by fragments per kilobase of transcript per million fragments mapped.

DEGs and transcription factor analysis

FPKM was applied to measure the expression level of a gene or transcript by StringTie using maximum flow algorithm (Florea, Song & Salzberg, 2013). Raw counts were input, and low expression was filtered. Differential expression analysis was performed using the edgeR (Robinson, McCarthy & Smyth, 2010). The analysis started by normalizing the input count. The false discovery rate (FDR) <0.01 & fold change ≥2 were set as the threshold for significantly differential expression. GO enrichment analysis of the DEGs was performed using the GOseq R packages based on Wallenius non-central hyper-geometric distribution (Young et al., 2010). KOBAS software (Mao et al., 2005) was used to test the statistical enrichment of DEGs in KEGG pathways. Transcription factor analysis was performed using iTAK 1.2 software. The genes sets of KEGG pathway and GO terms on molecular functions, biological processes and cellular components were employed in gene set enrichment analysis (GSEA). The log2FC of each differential group was used as the score of the background gene set to analyze the enrichment of the gene set. Enriched gene sets were identified as p-value<0.001 and FDR<0.05 (Khan et al., 2017).

Validation of transcriptome sequencing results

Template cDNA was generated using PrimeScript™ RT Master Mix (Perfect Real Time) (RR036A; TaKaRa Bio, Shiga, Japan) for quantitative real-time PCR (qRT-PCR). The PCR reactions for reverse transcription (10 µL total volume) contained 5 × PrimeScript RT Master Mix (2 µL), RNase Free dH₂O (7 µL), and total RNA (1 µL). The reaction conditions were 37 °C for 15 min, 85 °C for 5 s, and 4 °C for termination. Gene-specific primers (Table S2) were designed using Primer Premier 5.0. qRT-PCR was performed in three biological replicates. The gene primers used in the qRT-PCR experiments were listed in Table S2. 18SrRNA was used as a reference gene for qRT-PCR (Jin, 2018). qRT-PCR was performed using a StepOnePlus real-time PCR system. qRT-PCR was performed using TB Green^® Premix Ex Taq™ II (Tli RNaseH Plus) (TaKaRa Bio, Shiga, Japan). The PCR reactions (20 µL total volume) contained 10 µL 5 × TB Green Premix Ex Taq II (Tli RNaseH Plus), 6 µL RNase-free water, 0.8 µL upstream and downstream primers each, 2 µL cDNA template, 0.4 µL ROX reference Dye (50×). qRT-PCR was performed using three-step qPCR. The reaction conditions were 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, 55 °C for 30 s, and 72 °C for 30 s.

RNA isolation, library construction, and sequencing

RNA concentration of each sample ranged from 331.5 ng/ul to 925.2 ng/ul, OD260/280 ranged from 2.10 to 2.18, OD260/230 ranged from 1.16 to 2.21, RIN value ranged from 8.7 to 9.2, 28S/18S ranged from 1.94 to 2.09 (Table S3). The results showed that the quality of obtained total RNA was satisfying, and could meet the experimental requirements.

After quality control of sequencing data, 44.96 Gb clean data was obtained. Clean reads ranged from 22,462,256 to 27,048,846, clean bases ranged from 6,713,330,540 to 8,096,571,394, GC Content ranged from 45.89% to 46.03%, Q30 value ranged from 92.01% to 92.92% (Table S4). The results showed that the amount of data met the quality requirements for subsequent analysis.

Power analysis

The mean power of all DEGs detected with FDR <0.01 and fold change ≥2 was 56.85%. The power value of 19.54% DEGs was above 60.00%, and the power value of 93.33% DEGs was above 50.00% (Table S1).

Correlation assessment of biological replicates

The Pearson correlation coefficients among the three biological replicates of each experimental group were 0.986 (APt1/APt2), 0.977 (APt1/APt3), and 0.982 (APt2/APt3). The Pearson correlation coefficients among the three biological replicates of each control group were 0.973 (NPt1/NPt2), 0.983 (NPt1/NPt3), and 0.972 (NPt2/NPt3) (Fig. 1). The correlations between the samples were all greater than 0.9, demonstrating an excellent internal consistency, which met the requirements for further biological analysis.

The results of sequence alignment

Compared to the reference genome, the mapped rate of the experimental group samples ranged from 92.87% to 93.04%, and that of the control group samples ranged from 92.38% to 93.13%. The rate of reads map to ‘+’ and the rate of reads map to ‘−’ was similar (Table S5). The results showed that the quality and data volumes of transcriptome sequencing were relatively high, which met the requirements for subsequent data assembly and processing.

Gene function annotation

A total of 30,627 unigenes were annotated in least one database, accounting for 89.9% of the total unigenes. The highest annotation rate was obtained in the NR database, which was assigned 30,543 (89.7%) unigenes. In the GO database, 23,405 (68.7%) unigenes were annotated. In the eggNOG database, 23,357 (68.6%) unigenes were annotated. In the Pfam database, 21,616 (63.5%) unigenes were annotated. In the KEGG database, 18,941 (55.6%) unigenes were annotated. In the Swiss-Prot database, 18,634 (54.7%) unigenes were annotated. In the KOG database, 13,897 (40.8%) unigenes were annotated. In the COG database, 8256 (24.2%) unigenes were annotated (Table S6).

Differential gene expression analysis

A total of 1950 DEGs related to pistil abortion in P. sibirica were identified, including 1,000 upregulated genes and 950 downregulated genes (Fig. 2).

GO enrichment analysis of DEGs

The 1,422 identified DEGs were annotated and were assigned, 20 terms belong to biological process, 18 terms belong to cellular component, and 16 terms belong to molecular function. Biological processes included predominantly metabolic process (597 genes), cellular process (575 genes), single-organism process (468 genes), biological regulation (251 genes), response to stimulus (235 genes), localization (132 genes), signaling (116 genes), cellular component organization or biogenesis (87 genes), multicellular organismal process (56 genes), developmental process (46 genes), multi-organism process (45 genes), reproduction (38 genes), and reproductive process (38 genes). Cellular component included predominantly membrane (469 genes), membrane part (415 genes), cell (386 genes), cell part (386 genes), organelle (252 genes), organelle part (81 genes), extracellular region (47 genes), extracellular region (43 genes), membrane-enclosed lumen (18 genes), cell junction (17 genes), and symplast (16 genes). Molecular function included predominantly binding (739 genes), catalytic activity (704 genes), transporter activity (93 genes), molecular function regulator (26 genes), nucleic acid binding transcription factor activity (23 genes), molecular transducer activity (18 genes), signal transducer activity (14 genes), and structural molecule activity (11 genes) (Fig. 3, Table S7).

KEGG pathway enrichment analysis of DEGs

KEGG pathway analysis showed that 708 unigenes were annotated and assigned to 114 pathways. The 114 KEGG pathways were divided into five KEGG categories: metabolism, genetic information processing, organismal systems, environmental information processing, and cellular processes, including 401, 131, 118, 67, and 34 DEGs, respectively. In the metabolism pathway, there were additional genes involved in starch and sucrose metabolism (78 genes), phenylpropanoid biosynthesis (40 genes), and pentose and glucuronate interconversions (31 genes). In the genetic information processing pathway, additional genes were involved in spliceosome (34 genes), protein processing in endoplasmic reticulum (22 genes), and ribosome biogenesis in eukaryotes (15 genes). In the organismal systems pathway, the largest number of genes (104 genes) were involved in plant-pathogen interaction, and some (14 genes) were involved in the circadian rhythm-plant. In the environmental information processing pathway, the genes involved in plant hormone signal transduction were the highest (51 genes), followed by those involved in ABC transporters (14 genes). In the cellular processes pathway, the genes involved in endocytosis were the most prevalent (22 genes), followed by genes involved in peroxisome (eight genes) (Fig. 4, Table S8).

The top 30 KEGG pathways with the most number of annotated genes were plant-pathogen interaction, starch and sucrose metabolism, plant hormone signal transduction, phenylpropanoid biosynthesis, MAPK signaling pathway-plant, spliceosome, pentose and glucuronate interconversions, flavonoid biosynthesis, protein processing in endoplasmic reticulum, endocytosis, stilbenoid, diarylheptanoid and gingerol biosynthesis, amino sugar and nucleotide sugar metabolism, cyanoamino acid metabolism, ribosome biogenesis in eukaryotes, alpha-Linolenic acid metabolism, galactose metabolism, circadian rhythm-plant, ABC transporters, ribosome, homologous recombination, ascorbate and aldarate metabolism, RNA transport, carbon metabolism, cysteine and methionine metabolism, glycerophospholipid metabolism, biosynthesis of amino acids, ubiquitin mediated proteolysis, tyrosine metabolism, isoquinoline alkaloid biosynthesis, and monoterpenoid biosynthesis (Fig. 5, Table S9).

Clusters of orthologous genes (COG) annotation

After searching the COG database, 577 unigenes were annotated to 22 COG categories. General function prediction only had the largest number of DEGs (124 genes). The second most prevalent DEGs were secondary metabolites biosynthesis, transport, and catabolism (81 DEGs), carbohydrate transport and metabolism (72 DEGs), signal transduction mechanisms (58 DEGs), and lipid transport and metabolism (57 DEGs) (Fig. 6, Table S10).

Gene set enrichment analysis (GSEA)

The most significant four pathways of biological process in the GSEA-GO analysis were signal transduction (275 genes), DNA replication (45 genes), response to fungus (36 genes), and meristem maintenance (20 genes) (Fig. 7). The most significant four pathways of cellular component in the GSEA-GO analysis were nucleosome (43 genes), small-subunit processome (34 genes), extracellular space (27 genes), and cyclin-dependent protein kinase holoenzyme complex (14 genes) (Fig. 8). The most significant four pathways of molecular function in the GSEA-GO analysis were protein serine/threonine kinase activity (362 genes), ADP binding (308 genes), endonuclease activity (241 genes), and terpene synthase activity (14 genes) (Fig. 9, Table S11).

The most significant four pathways in the GSEA-KEGG analysis were starch and sucrose metabolism (345 genes), flavonoid biosynthesis (73 genes), propanoate metabolism (50 genes), and fatty acid biosynthesis (41 genes) (Table S12, Fig. 10).

Differentially expressed transcription factor genes

A total of 901 transcription factor (TF) genes were identified from all DEGs, which were divided into 51 TF families, including 434 downregulated DEGs and 467 upregulated DEGs (Table S13). NAC, bHLH and B3 contained the largest number of DEGs, and the gene expression pattern was shown in Figs. 11, 12 and 13. The NAC transcription factor family contained the largest number of DEGs, up to 115, of which 73 were upregulated and 42 were downregulated in flower buds with pistil abortion, respectively (Table S14). The bHLH transcription factor family contained 86 DEGs, of which 38 were upregulated and 48 were downregulated in flower buds with pistil abortion, respectively (Table S15). The B3 transcription factor family contained 65 DEGs, of which 34 were upregulated and 31 were downregulated in flower buds with pistil abortion, respectively (Table S16). Other transcription factor families contained more DEGs, MYB_Related contained 45 DEGs, WRKY contained 43 DEGs, C3H contained 40 DEGs, M-type contained 40 DEGs, HSF contained 34 DEGs, ERF contained 31 DEGs, bZIP contained 30 DEGs, C2H2 contained 29 DEGs, and MYB contained 27 DEGs (Fig. 14, Table S17).

DEGs related to starch and sucrose metabolism

In starch and sucrose metabolism, 82 related DEGs were identified, of which 58 were upregulated and 24 were downregulated, including 56 sucrose synthase, eight beta-glucosidase, three beta-fructofuranosidase, three glycogen phosphorylase, three trehalose 6-phosphate phosphatase, one alpha-amylase, one beta-amylase, one beta-carotene hydroxylase, one abscisate beta-glucosyltransferase, one endoglucanase, one xanthoxin dehydrogenase, one isoamylase, one abscisic acid 8′-hydroxylase, one sucrose-phosphate synthase. The expression patterns of these genes were shown in Fig. 15 (Table S18).

DEGs related to phytohormone biosynthesis and signaling pathways

There were 87 genes involved in the phytohormone biosynthesis and signaling pathways, of which 43 genes were upregulated and 44 genes were downregulated. There were 15 genes involved in the auxin signaling pathway, of which seven genes were upregulated and eight genes were downregulated. There were nine genes involved in cytokinin biosynthesis, of which three genes were upregulated and six genes were downregulated. There were two genes involved in the gibberellin signaling pathway, of which one gene was upregulated and one gene was downregulated. There were seven genes involved in the abscisic acid signaling pathway, of which four genes were upregulated and three genes were downregulated. There were four genes involved in the ethylene signaling pathway, of which two genes were upregulated and two genes were downregulated. There were 31 genes involved in the steroid signaling pathway, of which 15 genes were upregulated and 16 genes were downregulated. There were seven genes involved in the jasmonate signaling pathway, of which four genes were upregulated and three genes were downregulated. There were five genes involved in the salicylic acid signaling pathway, of which three genes were upregulated and two genes were downregulated. There were seven genes involved in the jasmonate acid signaling pathway, of which four genes were upregulated and three genes were downregulated (Table S19). The expression patterns of these genes were shown in Fig. 16.

Real-time quantitative PCR (qRT-PCR) validation

To confirm the accuracy and reproducibility of the transcriptome sequencing results, 16 representative genes were chosen to validate the expression levels between normal and abortive flowers by qRT-PCR. The relative expression of the selected genes was further compared with that of the transcriptome sequencing analysis. The relative trends in the expression patterns of the qRT-PCR results were all consistent with the transcriptome sequencing data, and the correlation coefficient was determined to be 0.797, supporting the reliability of the transcriptome sequencing results in this study (Fig. 17, Table S20, Table S21).