Jin668, an upland cotton (Gossypium hirsutum L.) line developed by the National Key Laboratory of Crop Genetic Improvement, Huazhong Agricultural University. We have described the transformation system of Jin668 previously (Jin et al., 2006; Li et al., 2019). The wild-type (negative control) and transgenic lines were planted in Wuhan, Hubei under normal farming practices or grown in the greenhouse during the winter in 2018 - 2020. The greenhouses were kept at a temperature of 28–35/20–28°C day/night.

We conducted a genome-wide assessment and chose sgRNA through the CRISPR-P 2.0 (http://crispr.hzau.edu.cn/cgi-bin/CRISPR2/CRISPR) (Liu et al., 2017). The process of vector construction refers to our previous report (Wang et al., 2018). The different vectors were transformed into Agrobacterium GV3101 which then was transformed into cotton Jin668. Refer to published articles for cotton transgene (Jin et al., 2006).

In order to detect the editing efficiency of transgenic lines, the targeted genomic DNA was amplified by PCR with a pair of site-specific primers at the 5’ end with common bridging sequences. Specific steps were similar to previous reports (Liu et al., 2019; Ramadan et al., 2021), and the primers used ware shown in Supplementary Table 1 . PCR products were sequenced on Illumina HiSeq platform (Illumina, USA) after recovery. Hi-TOM website (http://www.hi-tom.net/hi-tom/) was used to analyze the sequencing results. In order to detect the off-target situation in the transgene lines editing process, “sgRNAcas9_3.0.5” (Xie et al., 2014) was used to predict off-target sites, following the software default settings. The “extract_targetSeq.pl” script was used in the software package to extract the flanking sequence of the off-target site on the genome. We designed off-target site primers in batches through the “batchprimer3” website (http://batchprimer3.bioinformatics.ucdavis.edu), PCR products were sequenced on Illumina HiSeq platform (Illumina, USA). “CRISPResso2” (Clement et al., 2019) was employed for sequencing results to analyze the off-target sites. The off-target sites and sequences are shown in Supplementary Table 2 , 3 . The primers for amplification of off-target sites are shown in Supplementary Table 1 .

To detect pollen viability, the anther at 0 days post-anthesis (DPA) of WT and mutants was immersed in 2,3,5-Triphenyl tetrazolium chloride (TTC) solution (8 g TTC dissolved in 1 L phosphate buffer) according to a previous report (Min et al., 2013). After being cultured in a 37°C incubator for 30 min, the staining reaction was terminated with 2% (v/v) sulfuric acid solution. Pollen grains were placed on a microscope slide and the Zeiss (Oberkochen, Germany) Axio Scope A1 microscope was used to collect images.

In polyacrylamide gel electrophoresis (PAGE) separations, the 8% non-denaturing polyacrylamide gel (acrylamide: methylene bisacrylamide = 29:1) containing the PCR amplification products were placed in the electrophoresis chamber, and the driving force was set to 60 W. After 1 hour, the sample products were immersed in 0.2% silver nitrate solution for 10 minutes. After that, the products were washed twice with ddH₂O and then put them in the chromogenic solution (1.5% sodium hydroxide, 0.4% formaldehyde) for 5 minutes. Finally, protein band patterns could be visualized and subjected to adequate analysis.

Anthers from transgenic lines and wild-type at different developmental stages were immersed in 50% FAA (50% ethanol, 5% propionic acid, and 3.7% formaldehyde) and vacuum infiltrated for 2 h at 4°C, and placed at 4°C for 24 h to fix the tissue. For dehydration, a graded ethanol series (50, 70, 80, 95, and 100%) was used and samples were embedded in paraffin. The embedded tissues were sectioned into 10 μm sections. Anther sections were stained with toluidine blue solution (1%) and the Zeiss (Oberkochen, Germany) Axio Scope A1 microscope was used to collect images. TUNEL detection of apoptosis was performed similar to previous report (Min et al., 2013). Paraffin sections of the anthers were used for TUNEL analysis of the fragmented DNA of apoptotic cells using the DeadEnd™ Fluorometric TUNEL System (G3250, Promega). The analytical wavelengths of fluorescein and propidium iodide were 520 ± 10 nm and 640 ± 10 nm by a confocal microscope (TCS SP2; Leica), respectively.

Anthers from transgenic lines and wild-type were sampled and total RNA was extracted. The library preparations were sequenced on an Illumina Novaseq platform and 150 bp paired-end reads were generated. Raw data of fastq format were firstly processed through FastQC (Andrianov et al., 2010). Paired-end clean reads were aligned to the G. hirsutum genome using Hisat2 v2.0.5 (Kim et al., 2015). FeatureCounts v1.5.0-p3 was used to count the reads numbers mapped to each gene (Liao et al., 2014). And then fragments per kilobase of exon model per million mapped fragments of each gene was calculated based on the length of the gene and reads count mapped to this gene. Differential expression analysis of two groups was performed using the DESeq2 R package (1.16.1) (Love et al., 2014). Genes with Padj <0.05 and |log₂FoldChange| >1 were assigned as differentially expressed. The GO enrichment was performed by the R package ‘clusterProfiler’ (Yu et al., 2012).

Extraction and measurement of the endogenous IAA were as described by Miao et al. (Miao et al., 2019). Three replicates, each of 100 mg of anthers from transgenic lines and wild-type, were sampled at anther developmental stage 6, 7, 9 and 10, mixed with 750 μL of ice-cold 80% methanol containing ²H₅-IAA (OlChemlm Ltd, CAS: 76934-78-5, 10 ng ml^-1) as internal standard, and shake for 16 hours in the dark at 4°C. After centrifugation at 13,000 rpm for 5 minutes, the supernatant was dried with nitrogen, and the residue was reconstituted in 300 μL of 80% methanol. Finally, the IAA content was measured using an Agilent 4000Q-TRAP HPLC-MS system.

The Arabidopsis EXCESS MICROSPOROCYTES1 (AT5G07280, AtEMS1) controls somatic and reproductive cell development in anthers (Zhao et al., 2002). To determine whether the EMS1 gene participated in the reproductive cell development in cotton, we used AtEMS1 as a query to perform BLASTP searches and identified 11 EMS1 members in G. hirsutum. To get a better understanding of the phylogenetic relationships between EMS1, a phylogenetic tree was constructed based on these 11 G. hirsutum EMS1 and AtEMS1 protein sequences. Clearly, the EMS1s were classified into four branches, the Ghir_A08G010860 (GhEMS1_A08), Ghir_D08G010810 (GhEMS1_D08), and Ghir_A09G018830 (GhEMS1_A09) in the same branch with AtEMS1 ( Figure 1A ). GhEMS1_A08, GhEMS1_D08 and GhEMS1_A09 have leucine rich repeat N-terminal domain and leucine-rich repeat sequences ( Supplementary Figure 1 ). This result indicated that the three GhEMS1 genes may have similar functions with AtEMS1.

To further identify the biological function of the EMS1 genes involved in cotton-specific developmental processes, we have summarized the expression of EMS1 genes in different organs/tissues (including roots, leaves, anthers in different length buds) of G. hirsutum. As shown in Figure 1B , most of GhEMS1 genes have expressed in anthers, and Ghir_A08G010860 (GhEMS1_A08), Ghir_D08G010810 (GhEMS1_D08), and Ghir_A09G018830 (GhEMS1_A09) were all predominantly expressed in early-stage anthers (stage 4/5, bud length: 3~5 mm) and their expression gradually decreased with anthers development, implying that these genes may play crucial roles in identification of reproductive cell development in early stage anthers.

Due to combination of two sgRNAs with tRNA can improve the transcription and knockout efficiency (Xie et al., 2015; Wang et al., 2018). Thus, two sgRNAs targeting the same GhEMS1 genes were combined by overlap extension PCR and then ligated to the expression vector pRGEB32-GhU6.7, to produce polycistronic tRNA-gRNA genes PTG1 and PTG2 vectors ( Figure 1C ). The PTG1 contained sgRNA1 and sgRNA2, which was designed to knock out GhEMS1_A08, GhEMS1_D08, and GhEMS1_A09 ( Figure 1D ). The PTG2 contained sgRNA3 and sgRNA4, which was designed to knock out GhEMS1_A08 and GhEMS1_D08 at the same time ( Figure 1D ). The prepared vectors were transformed the cotton by Agrobacterium (GV3101) ( Supplementary Figure 2 ). We obtained 3 sterile plants (KO1~3) for PTG1, and 2 sterile plants (KO4~5) for PTG2 ( Figure 2A ). Five independent transformation plants showed different degrees of necrosis-like dark spots on the surface of anther and have different pollen viability ( Figures 2A-C ). KO1-KO3, showed obvious necrosis-like dark spots on the surface of the anthers ( Figures 2A, B ). There was no pollen in the anthers of KO1 and KO2, and only a few shriveled pollen grains in the KO3 anthers ( Figure 2C ). KO4 and KO5 showed a few anthers with dark spots, and the pollen quantity and viability were lower than the wild type (WT) but higher than those of KO1- KO3 ( Figures 2A-C ).

Furthermore, a three-month fertility assay was performed on WT and the five transgenic plants. We divided plant fertility into six levels: Grade 1 indicates that anthers with necrosis-like dark spots and all anthers without pollen; Grade 2 indicates that < 25% of the anthers have a few inactive pollen grains without dehiscence with a few anthers have dark spots; Grade 3, 4, and 5 indicate that 25%, 50%, and 75% of the anthers spread pollen, respectively; Grade 6 indicates that all anthers dehiscence and release active pollen. The phenotype of every flower was recorded, and the fertility of different plants was counted. The results showed that the fertility of the WT plants was relatively stable at grade 6, with pollen viability higher than 99.5% and anthers normal dehiscence and no dark spots; KO1- KO3 were below grade 2, with no pollen or very few inactive pollen grains, and dark spots on the anther surface ( Figures 2C, D ). KO1 was the most stable, with 100% anthers of all flowers having necrosis-like dark spots on the anther surface, and no pollen ( Figure 2D ). The fertility of KO4 and KO5 was above grade 2 and below grade 6, with a few anthers have dark spots ( Figures 2B-D ).

To check the mutation at the selected target site in KO1-KO5 lines, the sanger sequencing was performed. We found that three genes, GhEMS1_A08, GhEMS1_D08 and GhEMS1_A09, all were 100% edited in KO1-KO3 plants ( Figure 2E ), and the deletion length was in the range of 2 to 29 bp ( Supplementary Figure S3 ). In the KO4 and KO5 plants, GhEMS1_A08 and GhEMS1_D08 were successfully edited with 100% editing efficiency, and no editing in the GhEMS1_A09 ( Figure 2E and Supplementary Figure S3B ). Moreover, the expression level of GhEMS1_A09 was lower in the early stage anthers, compared with the expression of GhEMS1_A08 and GhEMS1_D08 in the same stage anthers, and the expression level of GhEMS1_A09 decreased earlier ( Figure 1B ). In all, the male fertility of KO4 and KO5 was better than that of KO1-KO3 ( Figure 2C ), indicating that the three GhEMS1 genes were essential for male fertility and played crucial roles in male fertility. In addition, 38 potential off-target genes of the two sgRNAs in KO1 (complete male sterility plant) were analyzed, and no off-target effects were found in KO1 ( Supplementary Tables 2 , 3 ), this result suggested the formation of male sterility of KO1 only caused by the mutations of GhEMS1_A08, GhEMS1_D08 and GhEMS1_A09 genes.

Whether the sterile phenotype can be stably inherited to the offspring is related to the successful application of sterile mutant to cross breeding. To test the genetic stability of the necrosis-like dark spots as a marker of sterility, we applied WT pollen to the stigma of sterile KO1 plants. Then, the target site fragments in PTG1 of GhEMS1_A08, GhEMS1_D08 and GhEMS1_A09 in WT, KO1, T1 generation (KO1×WT) were amplified by PCR, and polyacrylamide gel electrophoresis (PAGE) was used to identify the editing. The results showed that the T1 generation had the same fragment with T0 generation, but due to WT pollination, some new editing types were generated, such as in T1-3, for which two sgRNAs generated new editing types ( Supplementary Figure 3C ). The individual sgRNAs showed different efficiencies in the detection of PAGE, indicating that it is necessary to connect two sgRNAs in tandem with one vector ( Supplementary Figure 3C ). In addition, the necrosis-like dark spots on the anthers could also be observed from the T0 to T3 ( Figures 2A, B and 3A ), T3 plants (n=64) had 35.9% (23/64) complete necrosis anthers ( Figure 3A ). Gene editing types were identified by Hi-TOM (Liu et al., 2019), which revealed that the T3 generation had the same editing type as the T0 generation, such as the mutations (-1 bp, -5 bp, -6 bp) at the sgRNA1 target site ( Figures 4A, C, E, G ), and mutations (-2 bp, -20 bp, -1 bp) at the sgRNA2 target site ( Figures 4B, D, F, H ). However, due to WT pollination and the retention of Cas9, some new editing types have also been generated, such as the mutation (+1 bp) at the sgRNA1 target site ( Figure 4A ). Moreover, the pollen quantity and viability of T3 plants were analyzed, the results were consistent with the gene editing efficiency and the number of necrosis-like dark spots on the anthers in these plants, suggested the GhEMS1 mutants displayed the genetic stability of the necrosis-like dark spots as a marker of male sterility.

The necrosis-like dark spots on the mature anthers could be inherited, but when and how the necrosis-like dark spots appeared on the anther surface? To explore the formation period of necrosis-like dark spots, we obtained the stage 6 to stage 14 anthers of KO1. At stage 7, yellow spots appeared on the anthers of sterile phenotype plants, which gradually deepened at stage 10 and stage 14 to form necrosis-like dark spots eventually ( Figure 3B ). Therefore, the dark spots on the surface of the anthers started earlier, which can help us to screen sterile plants at the early stage.

To observe the microspore development of KO1, anther tissue cross-sections were made. At stage 6, the WT exhibited four complete anther cell layers (from outside to inside: epidermis, endothecium, middle layer, and tapetum layer) and microsporocytes. However, the sterile mutant anthers lacked middle layer and tapetum cells ( Figures 5A, C ). Furthermore, in the WT anthers, the microsporocytes completed nuclear division at stage 6 ( Figure 5A ), tetrads formed at stage 7 ( Figure 5A ), microspores were released from tetrads at stage 8 ( Figure 5A ), and then microspores developed into mature pollen during stage 9 to stage 11 with the tapetum development and degeneration ( Figure 5A ). However, the mutants could complete nuclear division but not cytoplasmic division, causing the enlarged microsporocytes and undetached microsporocytes at stages 7 to 9 ( Figure 5A ). At stage 10, the microsporocytes of mutants started to degrade, and they were completely degraded at stage 11 ( Figure 5A ). What’s more, the degree of DNA fragmentation in the anthers was detected by TUNEL (terminal deoxynucleotidyl transferase dUTP nick end labeling). In the WT, there were a few yellow fluorescence signals in the tapetum layer, a few microspores at stage 9, and the yellow fluorescence signals were observed enhanced at stage 10. However, no fluorescence signals were observed in the stage 9 anthers of the mutant, and only faint fluorescence signals were observed appearing in the degrading microsporocytes in the stage 10 anther locules in the mutant ( Figures 5B, C ). This result indicated that there was no normally developed pollen in the chamber of the mutant and pollen deformity was caused by the absence of tapetum formation and the PCD process of microspore mother cells.

To explore the molecular mechanisms of anther abortion in GhEMS1 mutants, we compared the transcriptomes of WT and Ghems1 anthers at four developmental stages (stages 6, 7, 9, 10). A total of 7,172 genes were found to be differentially expressed between Ghems1 and WT at four anther developmental stages ( Supplementary Table 4 ). Of these differentially expressed genes (DEGs), 2,070 (28.86%) were up-regulated and 5,102 (71.14%) were down-regulated (|log₂(fold change)|≥1 and padj < 0.05) in Ghems1 ( Figure 6A ). Compared with WT, Ghems1 had more up-regulated genes than down-regulated genes at stages 6 and 7 ( Figure 6A ). However, at stages 9 and 10, Ghems1 had more down-regulated genes than up-regulated genes ( Figure 6A ). Of these genes, 239, 1,250, 551 and 631 genes were unique at stages 6, 7, 9 and 10, respectively, and 23 genes showed differential expression in all four developmental stages ( Supplementary Table 5 and Figure 6B ). These 23 genes included three GhEMS1 genes that have been edited, and the expression of three GhEMS1 genes was downregulated in all four stage mutant anthers ( Figure 6E ). The expression level of GhEMS_D08 was higher than that of GhEMS_A08/A09, and gradually decreased following the anther development. GO enrichment analysis, “peroxidase activity” was highly enriched in stage 6 and 7 anthers of mutants, and “monooxygenase activity” is highly enriched in stage 6, 7, 9, which is related to peroxide in GO enrichment analysis ( Supplementary Figure 4 and Table 6 ). So, we detected the content of peroxide at four developmental stages of Ghems1 and WT anthers ( Figure 6G ). The results showed that compared with the WT, the mutants had lower peroxide content in stage 6 and 7, but higher peroxide content in stage 9 and 10. Black spots began to appear in stage 7 and appeared in large numbers in stage 9, might be responsible for the appearance of black spots on the surface of anthers. In mutant stage 6 down-regulated genes, “activation of MAPKK activity” and “MAP kinase kinase kinase activity” were enriched. In stage 6, 9 down-regulated genes, “pollen exine formation” was enriched, while in stage 9 and 10 DEGs was enriched in “plant- type cell wall modification”, possibly related to pollen formation ( Supplementary Figure 4 and Table 6 ).

To generate co-expression networks for all DEGs and biological samples, the weighted gene co-expression network analysis was performed. A total of 8 gene modules were identified ( Figure 6C ). EMS1 genes were distributed in the red module. We extracted the red module genes that may be related to GhEMS1 genes to do a further co-expression network analysis and showed that there are four leucine-zipper transcription factors TGACG9/10 (TGA9/10, Ghir_D01G020290, Ghir_A07G011740, Ghir_D07G011800 and Ghir_A09G010350) genes, one tetrapeptide alpha-pyrone reductase 1 (TKPR1, Ghir_D03G005390), and one indole-3-acetic acid-amido synthetase (GH3.6, Ghir_D01G006580), were associated with GhEMS1. TGA9 and TGA10 are expressed throughout early anther primordia, and mutations in TGA9 and TGA10 lead to male sterility and differential defects in abaxial (Murmu et al., 2010). In the genetic framework for control of anthers development (Wilson and Zhang, 2009), the expression trends of SPL/NZZ, EMS1, TGA9/10, DYT1 and AMS were the same, but MS1 in the stage 6 anthers of mutant was higher than WT ( Figure 6E ). Previously reported that TGA9/10 was located downstream of SPL/NZZ and upstream or in parallel with DYT1 in the genetic hierarchy that controls anther development (Murmu et al., 2010), similar to EMS1. Interestingly, we also found that there was co-expression of TGA9/10 and EMS1 ( Figure 6D ), it was suggested that TGA9/10 and EMS1 may interact with each other to regulate anther development.

In the co-expression network, we found TKPR1 which co-expressed with EMS1 ( Figure 6D ), TKPR1 involved in the biosynthesis of hydroxylated tetraketide compounds that serve as sporopollenin precursors (the main constituents of exine) was essential for pollen wall development (Tang et al., 2009; Grienenberger et al., 2010). And GO enrichment analysis shows that “pollen exine formation” was highly enriched in stage 6 and 9 anthers of mutants. To confirm whether the synthesis of sporopollenin precursors in the mutants had been affected, we measured the autofluorescence intensity of sporopollenin of Ghems1 and WT anthers at stages 6, 7, 9, 10 by microscopy with UV light illumination ( Supplementary Figure 5 ) (Ma et al., 2022). In the mutant, the autofluorescence of sporopollenin was not detected in the four stages, but in WT, the autofluorescence of sporopollenin was detected in the 9 and 10 stages. This showed that the synthesis of sporopollenin was affected in the mutant.

Phytohormones play an important role in the regulation of anther development. In the co-expression network, we found GH3.6 which co-expressed with EMS1 ( Figure 6D ), GH3.6 catalyzes the synthesis of indole-3-acetic acid (IAA)-amino acid conjugates, providing a mechanism for the plant to cope with the presence of excess auxin (Staswick et al., 2005). Compared with WT, the expression of GH3.6 in the mutant gradually decreased ( Supplementary Figure 6 ). So, we detected the dynamic changes of auxin (IAA) in anthers of Ghems1s male sterile line ( Figure 6F ). During the four anther developmental stages, the free IAA content of male sterile plants KO1 changed little with the development of anther and was at a lower level. However, the free IAA content of WT increased greatly between the stage 6 and stage 7 of anthers development, and maintained a high level after that. The lower IAA content in the stage 7, 9, 10 anthers of Ghems1s mutants may be closely related to the occurrence of male sterility.

From the above results, we found that after the editing of both GhEMS1_A08 and GhEMS1_D08, there were yellow spots on the surface of the anthers and a few fertile pollen grains in the chamber. Moreover, the results of simultaneously editing three GhEMS1 genes showed that there were necrosis-like dark spots on the surface of the anthers, which contained completely aborted pollen grains, thus the necrosis-like dark spots can serve as a marker of completely male-sterile cotton lines with GhEMS1s mutants, and can help breeders to screen sterile plants at the early anther stage.