The four completed genome assemblies of G. arboreum (A2, CRI_V1.0) (Du et al., 2018), G. raimondii (D5, JGI v2_a2.1) (Paterson et al., 2012), hirsutum acc. TM-1 (AD1, ZJU) (Hu et al., 2019), and G. barbadense acc.H7124 (AD2, ZJU) (Hu et al., 2019) were downloaded from the CottonGen database¹ (Yu et al., 2014). The protein sequences of AtDA2 and OsGW2 were acquired from the Arabidopsis Information Resource² and the Rice Genome Annotation Project³, respectively. The protein sequences of AtDA2 and OsGW2 were used as queries to identify GW2-like genes against four genome databases of Gossypium with e-values of 1e^––10. The candidate sequences were submitted to the NCBI batch CDD program⁴ and SMART database⁵ (Letunic et al., 2015) for further confirmation.

The chromosomal distribution and the loci of all GW2-like genes were determined based on the results of identification genome annotation data. Molecular markers of reported QTLs were download from CottonQTLdb⁶ (Said et al., 2013, 2015a,b). Subsequently, the location images of GW2-like genes were drawn using the MapChart software (Voorrips, 2002).

The phylogenetic tree of GW2-like genes was analyzed and generated using multiple sequence alignments of four cultivated cotton species using ClusterW⁷ (Larkin et al., 2007). A phylogenetic tree was constructed using MEGA 7.0⁸ (Kumar et al., 2016) with pairwise distance and the neighbor-joining (NJ) method.

Gene duplication events and genomic synteny were analyzed using MCScanX and visualized using the TBtools software (Wang Y. et al., 2012; Cheng et al., 2020). The synonymous (Ks) and non-synonymous (Ka) substitution ratios were calculated using TBtools software.

The structural information of GFF3 files of four cotton species were download from CottonGen database (see text footnote 1). The protein length, molecular weight, and isoelectric point of the encoded proteins were predicted using the ProtParam tool⁹ (Gasteiger et al., 2003). To clarify the conserved motifs and protein domains, the online program MEME¹⁰ (Bailey et al., 2009) and NCBI database were used, respectively. The structure of GW2-like genes was predicted using TBtools software (Cheng et al., 2020).

Wild-type cotton (TM-1) and Arabidopsis (Col-0) were planted in the experimental field and greenhouse of the Institute of Cotton Research of Chinese Academy of Agricultural Sciences under conventional field management conditions in Anyang, China. To detect the relative expression of GW2-like genes, different samples were collected from WT cotton at different stages of ovule development. All ovule samples were self-pollinated. Total RNA was isolated from samples using Fast Pure Plant Total RNA Isolation Kit (Vazyme, Nanjing, China). The qRT-PCR analysis primers were listed in Supplementary Table 1. Gh_D03G0370 (GhActin3) (Li et al., 2005) and AT3G18780 (AtActin2) (Wang et al., 2016) were used as internal controls for the qPCR experiments.

The 1,272 bp complete coding sequence (CDS) of GhGW2-2D was amplified using the primers 35S:GhGW2-2D-F and 35S:GhGW2-2D-R (Supplementary Table 1). The 35S:GhGW2-2D construct were cloned into a pBI121 vector digested with BamHI and SacI using a PCR-based homologous recombination system. The 35S:GhGW2-2D plasmid transferred into Col-0 plants using Agrobacterium GV3101 and the transformants were selected on medium supplemented with kanamycin (50 mg/L). Three homozygous T₃ generation lines were used for further analysis. The seed and cotyledon were conducted by scanning to generate a digital image and then measured by ImageJ software.

In order to identify all GW2-like proteins in tetraploid cotton (G. hirsutum and G. barbadense) and their diploid ancestors (G. arboreum and G. raimondii), the protein sequences of AtDA2 and OsGW2 in rice and Arabidopsis were used as queries. We identified three, three, six, and six GW2-like genes in G. arboreum, G. raimondii, G. barbadense, and G. hirsutum, respectively. The length of these cotton GW2-like gene sequences ranged from 1,260 to 1,290 bp, encoding polypeptides from 419 to 429 amino acids long. The predicted isoelectric point of proteins varied from 4.70 to 4.93. According to the results of evolutionary analysis, all members were divided into three subgroups. More details of the gene names, locus IDs, and other features are shown in Table 1.

All GW2-like genes were unevenly mapped onto chromosomes based on the genome annotation data of the four species (Supplementary Figure 1 and Table 1). In the G. arboreum genome, the three identified GaGW2s were three located on chromosomes A01, A05, and A13. In G. raimondii, the three GrGW2 genes were distributed on three chromosomes: chromosome 02 (D01), chromosome 09 (D05), and chromosome 13 (D13). In G. hirsutum, six GhGW2 genes were mapped on three At subgenome chromosomes (A01, A05, and A13) and three Dt subgenome chromosomes (D01, D05, and D13). The chromosomal distribution of GW2-like genes in the G. barbadense genome was resemble to that of GhGW2s. Therefore, it could infer that GW2-like genes were highly conserved in the evolution from diploid to tetraploid in cotton.

An efficient way to predict gene function is to analyze its gene structure and the encoded protein domains. In this study, the GW2-like gene structures of four cotton species were visualized to further understand their differences and potential functional diversity.

A phylogenetic tree of GW2-like genes in G. hirsutum and G. barbadense and their ancestral diploids was constructed using the NJ method (Figure 1A). All GW2-like proteins were divided into three subgroups. Each group contained six members, one from G. arboreum, one from G. raimondii, two from G. hirsutum, and two from G. barbadense. All GW2-like genes have eight exons that were highly conserved in the four species of cotton (Figure 1B).

The NCBI and SMART databases were used to predict the domains of the GW2-like proteins. All GW2-like proteins had a typical E3 ubiquitin ligase RING_Ubox domain (Figure 1C), which can bind to a specific E2 ubiquitin-binding enzyme and determines the specificity of the ubiquitinated substrate. The MEME online program was used to predict the motifs of the GW2-like family genes. The results showed that most members in the same subgroup had resemble gene structures and functional domains (Figure 1D).

Gene replication and subsequent functionalization are important driving forces for genome and species evolution. In addition, gene duplication events are important for amplification of gene families. To clarify the amplification mechanism and evolution of GW2-like genes in the four cotton species, we conducted interspecific collinearity analysis and gene repeat event analysis (Figure 2 and Table 2).

Collinearity analysis was performed for G. arboreum and G. hirsutum (A subgenome), G. arboreum and G. barbadense (A subgenome), G. raimondii and G. hirsutum (D subgenome), G. raimondii and G. barbadense (D subgenome), respectively (Figure 2). A total of three, six, five, and seven pairs of collinear genes were found, respectively. The GbGW2-like genes in G. barbadense had more collinear gene pairs with G. arboreum (A2) and G. raimondii (D5), while the GhGW2-like genes in G. hirsutum had fewer collinear gene pairs with G. arboreum (A2) and G. raimondii (D5), indicating that GW2-like genes in G. barbadense were more related to two possible ancestral species of G. arboreum and G. raimondii than GW2-like genes in G. hirsutum.

In order to elucidating the evolutionary dynamics and selection pressure of protein coding, the non-synonymous (Ka) and synonymous (Ks) replacement rates were calculated using full-length sequences between each pair of linear homologous genes. Generally, a Ka/Ks ratio>1, = 1 and <1 indicates diversified selection, neutral selection, and purification selection, respectively. The Ka/Ks ratio of all gene pairs was less than 1, suggesting that these genes had undergone purification selection (Table 2). Above results indicated that the GW2-like genes in G. hirsutum and G. barbadense were relatively conserved.

To investigate the tissue-specific expression profiles of GW2-like genes, published TM-1 expression data, including root, stem, leaf, torus, petal, stamen, pistil, calycle, fiber and ovule at various developmental stages were used. Most GhGW2 genes were extensively expressed in many tissues (Supplementary Figure 2). GhGW2-1A and GhGW2-1D were relatively highly expressed in pollen, petal, and pistil. GhGW2-2A was highly expressed in 20 DPA fiber and roots. GhGW2-2D was highly expressed in the stem, root and early stages of ovule. GhGW2-3A and GhGW2-3D had similar expression trends with GhGW2-1A that highly expressing in flowers and pistils (Supplementary Figure 2A). The expression patterns of GbGW2-like genes in G. barbadense and G. hirsutum shared alike (Supplementary Figure 2B), indicating that these genes may have similar functions. qRT-PCR of diffident ovule development stages showed that GhGW2-1A, GhGW2-2A and GhGW2-1D had a similar tendency with highly expressed in early stages and later stages, and less expressed in 10–20 DPA ovules (Figures 3A–C). GhGW2-2D expression was highest in 10 DPA ovule and 30 DPA ovule, suggesting it may involve in cotton ovules development (Figure 3D). GhGW2-3A revealed an extensive expression in all ovule development stages (Figure 3E). GhGW2-3D was expressed at very low levels in all ovule development stages except in 30 DPA ovule (Figure 3F).

GhGW2-like genes had distinct expression patterns in different tissues of TM-1, suggesting that GhGW2-like genes may have different functions in cotton growth.

Based on the sequence alignment and homology analysis, GhGW2-2D was selected for further study because of its similarity to OsGW2 and AtDA2 (Figure 4). Moreover, the expression of GhGW2-2D was relatively high in the stage of ovule development. Plants overexpressing GhGW2-2D were generated and identified using semi-quantitative PCR (Figure 5A), and a total of three transgenic lines with high expression of GhGW2-2D were identified. qRT-PCR was used to further assess the relative expression levels of the three diverse transgenic lines (Figure 5B).

To evaluate the applicability of transgenic breeding to alter seed size, different developmental stages of transgenic Arabidopsis were characterised. Overexpression of GhGW2-2D reduced seed size (Figures 6A–D). The seed areas of lines 1, 2, and 3 decreased by 25, 24, and 26%, respectively, compared with that of the wild type Col-0 (Figure 6I). The growth of nine-day-old seedlings was measured after germinating. The transgenic lines had smaller cotyledons than Col-0 (Figures 6E–H). Compared with that of wild type Col-0, the cotyledon area of lines 1, 2, and 3 decreased by 31, 33, and 42%, respectively (Figure 6J). These results further support the role of GhGW2-2D in limiting seed and organ growth.