Comparative transcriptome analysis reveals key insights into male sterility in Salvia miltiorrhiza Bunge

Plant materials

This male sterile plant was accidentally discovered in the F₁ progeny of a S. miltiorrhiza (breeding lines from Shandong, China; purple flowers) crossed with S. miltiorrhiza f. alba (breeding lines from Sichuan, China; white flowers) in 2017. We temporarily named the sterile plant SW-S and the fertile plant in the F₁ offspring SW-F. We recorded the dates of the main growth and young panicle differentiation stages in SW-S and SW-F. Overall, the anther’s developmental stages can be divided into four stages based on the microsporogenesis progression: pollen mother cells, meiotic division, tetrad, and mature pollen grains. In order to ensure the reliability of the RNA-Seq results, a randomized complete block design was used with three biological replications. Anthers with the same length (at the tetrad stage, anther length was about two mm) were mixed separately from 30 plants for each biological replication. Approximately 100 mg of material was used in RNA extraction. Anthers from the SW-S and SW-F were also harvested at the four stages for paraffin section observation. Moreover, the mature pollen grains of SW-F and SW-S were collected for observation using scanning electron microscopy (SEM).

Morphological analysis

Photographs of the flowers were taken using a Canon 80D digital camera (Canon, Tokyo, Japan). For light microscopic observation, the anthers of SW-S and SW-F from different developmental stages were placed in the Carnoy’s fixative (methanol/acetic acid 3:1) and stored in a 4 °C refrigerator. Paraffin sections were created following the method described by Wang et al. (2015). In brief, the specimens were dehydrated using a graded ethanol series, infiltrated with xylene (Tiangen Biotech (Beijing) Co. Ltd., China) and embedded in paraffin (Sigma-Aldrich LLC., USA). 5 µm thick transverse sections were placed onto gelatin-coated glass slides (Sigma-Aldrich LLC., USA) and stained with hematoxylin (Sigma-Aldrich LLC., USA). The sections of anthers at different stages were observed using a DP70 high-resolution camera mounted on an Olympus BX51 microscope (Olympus, Tokyo, Japan). To observe pollen morphology, the mature pollen grains of SW-F and SW-S were collected, fixed in 2.5% glutaraldehyde. The fixed plant materials were then dehydrated using a graded ethanol series (50% twice, 60%, 70%, 80%, 90%, and 100% twice, 10 min per grade) and subjected to conventional critical-point drying using liquid carbon dioxide (CO₂) as described by Karcz (2009). The material was subsequently secured to stubs before being coated with gold (Karcz, 1996). Finally, the specimens were examined using scanning electron microscopy (JEOL JSM-6360LV) at an accelerating voltage of 20 kV.

RNA-Seq library preparation and sequencing

Total RNA from each of the six samples (three biological replications of each material) was isolated using the RNA Easy Fast Plant Tissue Kit (Tiangen Biotech (Beijing) Co. Ltd., China) according to the instructions provided by the manufacturer. For each RNA sample, the degradation and contamination were monitored using 1% RNase-free agarose gel electrophoresis. RNA integrity was assessed using the RNA Nano 6000 Assay kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Messenger RNA (mRNA) was purified from total RNA using ploy-T oligo-attached magnetic beads, before breaking into short fragments (Chen et al., 2017). First-strand complementary DNA (cDNA) was synthesized using a random hexamer primer and M-MuLV Reverse Transcriptase (Roche Pharmaceutical Ltd., Switzerland). Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and M-MuLV Reverse Transcriptase. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of the 3′ ends of the DNA fragments, the NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization (Chen et al., 2017). In order to select cDNA fragments of 250–300 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then, 3 µL of USER Enzyme (NEB, USA) was applied to size-selected, adaptor-ligated cDNA at 37 °C for 15 min, followed by 5 min at 95 °C before PCR. Then, PCR was performed using Phusion High-Fidelity DNA polymerase, Universal PCR primers, and Index (X) Primer. Next, PCR products were purified (AMPure XP system, Beckman Coulter, Beverly, USA), and their library quality was assessed on the Agilent Bioanalyzer 2100 system (Agilent Technologies Inc. USA). Lastly, Illumina sequencing was performed with the Illumina HiSeq PE150 platform by Novogene Biotech Co., Ltd. (Beijing, China).

Data processing and differentially expressed gene (DEG) identification

Raw reads in fastq format were first processed using in-house perl scripts (Gaiarsa et al., 2019). In this step, the clean reads were obtained by removing adaptor sequences, low-quality reads (base number Q_phred ≤20 in more than 50% of the whole reads), and reads containing poly-N from raw reads. The Q20, Q30, and GC content of the clean data were calculated simultaneously. All downstream analyses were based on the clean data (Wu et al., 2017). The clean reads were then mapped to the S. miltiorrhiza reference genome (http://www.ndctcm.org/shujukujieshao/2015-04-23/27.html) using Hisat 2 v2.0.5 (Kim, Langmead & Salzberg, 2015). The mapped reads of each sample were assembled by StringTie (v1.3.3b) in a reference-based approach (Pertea et al., 2015). FeatureCounts v1.5.0-p3 was used to count the read numbers mapped to each gene (Liao, Smyth & Shi, 2014). Next, the fragments per kilobase million (FPKM) of each gene were calculated according to the length of the gene and read count mapped to the gene. The differentially expressed genes (DEGs) were calculated using the DESeq2 R package (1.16.1), with q < 0.05, —log₂ fold change— ≥2 (Wang et al., 2013).

Functional analysis of DEGs

We annotated the DEGs on the basis of a set of sequential Basic Local Alignment Search Tool (BLAST) iterations, designed to find the most descriptive annotation for each sequence. The DEGs were compared with sequences in the National Center for Biotechnology Information’s (NCBI) nonredundant (Nr) protein and nucleotide (Nt) databases (http://www.ncbi.nlm.nih.gov/). Gene Ontology (GO) enrichment analysis was performed using agriGO (http://systemsbiology.cau.edu.cn/agriGOv2/index.php) (Du et al., 2010). We detected which of the DEGs were significantly enriched in GO terms using a p < 0.05 as the threshold of significance.

Quantitative RT-PCR (qRT-PCR) validation

First-strand cDNA synthesis was conducted using 5 µg of total RNA from the anthers of SW-S and SW-F at the tetrad stage with a PrimeScript RT reagent kit (Takara, Dalian, China) according to the manufacturer’s instructions. The SmActin gene (DQ243702) was employed as an internal control for qRT-PCR (Yang et al., 2010). The primers for the selected genes and SmActin gene were designed using Beacon Designer 7, and are listed in Table S1. qRT-PCR assays were performed with SYBR Green Dye (Takara, Dalian, China) using a Bio-Rad CFX96 real-time PCR platform (Bio-Rad, Hercules, CA, USA) according to the following program: 95 °C for 5 min, followed by 40 cycles at 95 °C for 5 s and 55−63 °C (annealing temperature varied according to primers) for 30 s. All assays for a particular gene were performed three times synchronously under identical conditions, and RNA transcript fold changes were calculated using the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Phenotypic characteristics of SW-S and SW-F anthers

The male sterile line SW-S was discovered in the F₁ progeny of an S. miltiorrhiza (from Shandong, China) crossed with S. miltiorrhiza f. alba (white flowers). Twenty-one individuals of 40 F₁ progenies in S. miltiorrhiza × S. miltiorrhiza f. alba showed male sterility, and the remaining 19 individuals were fertile (almost 1:1). However, crossing S. miltiorrhiza f. alba as the mother and S. miltiorrhiza as the father, produced F₁ offspring that were all fertile. So, SW-S may belong to the CMS. To confirm the stability of sterility, the roots of sterile plants were used for asexual reproduction for 2 consecutive years (2018 and 2019), and the plants with asexual reproduction were planted in the experimental field of China West Normal University, Nanchong, Sichuan, China (30°49′ north (N), 106°3′ east (E)). The development of anthers in the sterile plants was observed in April of the following year. The results showed that the anthers of all sterile plants were abnormal and could not produce normal pollen, suggesting that the sterility was stable. The comparison of agronomic traits showed that there were no significant differences in overall morphology and growth except for the floral organs of the sterile line SW-S and its fertility control SW-F. The corolla of the fertility control SW-F was blue-purple, while the flower color of the sterile line SW-S was shallow (Figs. 1A, 1B); in addition, all parts of the flower organs in SW-S were smaller than those in the fertile plants (Figs. 1E–1F). The day that SW-F flowered, the anther’s normal dehiscing was accompanied by the release of a large number of mature pollen grains. By contrast, the SW-S anthers did not dehisce or had a very small amount of pollen, which was found to be inactive (Figs. 1C, 1D).

To determine the abortive morphological features of SW-S, we collected SW-S and SW-F anthers at four development stages: pollen mother cells, meiotic division, tetrad, and mature pollen grains. In order to accurately observe the cytological structure of the sterile anthers, we prepared paraffin sections of anthers from the first three stages. As shown in Fig. 2, there were clear differences between the fertile anthers and sterile anthers. The normal anther structure involved the epidermis, endothecium, middle layer, and tapetum from outside to inside. Compared with the SW-F, the sterile line SW-S at the pollen mother cell stage had an unclear cell arrangement layer, while the tapetum cells were irregular sharped (Figs. 2A, 2E). During meiosis, the tapetum cells near the pollen sac began to detach in the sterile line (Figs. 2B, 2F). After entering the tetrad stage, all tapetum cells in the pollen sac completely fell off and gradually degraded in the sterile line. Furthermore, there was little or no callose deposited around microspore cells. The tetrad microspore was shriveled and had an abnormal morphology (Figs. 2C, 2G). Hence, we considered the tetrad period to be the critical period for pollen abortion. With the degradation of the callose, microspores were released. The microspores of normal fertile plants were fusiform or subglobose, with six germ furrows (Fig. 2D). However, the microspores of the sterile line were fused together, irregular in shape, and uneven in size, with no obvious germ furrows (Fig. 2H). The results of SEM showed that the normal pollen grains were long and spherical, while the germ furrows were distinct (Fig. 3). In addition, the exine of normal pollen had clear ornamentation, a smooth murus, and obvious oval lumina (Fig. 3). Compared with normal pollen grains, sterile plants had withered pollen, which was sticky and could not be dispersed (Fig. 3). The surface of the sterile pollen had few germ furrows, the surface ornamentation was not clear, and there were no distinct lumina (Fig. 3).

Transcriptome sequencing

To obtain key insights into male sterility for S. miltiorrhiza at the transcriptional level, we employed the Illunina HiSeq PE150 platform for transcriptome sequencing analysis using the anthers from the sterile line SW-S and the fertile line SW-F at the tetrad stage. Anthers were analyzed three times with a total of three biological replicates, and the sequencing read lengths were 150 bp. After filtering out adaptor sequences, low-quality reads, and reads containing ploy-N from raw reads, 266,722,270 clean reads were obtained, with 130,305,496 reads from the sterile line and 136,416,774 from the fertile line. The average GC content was 47.2%, and the Q20 percentage exceeded 97%. The clean reads obtained from SW-S and SW-F were matched with the S. miltiorrhiza reference sequence, where the mapping rates ranged from 88.19% to 89.65% (Table 1). The productivity and quality of the sequencing indicated that the RNA-seq data were sufficient for further analysis.

Identification of differentially expressed genes

A total of 36,534 genes were detected from SW-S and SW-F, using RNA-Seq, of which 4,607 were novel genes (unannotated genes compared to the reference genome). To make it easier to distinguish these novel genes, we added “novel” before the gene identifier (ID). To determine significant differences in gene expression levels in SW-S and SW-F, we used the following thresholds: q <0.05, and —log₂ fold change— ≥2. In total, 2,571 genes were differentially expressed in SW-S and SW-F anthers (Fig. 4). These DEGs were comprised of 939 upregulated genes and 1,632 downregulated genes in the SW-S anthers compared with SW-F anthers (Table S2). Among these DEGs, 2,152 genes could be mapped to the reference gene, whereas 419 novel genes could not be mapped to the reference gene.

GO enrichment analysis of DEGs

Go is a database established by the Gene Ontology Consortium. It aims to establish a semantic vocabulary standard for defining and describing the functions of genes and proteins applicable to various species, which can be updated with the continuous development of research (Ye et al., 2006). After enrichment analysis, the DEGs found in SW-S and SW-F were annotated according to 56 functional groups (Fig. 5). Among the biological processes, the DEGs were mainly associated with the transport, localization, and regulation of biological quality. In terms of cellular components, the DEGs were mainly related to the membrane and plasma membrane. In addition, transporter activity, transmembrane transporter activity, substrate-specific transporter activity, substrate-specific transmembrane transporter activity, active transmembrane transporter activity, and ion transmembrane transporter activity constituted the related molecular functions. There were more downregulated DEGs in SW-S than upregulated DEGs for most GO terms (Fig. 5). The genes involved in microgametogenesis and developmental maturation were all downregulated in SW-S (Table 2). Among the top 20 enriched GO terms, only cell tip growth was considered to be significantly enriched (p <0.001) (Table S3). The DEGs found in different GO terms may provide valuable insight for studying anthers development in SW-S.

Validation of key DEGs by quantitative RT-PCR (qRT-PCR)

qRT-PCR was performed on 23 DEGs, including five random DEGs (three upregulated genes and two downregulated genes), eight DEGs related to microgametogenesis, and 10 DEGs involved in developmental maturation. We compared the results obtained from qRT-PCR with those generated from the RNA-Seq analysis of the transcripts. For these DEGs, the expression trends obtained using the two methods were consistent, with a correlation coefficient of R² = 0.7892 (Fig. 6A). Although quantitative differences were found between the two methods to some extent, the trends were consistent. Differences in the level of expression of these genes were also reasonable because of the different methods employed. The results of qRT-PCR analysis further confirmed that eight genes related to microgametogenesis and 10 genes related to developmental maturation were all downregulated in SW-S (Figs. 6B, 6C). These genes may be related to male sterility for S. miltiorrhiza.

Male sterility may occur at the stage of microgametogenesis

Illumina sequencing of SW-S and SW-F produced 130,305,496 and 136,416,774 clean reads, respectively, with a total of 36,534 genes. Compared with SW-F, 2,571 DEGs were detected, and 63.5% of these genes were downregulated. GO enrichment analysis of these DEGs showed that two GO terms, microgametogenesis (GO:0055046) and developmental maturity (GO:0021700), may be related to the male sterility of S. miltiorrhiza, as genes related to these two GO terms were all downregulated. The results of qRT-PCR also confirmed that, compared with SW-F, the expression of these genes in SW-S was significantly reduced. We attempted Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis on 2,571 DEGs; however, approximately 400 DEGs were enriched in terms of KEGG pathways.

Pollen is produced in the anther during two successive developmental phases, microsporogenesis and microgametophyte development. Microsporogenesis includes a series of progressive developmental stages from sporogenous cells to haploid unicellular microspores. During microsporogenesis, the diploid sporogenous cells differentiate into microspore mother cells which undergo meiotic division to form four haploid microspores (tetrad) enclosed in the callose envelope. Microgametogenesis consists mainly of the process of releasing haploid microspore from the callose envelop and producing mature male gametophytes (Tütüncü, 2014). The results of paraffin sections showed that the meiosis of SW-S pollen mother cells was normal, and a normal tetrad could be produced. Therefore, we inferred that male sterility likely occurs at the stage of microgametogenesis. The RNA-Seq results also further confirmed that all DEGs related to microgametogenesis and developmental maturity were downregulated.

DEGs involved in microgametogenesis

There were eight DEGs associated with microgametogenesis, including glycerol-3-phosphate acyltransferase 3 (GPAT3), E3 ubiquitin-protein ligase RHF1A, phosphatidylinositol 3-kinase, phosphoribosylformylglycinamidine synthase (PFAS), MYB96, MYB78, callose synthase 5, and anther-specific protein LAT52. These genes play an important role in the development of pollen. For example, a transfer (T)-DNA insertional mutant in the glycerol-3-phosphate acyltransferase (GPAT) gene exhibits abnormal tapetum morphology and defects in the formation of the pollen wall (Zheng et al., 2003). Studies have shown that the E3 ubiquitin-protein ligase RHF1 gene plays an important role in gametogenesis. Double mutants of rhf1a/rhf2a in Arabidopsis lead to defects in the formation of male gametophytes due to interphase retardation of the mitotic cell cycle during the microspore stage of pollen development (Liu et al., 2008). Phosphatidylinositol 3-kinase is essential for vacuolar reorganization and nuclear division during pollen development in Arabidopsis (Lee et al., 2008). Phosphoribosylformylglycinamidine synthase (PFAS) is an essential enzyme in the de novo synthesis of purine. It has been reported that a missense mutation in PFAS can lead to the death of embryo in cattle (Michot et al., 2017). However, there is no systematic report of the role of PFAS in pollen development in plants. We observed that the PFAS gene was downregulated in the anthers of sterile line SW-S. Therefore, we speculate that the PFAS gene may be related to male sterility of plants. As the largest family of transcription factors in plants, MYB s play an important regulatory role in growth and development. A large number of studies have shown that R2R3-type MYB transcription factors play an important role in the regulatory pathways of anther development and pollen formation, including tapetum development, callose deposition and degradation, transportation of photosynthetic products, anther cracking, and male gametophyte formation (Higginson, Li & Parish, 2010; Xin et al., 2016). MYB96 and MYB78 are R2R3-type transcription factors; thus, the downregulation of these two genes may be related to the male sterility of S. miltiorrhiza. The anther-specific protein LAT52 gene is related to late pollen development. The RNA of LAT52 appears after the mitosis of the microspore, with its expression gradually increasing with microspore development, reaching its peak in mature pollen (Twell et al., 1989). Therefore, LAT52 may be related to pollen maturation.

Among the eight DEGs related to microgametogenesis, we believe that the callose synthase 5 gene (Cals5) is the most closely related to the male sterility of S. miltiorrhiza. Callose, primarily composed of β-1,3 glucan, separates developing pollen grains, preventing their exine walls from fusing (Nishikawa et al., 2005). Cals5 encodes a β-1,3 glucan synthase required for exine formation during microgametogenesis and for pollen viability in Arabidopsis (Dong et al., 2005). Mutations in the Cals5 gene disrupt callose formation around developing microspores, causing aberrations in exine patterning, degeneration of developing microspores, and pollen sterility (Verma & Hong, 2001). We found that the Cals5 gene was downregulated in the male sterile line SW-S. Observation of paraffin sections also revealed little or no callose around the tetraploid gametes of SW-S. Therefore, we speculate that downregulation of the Cals5 gene results in a decrease in callose around the developing pollen, leading to male sterility.

DEGs related to developmental maturity

There were 10 DEGs associated with developmental maturity, including potassium channel AKT1 (containing two homologous genes), actin 3 (ACT3), putative translationally controlled tumor protein (TCTP) (containing two homologous genes), LRR receptor-like serine/threonine-protein kinase RPK2, myosin-17-like, xyloglucan 6-xylosyltransferase 2, dynamin-related protein 1C (DRP1C), and type I inositol polyphosphate 5-phosphatase 5-like. Among these DEGs, three genes were reported to be related to pollen maturation: ACT3, RPK2, and DRP1C. In Arabidopsis, expression of the ACT3 gene in mature pollen is higher than that in other tissues (An et al., 1996). This may indicate that the ACT3 gene is related to pollen maturation. RPK2 is a key regulator of anther development in Arabidopsis. Two RPK2 T-DNA insertional mutants (rpk2-1 and rpk2-2) displayed enhanced shoot growth and male sterility due to defects in anther dehiscence and pollen maturation (Mizuno et al., 2010). The dynamin-related protein (DRP) family is essential for gametophyte development in plants. Backues et al. (2010) demonstrated that both male and female gametes require at least one functioning DRP family member (DRP1 and DRP2) to progress beyond the single nucleate stage of development. Previous studies found that drp1C-1 mutants exhibited male gametophytic lethality in Arabidopsis, displaying small and shriveled pollen grains that do not germinate (Kang, Rancour & Bednarek, 2003). There were no other genes found to be directly related to pollen development, however, other genes may indirectly participate in the regulation of pollen maturity and requires further in-depth study.