Yongjie Zhng,Yng Hn,Meng Zhng,Xuexin Zhng,Liping Guo,Tingxing Qi,Yongqi Li,Junjun Feng,Hilin Wng,Huini Tng,Xiuqin Qio,Lingling Chen,Xitong Song,Chozhu Xing,*,Jinyong Wu,,*

a Zhengzhou Research Base,State Key Laboratory of Cotton Biology,School of Agricultural Sciences,Zhengzhou University,Zhengzhou 450001,Henan,China

b State Key Laboratory of Cotton Biology/Key Laboratory for Cotton Genetic Improvement,Ministry of Agriculture,Institute of Cotton Research of Chinese Academy of Agricultural Sciences,Anyang 455000,Henan,China

c State Key Laboratory of Cotton Biology(Hebei Base),College of Agronomy,Hebei Agricultural University,Baoding 071001,Hebei,China

Keywords:Cytoplasmic male sterility Mitochondrial genome Chimeric gene Reactive oxygen species ATP

ABSTRACT Plant cytoplasmic male sterility(CMS)is maternally inherited and often manifested as aborted pollen development,but the molecular basis of abortion remains to be identified.To facilitate an investigation of CMS in cotton,the complete sequence of cotton mitochondrial(mt)genome for CMS-D2 line ZBA was determined.The mt genome was assembled as a single circular molecule with 634,036 bp in length.A total of 194 ORFs,36 protein-coding genes,six rRNAs,and 24 tRNAs were identified.Several chimeric genes encoding hypothetical proteins with transmembrane domains were identified.Among them,a previously unknown chimeric gene,orf610a,which is composed of atp1 and a 485-bp downstream sequence of unknown nature,was identified.RT-PCR and qRT-PCR validation indicated that orf610a was expressed specifically in a sterile line.Ectopic expression of orf610a in yeast resulted in excessive accumulation of reactive oxygen species and reduction in ATP content,in addition to inhibition of cellular growth.Transgenic A.thaliana overexpressing orf610a fused with a mitochondrial targeting peptide displayed partial male sterility.Interaction between ORF610a and the nuclear-encoded protein RD22 indicated an association between ORF610a and pollen abortion.Positive feedback during transcriptional regulation between nuclear regulatory factors and the mt CMS gene may account for the male sterility of ZBA.

1.Introduction

Mitochondria(mt)are semiautonomous organelles of eukaryotic cells and are the main sites for intracellular oxidative phosphorylation and synthesis of adenosine triphosphate(ATP),providing the necessary energy for various life activities in the organism[1].There is vast diversity in the sizes of crop mt genomes.The largest and smallest mt genomes were reported as those of zucchini[2](982,833 bp)and rapeseed[3](221,853 bp),respectively.The mt genomes of plants have complex noncoding and repeated regions[4,5],and the latter may contain many chimeric genes and open reading frames(ORFs),which can lead to low fertility[6].Recombination can also occur,and leads to alteration in the location and orientation of mt genes and duplication and/or deletion of DNA sequences[7,8].

Cytoplasmic male sterility(CMS)is a ubiquitous phenomenon in crop plant and is maternally inherited.The molecular basis of CMS is largely attributed to the frequent recombination or rearrangement of intramolecular or intermolecular mt DNA that results in the formation of abnormal chimeric genes.Progress has been made in understanding the molecular mechanisms of CMS in several economically important crops including rice[9],maize[10],rapeseed[11-13],wheat[14],and pepper[15].Previous studies have generally confirmed that the occurrence of CMS is closely associated with mutations in the mt genome.Many CMS genes are derived from the rearrangement of the mt genome,including orf352 in rice CMS-WA[16],orf355 in maize CMS-S[17],orf224 in Brassica CMS-Pol[18],and orf256 in wheat CMS-AP[14].It is also commonly recognized that male-sterility traits in plant CMS systems are closely associated with mt dysfunction,as exemplified by orf79 and orfH79 in rice,where CMS is caused by the inability of mitochondria to provide the energy required for the normal development of pollen as the result of defective ATP synthase subunit genes[19].Compromised activity of ATP synthase encoded by the mt genome could be one of the key reasons underlying pollen abortion.In most cases,the ATP synthase genes have undergone mtDNA rearrangements,leading to the formation of new chimeric ORFs expressing novel peptides that are often cytotoxic[20].Male sterility in plants is usually caused by the interaction of CMS genes with nuclear-encoded mitochondrial factors.In the CMS-WA rice line,WA352,which is expressed specifically in the tapetum at the microspore mother cell(MMC)stage,interacts with the nuclear gene OsCOX11 to cause a burst of reactive oxygen species(ROS),while cytochrome c(Cyt c)oxidase is transferred from the mitochondria to the cytoplasm to the nucleus.As a result,the tapetum layer emits a programmed cell death(PCD)signal,causing male sterility[16].In maize,pollen abortion is caused by the mt CMS gene orf355,which mediates the upregulation of a nuclearencoded transcription factor,ZmDREB1.7[21].In cotton,despite its economic importance as a fiber and oil crop,it remains unknown how CMS genes induce male sterility and whether CMS induction involves interaction between CMS proteins and nuclear-encoded mt factors.

In cotton,CMS was first reported in CMS-D2,which contains a G.hirsutum nucleus and G.harknessii Brandegee cytoplasm.It is homonuclear and heterogeneous with the maintainer,but the molecular mechanism underlying CMS remains obscure.Present research on cotton CMS has focused mainly on genetic variation and marker development for breeding purposes[22],with little effort toward the identification and functional verification of CMS genes.In the present study,de novo sequencing of the CMS-D2 mt genome was performed,and by comparative analysis,the most promising candidate gene was identified for CMS in CMS-D2.Molecular mechanisms of candidate genes leading to pollen abortion were analyzed by RT-PCR,qRT-PCR,ectopic expression in yeast,biochemical analysis,A.thaliana overexpression and Y2H analysis.Hope the results and conclusions can provide a theoretical basis for the molecular abortion mechanism of cotton CMS.

2.Materials and methods

2.1.Plant materials and mt DNA preparation

The CMS-D2 line ZBA[an A line with genotype S(rf1rf1)],the maintainer line ZB[N(rf1rf1),B line],the restorer line ZBR[N(Rf1Rf1),R line]and the F1generation derived from the cross of the A and R lines were provided by the State Key Laboratory of Cotton Research Institute,Chinese Academy of Agricultural Sciences,Anyang,Henan,China[23].CMS in the cotton line ZBA with cytoplasm from G.harknessii and its isogenic maintainer ZB with G.hirsutum(AD1)cytoplasm were provided by the Institute of Cotton Research(ICR),Chinese Academy of Agricultural Science(CAAS),Anyang,Henan,China.Plants were grown in a greenhouse under standard growth conditions at ICR.Clean and healthy mature seeds of ZBA were selected and sown in a germination box following soaking in water for 24 h,and maintained in the dark at 25°C for 7 days.The roots of approximately five grams etiolated seedlings thus obtained were collected,washed,and snap-frozen in liquid nitrogen prior to storage at-80 °C.

2.2.Mt genome sequencing and genome assembly

Approximately five grams of roots of etiolated seedlings were collected and mt DNA was extracted using a modified cetyltrimethylammonium bromide(CTAB)method[24].Then one microgram purified mt DNA was fragmented using the Covaris M220 system and used to construct a 450-bp short insertion library using the NEBNext Ultra DNA Library Prep Kit for Illumina(New England Biolabs,Ipswich,MA,USA)according to the manufacturer’s instructions.Libraries were then sequenced with 150 bp paired-end reads length on the Illumina NovaSeq 6000 platform(Biozeron Co.,Ltd.,Shanghai,China)[25].In parallel,the size selected SMRTbell library(15-20 kb)was prepared following DNA size fractionation using BluePippin(Sage Science,Beverly,MA,USA)and sequenced on a PacBio Sequel II platform(Biozeron Co.,Ltd.,Shanghai,China)following the manufacturers’instructions.

Before assembly,the original Illumina reads were first filtered by removing reads with adapter sequences or Q values below 20(Q＜20),short fragments less than 75 bp in length after removing adapter and quality trimming,and those containing as many as 10% unknown(‘‘N”)and repetitive sequences.Both the PacBio Sequel data and Illumina NovaSeq data were used to reconstruct the mt genome as follows.First,the short Illumina sequencing clean data were de novo assembled using GetOrganelle 1.7.1[26].The assembled Illumina NovaSeq sequence reads were then aligned to the PacBio Sequel reads using BWA 0.7.17[27],and the valid PacBio clean data were retained.Then,the extracted Pac-Bio data were combined with Illumina data to assemble with SPAdes 3.14.1 software[28],and the candidate sequences with sequencing coverage depth more than 50×and assembly length greater than 500 bp were aligned against the plant mitogenome database(ftp://ftp.ncbi.nlm.nih.gov/refseq/release/mitochondrion/)to confirm the long mt scaffold sequences,and these sequences were then joined by overlapping.The clean reads were aligned back to the mt reference genome sequences(GenBank IDs:JX536494.1,JX944505.1,JX065074.1,and JX944506.1)and the bases were corrected using Pilon 1.23[29].Finally,the starting position and direction of the mt scaffold were determined from the cotton reference genome,based on which the final draft of the mt genome sequence of the CMS-D2 line ZBA was completed and made publicly available at the National Center for Biotechnology Information(NCBI)in the Sequence Read Archive(SRA),under accession number:SRR17406254.

2.3.Analysis and annotation of the mt genome

The genes in the assembled CMS-D2 line ZBA mt genome were annotated by both homology prediction and de novo prediction.The protein sequences were aligned with the cotton mt reference genome using GeneWise(https://www.ebi.ac.uk/Tools/psa/genewise/),and de novo gene prediction was performed with EVidenceModeler 1.1.1[30].The resulting gene sets was then manually corrected and integrated to obtain the final protein-encoding genes.Non-coding RNAs(ncRNAs)were predicted using the homology prediction methods of tRNAscan-SE 1.3.1[31]and rRNAmmer 1.2[32].The complete mt genome was then searched against the Kyoto Encyclopedia of Genes and Genomes(KEGG)[33-35],Clusters of Orthologous Groups(COG)[36,37],Non-Redundant Protein(NR),Swiss-Prot[38],and Gene Ontology(GO)[39]databases using the BLAST+2.7.1 tool with a typical cut-off E-value of 1e-5[40].A circular diagram of the mt genome was drawn with OrganellarGenomeDRAW 1.2[41].

2.4.Identification of candidat e CMS genes

Based on previous reports on CMS genes[16,20,21,42,43],sequences of chimeric genes were extracted from the mt genome and aligned with the mt genome of the G.hirsutum maintainer line 2074B(GenBank ID,JX065074.1)[1]using BLASTN(https://blast.ncbi.nlm.nih.gov/Blast.cgi)with an identity of 99% and an E-value cutoff of the 1e-5.The presence of mt chimeric ORFs and a transmembrane domain were used as the two key screening criteria.The transmembrane domains of all candidate ORFs were predicted using TMHMM Server 2.0(https://www.cbs.dtu.dk/services/TMHMM/).

2.5.RT-PCR and qRT PCR analysis of candidate ORFs

The flower buds of CMS-D2 line ZBA,the maintainer line ZB,the restorer line ZBR and the three-line hybrid F1derived from the cross of ZBA and restorer line at the 3-mm stage were used for measuring expression of candidate genes[44].Total RNA from flower buds was isolated with a Spectrum Plant Total RNA Kit(Sigma-Aldrich,St.Louis,MO,USA)according to the manufacturer’s protocol.A total of 1μg total RNA was first used for firststrand complementary DNA(cDNA)synthesis using a PrimeScript RT reagent Kit for Perfect Real Time(RR037A,Takara,Japan)according to the manufacturer’s instructions.Semiquantitative RT-PCR was performed with the following program:3 min at 94 °C,followed by 30 cycles of 30 s at 94 °C,30 s at 55 °C,and 90 s at 72°C.qRT-PCR was performed with a CFX96 Real-Time System(Bio-Rad,Hercules,CA,USA)using TransStart Top Green qPCR SuperMix(Transgen Biotech).G.hirsutum His3 was used as the internal reference gene for qRT-PCR.The relative gene expression level was calculated by the 2-ΔΔCT method[22,45-47].Table S1 lists the primers used for semiquantitative RT-PCR and qRT-PCR.

2.6.Expression of orf610a in yeast cells

Gene-specific primers were used to amplify the full-length and transmembrane/nontransmembrane domain-coding region of orf610a from the cDNAs of the ZBA by PCR.PCR products were cloned into the pYes2 vector(Invitrogen,Waltham,MA,USA)pre-digested with HindIII and BamHI.The expression plasmid thus obtained was transformed into bakers’yeast(Saccharomyces cerevisiae)strain INVSc1.For inducing expression,yeast cells were grown on supplemented synthetic complete(SC)-ura medium consisting of 0.67 g L-1yeast nitrogen base,20 g L-1glucose,and appropriate amino acids without uracil to select for positive transformants on YPGal medium consisting of 0.67 g L-1yeast nitrogen base,20 g L-1galactose,and appropriate amino acids.The transformed yeast cells in 1 mL YPGal were cultured overnight in an incubator at 29 °C and their total RNA was extracted with a Yeast RNA Extraction Kit(Nanjing Jiancheng Bioengineering Institute,Nanjing,Jiangsu,China).First-strand cDNA was then synthesized using a RNA reverse transcription kit(RR037A,Takara,Japan)and used for RT-PCR as previously described[22,45-47].Yeast actin was used as the internal control.Table S1 lists the primers used for semi-quantitative RT-PCR.

The viability of the transgenic yeast expressing orf610a was determined based on growth density.The overnight cultures of yeast grown at 30 °C with shaking(220 r min-1)were adjusted to OD600of 1.0,and tenfold serial dilutions were spotted onto SC-ura and YPGal agar plates.These plates were incubated at 30 °C for 48 h.To test the growth of yeast cells,it was necessary to culture the liquid transformed yeast overnight at 30 °C with shaking(220 r min-1)and monitor the OD600.

2.7.Detection of ROS and ATP content in yeast cells

After culture for 10 h in SC-ura or YPGal medium,the recombinant cells were collected and resuspended in Tris/HCl(pH 8.0,1 mmol L-1)buffer.The cells were incubated with 10μg of 2′,7′-dichlorodihydrofluorescein diacetate(H2DCFDA,an oxidantsensitive probe)for 2 h at 30 °C.After the yeast cells were briefly rinsed with Tris/HCl buffer for three times to remove the dye,images were collected with a Leica confocal microscope(Leica Micosystems,Wetzlar,Germany).The ATP content was determined using a luciferin-luciferase method-based ATP assay kit(Beyotime,Beijing,China).After culture for 10 h at 30 °C in SC-ura and YPGal medium,the yeast cultures were adjusted to OD600of 1.0.Then,2 mL of yeast cells were used to determine the ATP content following the manufacturer’s instructions.A Synergy HT microplate reader(BioTek,Winooski,VT,USA)was used to immediately measure the luminous intensity,and the result was calculated against the standard curve of the serial dilution of the ATP standard solution.

2.8.Vector construction and genetic transformation of A.thaliana

The atp31-231encoding a mitochondrial targeting peptide of A.thaliana(GenBank accession NM_128864.4)was fused with the ORF of orf610a and inserted into the modified pCAMBIA2300 vector.The construct was then introduced into Agrobacterium tumefaciens strain GV3101,which was used to transform A.thaliana as previously described[48].

2.9.Yeast two-hybrid(Y2H)assays

The restorer line ZBR anther cDNA library was constructed by cloning cDNAs derived from developing anthers of ZBR plants into the pGADT7 vector using the CloneMiner cDNA Library Construction Kit(Invitrogen).The vectors were transformed into S.cerevisiae strain Y187.In parallel,the ORF of orf610a(1-610 aa)was amplified and cloned into the pGBKT7 vector,and the resulting pGBKT7-orf610a construct was transfected into the S.cerevisiae strain Y2HGold to generate a bait clone.The constructed yeast library was then screened with the bait ORF610a via the yeast mating method per the instructions of the Matchmaker Gold Yeast Two-Hybrid System(Clontech,Mountain View,CA,USA).For verification,the plasmids of positive clones were extracted and sequenced.To confirm the specific protein-protein interactions and the region of interaction,ORF610a(1-610 aa,1-449 aa,399-499 aa and 480-610 aa)and RD22 were amplified by PCR and cloned into both the bait vector pGBKT7 and the prey vector pGADT7.A pGBKT7/bait plasmid and a pGADT7/prey plasmid were co-transfected into the yeast strain Y2HGold and verified on selective media DDO(SD-Trp-Leu)and DDO(SD-Trp-Leu-His-Ade)following the instruction provided for the Matchmaker Gold Yeast Two-Hybrid System(Clontech).

3.Results

3.1.Sequencing and assembly of the mt genome in ZBA

The mt genome of ZBA was sequenced by both short-read(Illumina)and long-read sequencing(PacBio Sequel II)technologies.For the Illumina sequencing,6822.2 Mb raw data and 5822.9 Mb clean data were obtained.The Q20 and Q30 values were 96.61%and 90.19%,respectively.The GC content of the sequences was 42.81%.For PacBio sequencing,the number of subreads after filtering was 33,789,the total length of subreads after filtering was 578,351,143 bp,and the lengths of subreads with N50 and N90 were 19,551 and 10,940 bp,respectively.The mean length of the subreads was 17,117 bp,although the longest one was recorded as 136,673 bp.A draft of the ZBA mt genome 634,036 bp in length and 44.9%GC content was constructed by de novo assembly(Tables S2,S3).

3.2.Structure and content of the mt genome of ZBA

Fig.1.Circular map diagram of the mt genome of the cytoplasmic male sterility(CMS)line ZBA.The details of the DNA strands transcribed clockwise(+)and counterclockwise(-)are displayed on the inside and outside of the circle,respectively.The color indicates the function of the gene,as shown in the legend at the bottom left of the figure.

Table 1Features of the mt genome assembly.

Following annotation,the ZBA mt genome sequences were assembled into a single circular molecule free of gaps(Fig.1).A total of 194 genes,24 transfer RNAs(tRNAs),and six ribosomal RNAs(rRNAs)were identified(Table 1),all of which encoded proteins of more than 100 amino acids.Gene annotation was performed using five reference databases,among which NR and Swiss-Prot accounted for the largest number of annotated genes,whereas GO and KEGG,which accounted equally for the fewest annotated genes(Fig.2A).Eleven genes were common to the five reference database annotations.The 194 mt genes could be grouped into six classes based on the number of amino acids(Fig.2B).Some 76% of the genes showed lengths in the range of 100-200 aa and fewer in the ranges of 400-500 and 600-700 aa.Thirty-six protein-coding genes were identified(Table S4),among which three were present in multiple copies(5 of nad1,5 of nad2,and 3 of nad5).It appeared that most of the genes were associated with the oxidative phosphorylation system and electron transport,encoding nine NADH-ubiquinone oxidoreductases(complex I),two succinate dehydrogenases(complex II),one cytochrome bc1(complex III),three cytochrome c oxidases(complex IV),and five ATP synthases.Genes encoding six ribosomal small subunit proteins and four ribosomal large subunit proteins were also identified.The numbers of genes encoding ribosomal proteins in the mt genome vary among plant species.Other genes identified,included the four cytosine C genes ccmB,ccmC,ccmFC,and ccmFN,a maturase gene matR and a gene mttB encoding a transporter.Nine genes in the entire mt genome contained introns,with four introns in nad5,nad1,nad2,and nad7,three in nad4,and one in rps3,rps10,and ccmFc.

Fig.2.Distribution of annotated genes in five reference databases and by amino acid size.(A)Venn diagrams showing the number of genes annotated with five reference databases.(B)Distribution of genes by amino acid number among 194 genes of the mt genome.

Fig.3.Distribution of three types of ORFs in the mt genome.The bar graph(A)and the Venn diagram(B)show the distribution of three types of ORFs.CM-ORFs,numbers of chimeric ORFs in ZBA.TM-ORFs,numbers of ORFs with transmembrane domains in ZBA;DR-ORFs,numbers of ORFs with different sequences in ZBA and the maintainer line.

3.3.Mining of CMS genes in ZBA

Many CMS genes result from rearrangements of the mt genome,displaying the following three characteristics:1.recombination with a functional gene to form chimeric ORFs;2.encoding transmembrane proteins;3.distinct from the maintainer line[42].In the ZBA mt genome,six ORFs:orf186a-2,orf142a,orf116a,orf175a,orf280a,and orf610a were identified as distinct from those in the maintainer line 2074B(Fig.3).At least one transmembrane structure was found in each of the six unique ORFs,while two or more were found in orf116a,orf175a,and orf280a(Fig.S1A-F).

Among the above six genes of interest featured chimeric ORFs,unknown sequences were located 3′end of a functional gene in orf186a-2,orf116a,orf175a,and orf610a,both 5′and 3′ends of a functional gene in orf142a,and 5′end of two functional genes in orf280a(Fig.4A).These sequences of the unknown origin were reconstituted into new chimeric ORFs with their adjacent functional genes and co-transcribed with them,with the potential to perform novel functions.Remarkably,sequence alignment analysis found that the 5′end of orf610a was the same as the 1-1348 bp of the known gene atp1,and the 3′end was a unique unknown sequence of 485 bp(Fig.S2).

For verification,RT-PCR analysis of these six chimeric ORFs was performed with the RNAs derived from the flower buds of four different genetic materials,including B,A,R,and F1.Except for orf610a,the expression levels of the other five chimeric ORF genes did show no obvious difference among the four genetic materials.Specifically,the expression level of orf610a in the sterile line was significantly higher than that in the restorer line and hybrid F1,whereas it was hardly expressed in the maintainer line(Fig.4B).

Fig.5.Expression patterns of six candidate ORFs.(A-F)qRT-PCR analysis of orf186a-2,orf142a,orf116a,orf175a,orf280a,and orf610a expression in B,A,R and F1 floral buds(3 mm).GhHis3 was used as an internal control.Error bars indicate the standard deviation(SD,n=3).B,the maintainer line;A,the sterile line;R,the restorer line;F1,three-line hybrid.

3.4.qRT-PCR analysis of the six chimeric ORFs in flower buds of four different materials

The relative gene expression analysis of six chimeric ORF genes mentioned above among the four genetic materials were also performed via qRT-PCR(Fig.5).Obviously,only orf610a was highly expressed in the sterile line but hardly discernible in maintainer line(Fig.5F),a finding consistent with the RT-PCR analysis described in section 3.3.In addition,orf610a was constitutively expressed in various organs of the ZBA(Fig.S3).

3.5.Expression of orf610a inhibits the growth of yeast

Fig.6.Effect of orf610a expression on yeast growth.(A)Growth of transformed yeast cells.Yeast cells with pYeas2 and pYesorf610a plasmids were cultured in liquid SC-ura medium at 30 °C overnight.Yeast cells with pYeas2 and pYesorf610a plasmids were cultured in SC-ura liquid expansion and YPGal expression induction media at 30 °C overnight.The overnight cultures were adjusted to an OD600 of 1.0 and then diluted 1-,101-,102-,103-,104-,and 105-fold with SC-ura and YPGal liquid media,after which each serial dilution was dripped on SC-ura and YPGal agar plates.The yeast cells were cultured in an incubator at 30 °C for 3 days.Glu,glucose;gal,galactose.(B)Growth curve of transformed yeast cells in liquid medium.The density of the transformed yeast cells was determined every hour by measuring the OD600.Cells were cultured for 13 h in SC-ura and YPGal liquid media at 30 °C.

Fig.7.ROS and ATP detection in yeast cells.(A)After confocal observation,ROS accumulation in yeast cells grown in YPGal-containing medium was detected.After culture for 6 h in liquid YPGal medium,the transformed yeast cells were dyed with H2DCFDA and observed with a confocal microscope.Scale bars,100μm.Magnification,400×.(B)The ATP concentrations of transformed yeast cells were determined by ATP bioluminescence assay.

For investigation of its function,orf610a was ectopically expressed in S.cerevisiae(Fig.S5A).In the process of inducing orf610a expression,it was found that compared with the empty vector control,the growth of orf610a transformants was somewhat inhibited.The presence of the transgene in the orf610a transformants was verified by RT-PCR(Fig.S5B).In the cell growth and viability test,the transgenic yeasts expressing either pYes2 or pYesorf610a showed similar viability on the SC-ura expansion medium plate(Fig.6A).In contrast,the pYesorf610a transformant showed reduced cell viability compared with the pYes2 transformant on the YPGal expression induction medium plate upon the induction of orf610a expression.The inhibitory effect on the growth of yeast cells as the result of expressing orf610a was also observed in liquid culture of transgenic yeast(Fig.6B),corroborating the effect observed in yeast grown on solid media.

3.6.ORF610a causes ROS accumulation and inhibits ATP production

Fig.8.Ectopic expression of orf610a causes male sterility in A.thaliana.(A)Schematic diagram of vector construction used to transform A.thaliana.(B)RT-PCR analysis of orf610a expression levels among the five transgenic lines.AtActin was used as a control.(C)qRT-PCR analysis of orf610a expression levels among the five transgenic lines.Among the five transgenic lines,line 5 showed the highest expression,lines 2,3 and 4 moderate expression,and line 1 the lowest expression.AtActin was used as an internal control.Error bars,SD(n=3).(D)Compared with the wild type,the orf610a transgenic plant showed partial abortion(white arrow).(E)Homozygous transgenic plant siliques were shortened relative to the wild type.(F)orf610a exhibited shortened siliques with fewer seeds relative to the wild type.Scale bars,200μm.(G)orf610a displayed shortened pollen tube and filament relative to the wild type.Scale bars,1 mm.(H)Pollen grains were stained with 1%potassium iodide;dark staining indicates viable pollen.Scale bars,200μm.

The fluorescence intensity of H2DCFDA as an indicator of cellular ROS showed significant increase in pYesorf610a-transformed yeast cells upon YPGal induction for 10 h,in sharp contrast to the response of the pYes2 transformants,where the H2DCFDA fluorescence was hardly discernible(Fig.7A).ATP content of the transgenic yeast expressing pYesorf610a was significantly reduced following gal induction,compared to the transgenic yeast cells transformed with either pYesorf610a or pYes2,but without gal induction(Fig.7B).Several genes associated with cellular ROS were selected for expression analysis based on our previous studies[46,49],revealing that most ROS-production-related genes were up-regulated while those involved in ROS scavenging were down-regulated in ZBA(Fig.S6).This finding supports the notion that the unique presence of orf610a in ZBA could lead to excessive ROS accumulation and ATP reduction.

3.7.Functional analysis of orf610a by overexpression in A.thaliana

To verify the association of orf610a with CMS,five stably inherited independent transgenic lines overexpressing orf610a fused with mitochondrial targeting peptide atp31-231were obtained(Fig.8A).Among them,line 5 and line 1 showed the highest and lowest expression level of orf610a,respectively,as revealed by both RT-PCR and qRT-PCR(Fig.8B,C).It was evident that the high level of orf610a expression led to compromised development phenotypes,especially in the reproductive organs,including partial pollen abortion,shortened siliques,shortened filaments and pollen tubes(Fig.8D-G).Pollen viability was also reduced in orf610aoverexpressing plants relative to wild-type controls(Fig.8H).

3.8.ORF610a interacts with the nuclear-encoded protein RD22

To investigate how orf610a regulates its target genes,yeast twohybrid(Y2H)assay was conducted to screen for the potential interactors.The full-length CDS of orf610a encoding 610 amino acids was used as a bait for the Y2H assay.A total of 167 positive clones were obtained,from which the nuclear-encoded protein RD22,which is known to be associated with dehydration and drought resistance,was selected for further analysis.The interaction between RD22 and ORF610a is illustrated in Fig.9A.To identify the region of gene interaction,DNA fragments derived from ORF610a encoding three truncated proteins of 1-449 aa,399-499 aa and 480-610 aa were selected and subjected to further Y2H assay(Fig.9B).The truncated protein consisting of 399-499 aa of ORF610a with a transmembrane domain was found to interact with RD22,whereas the other two truncated proteins without transmembrane domain did not.For verification,qRT-PCR analysis of RD22 expression in the buds of the three genetic materials was conducted,confirming that RD22 was the most highly expressed in ZBA,relative to all other genetic materials(Fig.9C).The spatial expression pattern of RD22 in ZBA was further demonstrated,showing the highest expression in flower buds and stems(Fig.9D).

4.Discussion

4.1.Overview of global characteristics of the high-quality mt genome for CMS-D2 line ZBA

Compared with the mt sequences of cotton CMS lines such as 2074A and 2074S that have been assembled into linear maps[50],here we combined Illumina NovaSeq data with PacBio Sequel data to assembly of the ZBA mitochondrial genome into a circular DNA molecule with more completeness and better assembly quality(Fig.1;Table S2).Specifically,ZBA mt genome shares high sequence similarity with other sequenced cotton mt genomes,i.e.,2074A and J4A-1,both exceeding 600 kb,and with sequence similarities higher than 99%[50,51].Comparatively,the mt genomes between pepper CMS line 138A and its maintainer line 138B were found to be highly distinct[52].This difference may be associated with the occurrence of gene recombination and rearrangement and the evolutionarily low rates of mutation in mt genomes[53].Additionally,mt genomes in the same genus have been reported to be very similar[54].Consistently,totally 194 ORFs were identified in the mt genome of ZBA,with more variations in only six chimeric ORFs relative to the sequenced mt genome of upland cotton line 2074B(Fig.3;Table 1)[50],indicating the evolutionary conservation of mt genomes in Gossypium.

Fig.9.Interaction of ORF610a with RD22 and expression pattern of RD22 gene.(A)Y2H assay of the interaction of ORF610a with RD22.The cells were diluted 1-,101-,102-,103-,and 104-fold with DDO and QDO liquid media,and then each serial dilution was dripped on DDO and QDO solid media.(B)Y2H assay mapping of the ORF610a region(transmembrane structure area)interacting with RD22.TM,predicted transmembrane segment.Aa,amino acids.Other abbreviations are as in(A).(C)qRT-PCR analysis of RD22 expression in flower buds of three different materials.(D)qRT-PCR analysis of RD22 expression in organs of ZBA.Error bars,SD(n=3).

4.2.Identification of unique CMS associated ORFs in the mt genome of ZBA and their expression patterns

CMS is a ubiquitous phenomenon in higher plants,and the associated pollen abortion is generally associated with abnormal mt ORFs[16].In most plants,a common feature of the ORFs associated with CMS is that they are recombined from essential mitochondrial genes and a sequence of unknown origin,resulting in incorrect transcription and translation and causing pollen abortion[6,16,20].In this study,orf610a was the only gene that was highly expressed in a CMS line ZBA,but was barely discernible in its isogenic maintainer line(Figs.4B,5F),in agreement with the expression patterns of CMS genes in other crops,such as orf224 in Brassica CMS-Pol[18]and orf346 in Nsa CMS[13].

The question of how CMS genes affect only male organ development remains unanswered.In the I-12 CMS line of wild beet(3)and the CMS-WA line of rice,mt genes associated with CMS were expressed throughout the plant but caused a phenotype only in the male reproductive organs[16,55].In agreement with other plant CMS studies,our study found that orf610a was expressed in various organs in ZBA,but it may function only in flower buds.In fact,most identified CMS genes were produced by mt genome rearrangement and co-transcribed with functional genes in mitochondria to form new chimeric ORFs.It is essential that these chimeric ORFs all encode transmembrane proteins[42].For example,orf352 in rice CMS-WA[16],orf355 in maize CMS-S[17],orf224 in Brassica CMS-Pol[18],and orf256 in wheat CMS-AP[14]all encoded transmembrane structural proteins.In the present study,orf610a was co-transcribed with the atp1 gene(Fig.4A).The normal atp1 gene does not encode a transmembrane protein per se,but the chimeric orf610a contains a transmembrane region that is composed of a partial atp1 sequence and an unknown downstream region produced by the rearrangement of the source sequence(Figs.S1F,S4).Thus,a transmembrane domain is an indispensable part of a plant CMS gene.

4.3.ORF610a is a cytotoxic membrane protein

Many CMS proteins have been reported as toxic and hindering the growth of yeast,including rice BT CMS ORF79[43,56],rice WA CMS WA352[16],maize S CMS ORF355[21]and rapeseed Nsa CMS ORF346[13].It has been speculated[47]that the CMS protein directly kills cells in the cytotoxic model.However,nearly all toxin-based experiments have been performed in eukaryotic systems(yeast)or prokaryotic systems(E.coli),and there is no convincing evidence that the CMS protein would cause cytotoxicity in anther cells of a higher plant.The mechanism underlying the action of CMS toxic proteins associated with male sterility remains largely obscure.Previous studies[13,16,43]suggested that the cellular toxicity of CMS related proteins to E.coli or yeast was associated with a transmembrane protein whose expression in transgenic plants led to male sterility,as exemplified by transgenic expression of a truncated transmembrane region of WA352 causing male sterility in transgenic plants.Our study is in agreement with previous study[43]indicating that the ZBA CMS-related protein ORF610a was toxic to yeast cells(Fig.6).Overexpression of orf356 in Nsa CMS resulted in male sterility in transgenic A.thaliana,but not in transgenic rapeseed[13].In our study,overexpression of orf610a in A.thaliana impeded the development of reproductive organs,leading to partial pollen abortion,shortened filaments and pollen tubes,and small siliques(Fig.8D-H).Whether the toxicity of ORF610a is the direct cause of pollen abortion in ZBA invites further investigation.

4.4.Molecular mechanism of ZBA pollen abortion

Fig.10.Schematic diagram of the molecular mechanism underlying pollen abortion in ZBA.ORF610a may cause a decrease in ATP content and a burst of ROS in pollen development and eventually lead to pollen abortion.The anterograde and retrograde regulation of orf610a and RD22 could lead to cytoplasmic male sterility of ZBA.

It is known[42]that retrograde regulation of nuclear gene expression plays an important role in the determination of male sterility in CMS plants,by triggering the expression of nuclearcytoplasmic gene modules through signals mediated by mt ROS.Nuclear genes with retrograde regulation are thus essential for mt function and pollen development.In CMS-S maize,orf355 retrogradely regulates the expression of the nuclear transcription factor ZmDREB1.7 via ROS-mediated signaling,leading to pollen abortion[21].As a member of the maize DREB1 family,ZmDREB1.7 functions mainly in dehydration and drought resistance[57].In the present study,we showed that ORF610A interacts with the nuclear-encoded RD22(Fig.9A),in agreement with a previous study[21]conducted in CMS-S maize.RD22 belongs to the BURP family,which has eight members and functions primarily in dehydration and drought tolerance,in addition to its role in the development of anthers,microspores,seeds,and other reproductive organs as well as in pectin metabolism in the cell wall[58].RD22 may act as an abscisic acid response signal molecule that participates in a variety of plant growth processes and stress responses and positively modulates ROS metabolism and potential antioxidant enzyme activities[59].In A.thaliana,overexpression of RD22 and AtDREB1A led to severe growth retardation and a reduction in fruit setting[60].We found that the interaction between ORF610a and RD22 depends on the transmembrane domain(Fig.9B),in agreement with the finding of a previous study[16]of rice CMS-WA lines,in which the interaction between WA352 and OsCOX11 occurred within a specific region.Dynamic variations in mt gene expression and metabolic activity play an important role in nucleus-mitochondria communication.We have also shown that the nuclear-encoded CMS regulatory gene RD22 was highly expressed in the flower buds of sterile lines,and its coexpression with orf610a is a prerequisite for CMS as the result of cytoplasm-nuclear interaction in ZBA anthers.It is conceivable that the expression of orf610a activates mt retrograde signaling,which in turn induces RD22 expression.Whether the reverse regulation of orf610a and RD22 is the factor causing pollen abortion in ZBA awaits further study.

There is growing evidence[13,16,21,43]that PCD and ROS are involved in the CMS pathway,and the structural characteristics of CMS proteins tend to function in two ways:toxic proteins and energy deficits.In this research,the expression of orf610a in yeast showed that it could hinder the growth of yeast cells and give rise to surplus ROS and reduction in ATP content in transformed cells(Fig.7).These findings support an inference that energy deficiency and ROS accumulation are essential for pollen abortion.The present study,together with previous studies[47,61],lead us to propose a working model for the molecular modulation of CMS in cotton.In this model,the retrograde regulation of ORF610a interacting with RD22 disrupts the homeostasis of ROS production and ATP synthesis in mitochondria,triggering a burst of ROS that is transferred to the nucleus to initiate premature PCD,leading ultimately to pollen abortion(Fig.10).

5.Conclusions

The ZBA mt genome contained 194 ORFs,36 protein-coding genes,six rRNAs,and 24 tRNAs.In the ZBA mt genome,six chimeric genes encoding transmembrane structures were identified.The expression of one of these,orf610a,caused ROS accumulation and ATP reduction,suggesting an association between the newly identified mt CMS gene orf610a and pollen abortion in ZBA.The interaction of ORF610a with the nuclear-encoded protein RD22 and their co-expression patterns may determine cotton ZBA male sterility.Identification of this module provides new insights to understand the signal transduction between the nucleus and cytoplasm underpinning cotton CMS.

CRediT authorship contribution statement

Yongjie Zhang:Formal analysis,Investigation,Writing-original draft.Yang Han:Formal analysis,Investigation.Meng Zhang:Writing-review&editing,Investigation.Xuexian Zhang:Investigation,Resources.Liping Guo:Resources.Tingxiang Qi:Investigation.Yongqi Li:Investigation.Juanjuan Feng:Investigation.Hailin Wang:Resources.Huini Tang:Resources.Xiuqin Qiao:Resources.Liangliang Chen:Investigation.Xiatong Song:Investigation.Chaozhu Xing:Conceptualization,Project administration,Writing-review & editing.Jianyong Wu:Conceptualization,Project administration,Funding acquisition,Writing-review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by funds from the National Natural Science Foundation of China(31871679),the Tianshan Youth Program(2018Q010),and the Central Public-interest Scientific Institution Basal Research Fund(1610162021015).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.02.008.