Comparative transcriptome analysis of isonuclear-alloplasmic lines unmask key transcription factor genes and metabolic pathways involved in sterility of maize CMS-C

Plant materials

The maize CMS-C line C48-2 (C means the C cytotype) was obtained by successive backcrossing using N48-2 (N means the normal cytotype) as the male donor for more than eight generations (Fig. S1A). The two lines were cultivated in the same experimental plot in Sichuan Agriculture University (Chengdu, Sichuan Province, China). Three biological replicates of each line were planted at every two days, and each replicate was planted in three-row plots with 42 individual plants. Maize inflorescences are composed of spikelets which contain two florets, an upper and a lower one. The anthers in the upper floret develop 2–3 days ahead of the 3 anthers in the lower floret (Hsu, Huang & Peterson, 1988). In order to reconcile this developmental discrepancy, the lower florets were discarded and the rest of the spikelet were used for sample preparation.

Samples preparation for qRT-PCR of transcriptional factors

Upper florets at different development stages were crushed by dissecting needles and dyed with carbol fuchsin stains, then observed by optical microscope (Leica DM1000) (Fig. S1C). Spikelets of ten independent plants at different microspore developmental stages, i.e., pollen mother cell stage (PS), dyad stage (DS), tetrad stage (TS) and mononuclear stage (MS), were collected and pooled respectively, and snap-frozen in liquid nitrogen and kept at −80 °C for RNA extraction.

Samples preparation for RNA-seq

The spikelets at PS stage and MS stage were sampled and pooled from ten individual tassels of each plot respectively, and snap-frozen in liquid nitrogen and kept at −80 °C for RNA extraction. Three biological replicate samples were collected respectively and used for RNA-seq.

RNA extraction

Total RNA was isolated using TRizol kit (invitrogen, USA) and purified using mRNA purification kit (Promega, USA) following the manufacturer’s protocol. An Ultra-microspectrophotometer NanoDrop 2000 (Thermo Fisher Scientific, USA) was used to detect total RNA concentration and purity.

cDNA library construction and RNA-seq

The mRNA-enriched RNAs were broken into short pieces with the fragmentation buffer. The short pieces and random hexamer primers were used for the synthesis of the first strands of cDNAs. The second-strand cDNAs were synthesized using the Second Strand Synthesis Reaction Buffer and the Second Strand Synthesis Enzyme Mix from NEBNext^® Ultra™ RNA Library Prep Kit for Illumina (NEB, USA). Short fragments were purified with AMPure XP Beads (Agencourt, USA) for end repair and tailing A. A-tailed fragments were mixed with NEBNext Adaptor and Blunt/TA Ligase Master Mix for adaptor ligation. Finally, the sequencing libraries were constructed following PCR amplification with these adaptor-ligated fragments. Agilent 2100 Bioanalyzer and ABI StepOnePlus™ Real-Time PCR System were used for the quantification of the cDNA library. Paired-end sequencing was performed on the Illumina HiSeq 4000 sequencing platform. RNA-seq were completed by Beijing Annoroad Gene Technology Co., LTD (China).

Read processing and mapping

All sequencing data were first processed by trimming adapter sequences, removing higher N rate (>10%) sequences and filtering low-quality reads (the ratio of the bases with quality value Q ≤ 5 greater than 50%) using the SOAPnuke software (http://soap.genomics.org.cn/) developed by BGIA. The remaining clean reads were mapped to maize B73 reference genome (RefGen_v3) and genes using BWA (version 0.7.10-r789) (Li & Durbin, 2009) and Bowtie2 (version 2.2.5) (Langmead & Salzberg, 2012), respectively.

Gene expression analysis and identification of DEGs

The Bowtie2 aligner was used to map the clean reads on the genes and all the uniquely mapped reads were transformed into FPKM (fragments per kilobase of exon model per million mapped reads) using RSEM tool to estimate transcript expression levels in all samples (Li & Dewey, 2011). DESeq2 package (version 1.14.1) was applied to detect DEGs between each chosen sample pairs (Love, Huber & Anders, 2014). The results of all statistical tests were corrected for multiple testing using the Benjamini–Hochberg False Discovery Rate (FDR ≤ 0.001) (Benjamini & Yekutieli, 2001) and |Log₂FC| ≥ 1 (FC = fold change).

GO annotation and KEGG enrichment pathway analysis

Gene Ontology (GO) and functional enrichment analysis were conducted on all identified DEGs as described (Ashburner et al., 2000). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis utilizes hypergeometric test to search DEGs that are significantly enriched in certain KEGG pathways compared with other genes in the genome. GO::TermFinder (Module Version: 0.86, http://search.cpan.org/dist/GO-TermFinder/lib/GO/TermFinder.pm) was used for accessing GO information and finding signicicantly enriched GO terms (Boyle et al., 2004). KOBAS 2.0 (Xie et al., 2011) and KEGG database (http://www.genome.jp/kegg/genes.html) were used for KEGG analysis. GO categories with corrected P value ≤0.05 and KEGG pathways with Q value (P value after multiple hypothesis testing correction with a range between 0 and 1) ≤0.05 were defined as significantly enriched.

qRT-PCR analysis

cDNA was synthesized using PrimeScript RT Master Mix Kit (Takara, Japan) following the protocols. Gene-specific primers were designed with Primer-BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The primer sequence of DEGs and the control gene were listed in Table S1. qRT-PCR were carried out using an CFX96 Real-Time PCR Detection System (Bio-Rad, USA). The mixed solution of each qRT-PCR reaction contained 1.0μL of cDNA template, 1.0μL of 10μmol/L forward and reverse primers each, and 5μL of SYBR Premix Ex Taq II (Takara, Japan) in a final volume of 10μL. PCRs were implemented as follows: 95 °C for 2 min, followed by 40 cycles of 95 °C for 10 s, 60 °C for 30 s. qRT-PCR was performed with three replicates for each sample. Quantification cycle (Cq) was used for the determination of relative expresssion levels normalized to reference genes. 18S was used as a reference gene for qRT-PCR of transcriptional factors, and Actin was used as a reference gene for verification of RNA-seq results.

Phenotypic characterization of fertile and sterile spikelets

Because C48-2 and N48-2 are a pair of isonuclear-alloplasmic lines, no phenotypic variation were founded between the two lines (Fig. S1A). No emergence of anthers and lacking mature pollen grain lead to the complete male sterility of C48-2 line (Fig. S1B). With the development of microspores, no obvious variations were detected between the microspores of C48-2 and N48-2 plants until the MS stage (Fig. S1C). It is consistent with the previous cytological study results, i.e., the premature microspore abortion and tapetum degradation at MS stage is one of the prominent characteristics of CMS-C in maize (Chen & Duan, 1988). Moraes, Bered & Kaltchuk-Santos (2008) suggested spikelet lengths are reliable parameters for predicting the specific microspore developmental stage of maize under the same planting condition. As shown in Fig. S1D, the microspores of the upper florets will be at the PS stages and MS stages when the spikelet lengths are in the 5.9 ∼ 7.2 mm range and 9.2 ∼ 10.6 mm range, respectively. For more accurate microspore development stages, the spikelets at the lengths of 6.0 ∼ 7.0 mm and 9.5 ∼ 10.5 mm were sampled and used for the RNA-seq.

Production of ready-to-analyze gene expression data

With the removal of low quality reads, a total of 326.75  ±  6.11, 322.23  ±  1.85, 327.95  ±  6.94 and 323.83  ±  4.02 million clean reads and 408.43  ±  7.64, 402.79  ±  2.31, 409.94  ±  8.68 and 404.79  ±  503 MB nucleotides were acquired for C48-2 at the pollen mother cell stage (PS-C), N48-2 at the pollen mother cell stage (PS-N), C48-2 at the mononuclear stage (MS-C) and N48-2 at the mononuclear stage (MS-N), respectively. A high proportion of the clean reads, 67.45% from PS-C, 68.69% from PS-N, 69.04% from MS-C and 69.99% from MS-N, were mapped to the maize reference genome sequence, and 73.94%∼75.82% of the clean reads were assigned to known genes via BLAST analysis (Table 1). Pearson correlation coefficient analysis revealed high correlations between biological replicates, indicating that the sequencing data were reliable (Fig. S2).

Identification of DEGs

Gene expression of PS-C, PS-N, MS-C and MS-N were quantified and compared (Fig. 1A). At PS stage, 31,892 genes were detected in the two lines, 1,527 genes and 1,677 were uniquely expressed in C48-2 and N48-2, respectively. At MS stage, there were 1,920 genes specifically expressed in C48-2 and 1,747 in N48-2, and 32,063 genes were common expressed in the two lines. For the comparison between the two stages of the same line, C48-2 had 1,646 genes and 2,210 genes uniquely expressed at PS stage and MS stage, respectively, and 31,773 genes were commonly expressed at the two stages. As to N48-2, 1,805 genes were only expressed at PS stage and the number is 2,046 at MS stage, while 31,764 genes were found to be expressed at the two stages.

FDR ≤ 0.001 and ∣Log₂FC∣ ≥ 1 were used as the threshold to screen the DEGs between the two lines at the same stages or between the different stages of the same lines. Compared with N48-2, 983 DEGs including 454 up-regulated and 529 down-regulated were identified at MS stage of C48-2, and only 11 DEGs, 1 up-regulated and 10 down-regulated, were found in C48-2 at PS stage (Fig. 1B). A large number of DEGs that is 3,553 and 2,069 showed significant expression difference between the two stages of the C48-2 line and N48-2 line, respectively (Fig. 1B).

qRT-PCR verification on RNA-seq results and DEGs identification

qRT-PCR of 16 genes in C48-2 and N48-2 at the MS stage were applied to verify the filtering of DEGs using the same samples that used in RNA-seq. The qRT-PCR expression patterns of 14 genes were consistent with the RNA-seq results, and the coincidence rate was 87.50% (Fig. 2). As to the accuracy of filtering of DEGs, all the 9 genes except LOC103638691 displayed as DEGs by qRT-PCR were identified by RNA-seq analyses, and the coincidence rate was 88.89% (Fig. 2). These results indicated the RNA-seq data in the present study were reliable.

DEGs expression patterns between different comparisons

Further, the normalized reads number of DEGs between C48-2 and N48-2 were 57.30  ±  89.82 for PS-C, 161.44  ±  197.18 for PS-N, 194.85  ±  817.49 for MS-C and 608.83  ±  3136.55 for MS-N (Table 2). The mean values for N48-2 were 2.82 and 3.12 times larger than that for C48-2 at PS stage and MS stage, respectively. It suggested that more gene expressions were down-regulated along with the substitution of C cytotype for N cytotype. Meanwhile, the normalized read numbers of DEGs between PS stage and MS stage were 161.09  ±  534.04 for PS-N, 495.37  ±  2458.30 for MS-N, 155.20  ±  406.47 for PS-C and 229.25  ±  940.67 for MS-C (Table 2). The mean value of normalized read numbers of DEGs in MS stage was 3.08 times larger than that in PS stage in N48-2. But the mean value was only increased about 48% from PS stage to MS stage in C48-2. This suggested that, along with spikelet development, more genes need to be activated in N48-2 and the activation of these genes expression was inhibited in C48-2. Collectively, in combination with the large standard deviations of the normalized read numbers, the results not only suggest reasonable coverage for the DEGs, but it’s also a reasonable indication that the inhibition of gene expression by C cytotype causes the microspore abortive at the whole genome level.

Key DEGs selection

Because the C48-2 and N48-2 were isonuclear-alloplasmic lines, only 11 DEGs were detected between the two lines at PS stage, but the number of DEGs increased greatly to 983 at MS stage (Fig. 1B). The number of differentially expressed genes was consistent with the observed phenotypic differences between N cytotype and C cytotype, anther development at the early and late developmental stages. These results demonstrate the important role that mitochondrial retrograde regulation (MRR) at MS stage plays in the maize anther development, and also made it possible to draw the conclusion that there are no obvious changes at PS stage between the transcriptome levels of C48-2 and N48-2, and focusing the DEGs at MS stage is more reasonable.

In N48-2 line, 2,069 DEGs between PS stage and MS stage were detected (Fig. 1B). It suggests that the 2,069 DEGs should be related to the spikelet normal developments from PS to MS stage. Among the 2,069 DEGs, 453 DEGs were found to be significantly differentially expressed between C48-2 and N48-2 at the MS stage (Fig. 1C and Table S2). The 453 genes were thought to be the key DEGs that participated in the process or the causes of microspore abortion. Among the 453 DEGs, 385 (84.99%) genes were down-regulated and only 68 (15.01%) genes were up-regulated in C48-2 at MS stage (Table S2). It is worth noting that 92 (23.90%) of the 385 down-regulated genes and only 1 (1.47%) of the 68 up-regulated genes were differentially expressed more than ten times (Table S2). The utilization of the 453 DEGs for the following analyses will provide a narrow range for the identification of potential pollen sterility related genes and/or metabolic pathways.

GO and KEGG enrichment analyses of the 453 DEGs

GO and KEGG analysis are widely used in the description of gene product functions (Ashburner et al., 2000; Kanehisa et al., 2004). As shown in Fig. S3, the GO enrichment of the 453 DEGs was analyzed. Six GO categories, i.e., “polyamine metabolic process” of biological process (BP), “Cytoplasmic part”, “Cytoplasm”, “Cytoplasmic vesicle” and “Vesicle” of cellular component (CC) and “Aspartate oxidase activity” of molecular function (MF) were significantly annotated (Table 3). The 453 DEGs were enriched in 78 pathways by KEGG analysis and the top 20 pathways enrichment was stated in Fig. S4. Seven pathways including “Nitrogen metabolism”, “Alanine, aspartate and glutamate metabolism”, “Fatty acid elongation”, “Plant-pathogen interaction”, “Cutin, suberine and wax biosynthesis”, “Biosynthesis of unsaturated fatty acids” and “Arginine and proline metabolism” were significantly enriched (Table 4).

The analyses of DEGs encoding transcription factors

As shown in Table 5, 18 of the 453 DEGs were identified as transcription factor genes. The expression levels of all the 18 transcription factors have no significant difference between C48-2 and N48-2 at PS stage, but all of them were down-regulated in C48-2 at MS stage except LOC100279991 and LOC100383937. The mRNA amount of five transcription factors, including four genes LOC100272446, LOC100216895, LOC100191533 and LOC100274206 encoding MYB protein, and one gene LOC732739 encoding LBD protein in C48-2 were estimated to be only 0.28% ∼17.63% of that in N48-2 at MS stage (Table 5).

With the development of microspores in N48-2, the five genes displayed similar expression patterns, i.e., extreme low or even no expression at PS and DS stages then increased at the TS stage and reached a much higher level at the MS stage, but the expression levels were significantly lower than that in N48-2 at both TS and MS stages of C48-2, especially for the expression of the four MYB genes that were inhibited extremely at the two stages (Fig. 3). Phylogenetic analysis shows that LOC100272446 and LOC100274206 were homologous to MYB33 (at5g06100) and MYB65 (at3g11440); LOC100216895 and LOC100191533 were homologous to MYB42 (at4g12350), MYB85 (at4g22680), MYB20 (at1g66230) and MYB43 (at5g16600) (Fig. 4).

Different stages of microspore development have different MRR performance in the CMS line

CMS traits were generated by the changes in mitochondrial function due to the rearrangement of the mitochondrial genome (Hu et al., 2014). The changes trigger off the alterations of the nuclear gene expression through MRR (Rhoads & Subbaiah, 2007; Fujii, Komatsu & Toriyama, 2007). In the present study, the comparative transcriptome analysis between the CMS-C line C48-2 and its isonuclear-alloplasmic maintainer line N48-2 were analyzed. No obvious variations were detected between the microspore of C48-2 and N48-2 plants until the MS stage (Fig. S1C). Meanwhile, the differential expression analyses only detected 11 DEGs at PS stage and the number increased greatly to 983 at the MS stage between the isonuclear-alloplasmic lines (Fig. 1B). As we know, the heritable alterations of phenotype result from the changes of gene expression. The huge difference of DEG numbers between the PS and MS stages not only indicated the rationality and reliability of the RNA-seq data, but also showed the space–time specificity of the MRR in maize CMS-C.

Differentially expressed patterns of most DEGs that involved in the causes or process of microspore abortion were characterized as down-regulation

The cross-comparison of the DEGs from the two comparisons of “MS-N vs PS-N” and “MS-C vs MS-N” detected 453 DEGs that provided a narrow range for the identification of potential pollen sterility related genes and/or metabolic pathways (Fig. 1C). At MS stage, 84.99% of the 453 DEGs were down-regulated in C48-2 line (Table S2). Among the 453 DEGs, 16 out of 18 genes characterized as transcription factor genes were down-regulated in C48-2 at MS stage (Table 5). GO and KEGG enrichment analyses of the 453 DEGs found that the DEGs involved in “Polyamine metabolic process”, “Aspartate oxidase activity”, “Fatty acid elongation”, “Cutin, suberine and wax biosynthesis”, “Biosynthesis of unsaturated fatty acids” and “Arginine and proline metabolism” were all significantly down-regulated in C48-2 at the MS stage (Tables 3 and 4).

The significantly down-regulated expression of MYB was essential for male sterility

The alteration of transcription factors expression often leads to the abnormal pollen development and the male sterility during plant growth. Millar & Gubler (2005) found that the myb33 myb65 had a premeiotic abortion of pollen development because the tapetum undergoes hypertrophy at the pollen mother cell stage in Arabidopsis. A battery of SND1-regulated transcription factors, including MYB20 (at1g66230), MYB43 (at5g16600), MYB42 (at4g12350) and MYB85 (at4g22680), is required for normal secondary cell wall biosynthesis and secondary cell wall thickening of the endothecium is necessary for the release of pollen grains (Zhong et al., 2008; Keijzer, 1987; Millar & Gubler, 2005; Dawson et al., 1999; Steiner-Lange et al., 2003). In the present study, four genes encoding MYB transcription factors were homologous to Arabidopsis MYB33, MYB65, MYB42, MYB85, MYB20 and MYB43 and found to be absolutely suppressed or significantly down-regulated in C48-2 along with the microspore development (Figs. 3 and 4). We suggest that the dramatic decreases of these MYB mRNA amount at the TS and MS stages might be involved in the male sterility of C48-2 line.

Down-regulation of DEGs involved in energy metabolism potentially related to CMS in maize

High respiration rate and great energy demand are usually observed during microspore development (Tadege & Kuhlemeier, 1997). ATP, the energy currency for the Cell, in male sterile lines was significantly decreased (Bergman et al., 2000; Teixeira, Knorpp & Glimelius, 2005; Li et al., 2015). NAD is essential for energy metabolism and electron transfer, de novo biosynthesis of NAD starts with the oxidation of L-Asp to α-iminosuccinate and catalyzed by L-Asp oxidase (AO) in plants (Katoh et al., 2006). In this study, two DEGs, LOC103644300 and LOC100279867, encoding AO protein were found involved in “Aspartate oxidase activity” and the mRNA amount of LOC103644300 and LOC100279867 in C48-2 were estimated to be only 4.60% and 5.24% of that in N48-2 (Table 3). We speculated that the suppressed expression of the genes that related to the energy supply might result in a shortage in energy required for microspore development and ultimately led to the male sterility of C48-2.

Down-regulation of DEGs involved in polyamine metabolic pathway might resulted in the accumulation of ROS in C48-2 line

For the male fertility of higher plant, the importance of polyamines, i.e., putrescine, spermidine and spermine, has been widely demonstrated (Falasca et al., 2010; Aloisi et al., 2016; Song & Tachibana, 2007). In plants, the ornithine decarboxylase (ODC), arginine decarboxylase (ADC) and spermidine synthase (SPDS) were essential for the synthesis of putrescine (Aloisi et al., 2016). ADC and ODC enzymes were expressed at the mononuclear stage of Nicotiana tabacum microspore and inhibited activity of SPDS led to abnormal pollen grains in kiwifruit (Bokvaj, Hafidh & Honys, 2015; Falasca et al., 2010). 5 DEGs (LOC100501307 and LOC103633099 encoding ODC protein, LOC100501649 and LOC103638134 encoding ADC protein, and LOC100279222 encoding SPDS protein) that significantly enriched in “polyamine metabolic process” were all down-regulated in C48-2 at MS stage (Table 3). Particularly, the mRNA amount of LOC100501307 and LOC103633099 in C48-2 were estimated to be only 4.30% and 5.87% of that in N48-2 (Table 3). Therefore, we hypothesized these genes’ significant down-regulation in C48-2 might lead to the reduction of cellular polyamine levels. Polyamines can protect cellular constituents from reactive oxygen species (ROS) because of the powerful properties for hydroxyl radical scavenging and singlet oxygen quenching (Das & Misra, 2004). ROS induce programmed cell death (PCD) and high levels of ROS indiscriminately attack cellular constituents then leads to male sterility (Xie et al., 2014; Song et al., 2016). Our previous study found that more ROS were produced in C48-2 than that in N48-2 at MS stage (Huang et al., 2012). Collectively, it is reasonable to draw a conclusion that the sharp down-regulated expression of the polyamine metabolism related genes might result in the reduction of cellular polyamine levels and the balance between ROS production and scavenging were disturbed, then the higher concentrations of ROS were accumulated in C48-2 line at MS stage, eventually led to the male sterility.

Down-regulation of DEGs involved in fatty acid synthesis might disrupt the pollen wall development

In most cases, the male sterility traits were on account of the defective pollen wall (Wu et al., 2015). The pollen wall consists of the exine and the intine (Dobritsa et al., 2011). The exine was derived from the tapetum and consists of sporopollenin (Jiang, Zhang & Cao, 2013). Sporopollenin biosynthesis is closely associated with fatty acid metabolism for wax and cutin formation (Morant et al., 2007; Dobritsa et al., 2009; Li et al., 2010). Here, six DEGs were significantly enriched in the fatty acid metabolism related “biosynthesis of unsaturated fatty acids” and “fatty acid elongation” (Table 4). The synthesized C16 and C18 fatty acids are transported to the endoplasmic reticulum and lengthened into very long chain fatty acids (VLCFAs) that are used for sporopollenin and/or wax biosynthesis (Jiang, Zhang & Cao, 2013). In the process of fatty acid elongation, the four enzymatic reactions catalyzed by β-ketoacyl-CoA synthase (KCS), β-ketoacyl-CoA reductases (KCR), β-hydroxyacyl -CoA dehydrase and trans-2-enoyl-CoA reductases (ECR) are the core biosynthetic activities (Costaglioli et al., 2005). In this study, LOC100286157 that encoding KCS protein, LOC103653730 (GL8B) and LOC542359 (GL8) that encoding KCR protein, LOC100274548 (PAS2) that encoding β-hydroxyacyl -CoA dehydrase, and LOC100191179 that encoding ECR protein were all down-regulated in C48-2 at MS stage (Table 4). In maize, studies had demonstrated GL8 and PAS2 are two essential enzymes in VLCFAs synthesis (Xu et al., 1997; Xu et al., 2002; Bach et al., 2008). Meanwhile, all the DEGs significantly enriched in “Cutin, suberine and wax biosynthesis” were down-regulated in C48-2 at MS stage (Table 4). Based on these results, it is also reasonable to draw another conclusion that the down-regulation of the genes involved in the fatty acid synthesis might disrupt the pollen wall development and then contribute to the pollen abortion.

Analysis of DEGs involved in proline metabolism of pollen development

In a number of plant species, massive proline accumulation in anthers and pollen had been reported by different authors (Mattioli et al., 2012; Biancucci et al., 2015). It is not very clear, to date, the reason for such a massive proline accumulation in pollen. In Arabidopsis, aberrant and infertile pollen grains were found in the proline-deficient mutant p5cs1 p5cs2/P5CS2 that defective in proline biosynthetic enzymes, and exogenous proline supplied can partially complement both morphological and functional defects of pollen (Mattioli et al., 2012; Funck et al., 2012). This supports the idea that proline is required for pollen development at the genetic and molecular level. In maize, Zou et al. (2009) found that along with the maize microspore development, the proline content was gradually increase in the maintainer line, but the proline contents were significant reduced in the CMS-C line. In this study, seven DEGs were significantly enriched in “arginine and proline metabolism” pathway and all were down-regulated in C48-2 line compared with that of N48-2 line (Table 4). These results showed that the down-regulation of these genes expression might be correlated with the proline content reduction and then made the defects in pollen development and male fertility in C48-2 line.