The CPP gene and encoded proteins in various species like G. max, O. sativa, Z. mays, A. thaliana, Brassica rapa, G. soja, Cajanus Cajan, Chlamydomonas reinhardtii, and Selaginella moellendorffi were downloaded from plant transcription database.¹ To confirm the retrieved CPP proteins, local BLASTp, NCBI Batch CD-search, Interproscan V. 63² and SMART³ were also used. Non-redundant gene members were selected and the rest were excluded for further analyses. Other biophysical characteristics i.e., protein length, molecular weight (MW), isoelectric point (pI), and gravity values for GmCPPs were extracted using ExPASy Protparam Tool.⁴ Furthermore, sub-cellular localization of GmCPPs was also identified using the Softberry⁵ and CELLO V2.5 web-tool.⁶

Full-length amino acid sequences of all studied species were aligned and two phylogenetic trees were generated using MEGA X using the maximum likelihood method (ML) following parameters as reported by Kumar et al. (2018). The bootstrap method (1,000 replications) was used to determine the dependability of clades. Graphical representation of multiple sequence alignment of conserved CPP amino acid residues in G. max and G. soja was performed separately by the Clustal W program (Tamura et al., 2011) and WEBLOG webtool.⁷

To explore exon/intron structure, bed-files from databases were obtained and analyzed using GSDS 2.0.⁸ Protein motif distributions were determined using the online MEME tool.⁹ For cis-element analyses, ∼2 kb upstream regions were analyzed in the PlantCARE database (Lescot et al., 2002) and the elements were characterized on the basis of their predicted biological functions, and graphical representation was done by using TBtool software.

To determine the chromosomal distribution of GmCPPs, extracted gff3-files of soybean genome annotation were downloaded from SoyBase (SoyBase.org). Gene duplication analyses were performed following the methods as reported previously (Yang et al., 2017). CIRCOS was used to create the figure and Ka/Ks values were calculated using PAL2NAL (Suyama et al., 2006; Krzywinski et al., 2009).

A set of 46 soybean accessions were planted at MNS University of Agriculture, Multan in the springs of 2021 (2021-UAM). Field experiments were carried out under WL and WW experimental units following Augmented design (Check = UAMSB200). The WL experimental units were subjected to drought especially at the flowering stage, whereas, WW experimental units were irrigated after every fortnight (depending upon water requirement). Each soybean genotype was planted on two beds (length × width = 15 × 2.5 ft) on both sides. Plant-to-plant distance was maintained at a distance of 1 ft. Phenotypic data were recorded for plant height, thousand seed weight, pods^–1 plant, seeds^–1 pod, seed weight^–1 pod, seed length, seed thickness, seed width, and pod length from both water regimes.

Genomic DNA of the investigated soybean germplasm was extracted from young seedling leaves using the CTAB method (Lecharny et al., 2003). The quality of extracted DNA was initially checked by using NANO-Drop (K5800C Micro-Spectrophotometer) followed by running the extracted DNA on 1.0% agarose gel.

The whole genome sequence of three cultivars of soybean (Williams-82, Lee, and Zhonghuang-13) was downloaded from SoyBase.¹⁰ Local BLAST was performed to identify the sequences of GmCPPs in the aforementioned soybean genotypes. For the identification of sequence polymorphism, multiple sequence alignment was performed using the Seqman program in the DNASTAR Lasergene package. Standard kompetitive allele-specific PCR (KASP) guidelines¹¹ were followed for the development of KASP primers on the identified single nucleotide polymorphism (SNP) of GmCPP6. Allele-specific primers were developed having standard HEX and FAM tails with a targeted SNP at three prime ends. Two reverse primers (allele-specific) and one common forward primer were designed so that the total fragment length was less than 100 bp. The standard KASP reaction mixture, KASP assay, and PCR conditions were followed as reported by Rasheed et al. (2016), Majeed et al. (2018), Ur Rehman et al. (2019, 2021), and Irshad et al. (2019, 2021).

Phenotypic data were analyzed with XLSTAT Software 2014. Student’s t-test at p less than 0.05 was used to check the effect of each allelic variation on the recorded phenotypic traits.

We identified a total of 81 CPP genes in nine investigated species including chlorophytes (C. reinhardtii), lycophytes (S. moellendorffii), Brassicaceae (A. thaliana and B. rapa), Fabaceae (G. max, G. soja, and C. cajan), and Poaceae (Z. mays and O. sativa). Among these, 12 CPP genes were shortlisted in G. max, 10 in G. soja, 16 in B. rapa, 11 each in Z. mays and O. sativa, eight in A. thaliana, six in C. cajan, four in S. meollendorffii, and three in C. reinhardtii. A higher number of CPPs were identified in G. max as compared to chlorophytes and lycophytes indicating a duplication effect on GmCPPs in G. max. These findings also signify that CPPs experienced extension in higher plants. The transcription factor ID, taxonomic ID, and predict sub-cellular localization are presented in Supplementary Table 1. These results showed that the GmCPP coding sequence ranged from 1,656 to 2,715 bp for GmCPP-6 and GmCPP-11, respectively. Similarly, an amino acid number of GmCPP genes ranged from 483 to 904 for GmCPP-5 and GmCPP-11, respectively. Molecular weight ranged from 54,034.67 to 98,476.32 kDa for GmCPP12 and GmCPP-4, respectively. The isoelectric point of GmCPP4 was the highest (9.2) and that of GmCPP-1 was the lowest (5.41). The grand averages of hydropathicity values of all GmCPPs were less than zero and ranged from -0.734 for GmCPP-11 to -0.511 for GmCPP-4. In addition, all GmCPPs are localized in the nucleus, except GmCPP-1 and 10.

The phylostratum analyses of CPP genes identified the primitive lineage as CPP genes, which were also identified in chlorophyte (C. reinhardtii) (Figure 1). Further, the CPP genes were identified in lychophytes, dicots, and monocots. These outcomes signified that CPPs originated from early plants phylostratum and possible orthologs are present throughout kingdom Plantae. An evolutionary tree was generated to determine the phylogenetic relationship among the studied CPPs. To indicate the CPPs from C. reinhardtii, S. moellendorffii, A. thaliana, B. rapa, G. max, G. soja, C. cajan, Z. mays, and O. sativa the prefixes Cr, Sm, At, Br, Gm, Gs, Cc, Zm, and Os were used correspondingly. The phylogenetic analyses divided 81 genes into six clades based on sequence similarities (Figure 1). Clade-I comprised 16 members, Clade-II possessed 13 members, Clade-III contained 14 members, Clade-IV, and V had 15 members each, and Clade-VI contained eight members. Clade-I lacks genes from chlorophytes which suggested the evolution of CPP genes after the split of chlorophyte. Interestingly, CPP genes from monocot and dicot species were unsystematically distributed to all clades. Further, phylogenetic analyses indicated that G. max and B. rapa experienced gene family expansion since both have more CPP genes compared to other studied organisms.

To explore the conservation of each amino acid residue in GmCPP and GsCPP, multiple sequence alignment was executed to generate sequence logos in G. max and G. soja. The outcomes showed that the amino acid residue distribution was highly similar at most of the loci among the G. max and G. soja. For example, some amino acid residues such as C [6], L [7], Y [8], C [9], C [11], F [12], A [13], N [29], A [34], and so on were found to be highly conserved (Figure 2). Phylogenetic analyses also highlight that GmCPP and GsCPP members lie in close proximity to the evolutionary tree.

It has been well-documented that intron-exon distribution arrangement in a gene is related to its biological function. The intron number in GmCPP ranged from 7 to 9 (Figure 3). Conserved domains in each sequence were identified using the CDD tool of NCBI.¹² All members of the GmCPP gene family contain the TCR domain (Supplementary Figure 1). MEME tool was used to explore the conserved motif distributions of GmCPPs. The outcomes indicated that most of the GmCPP proteins exhibited similar motif distribution patterns such as motifs one, two, and eight exist in almost all proteins (Figure 4). We also identified cis-acting elements in the upstream regions of GmCPPs and grouped them on the basis of their functional relevance. All GmCPPs had cis-acting elements related to plant development, stress, and light responses (Supplementary Table 2).

The 12 GmCPP genes are scattered on five chromosomes, including three of each gene on chromosomes one and four, respectively. Two each gene are on chromosomes five and 10 while one gene is present on chromosome seven (Supplementary Table 3). To investigate the relationship of gene pairs, we explored the gene locus on a chromosome and executed synteny analyses. Their, synteny analyses showed that GmCPP genes were highly conserved among five chromosomes (Figure 5). Whole genome duplication, segmental duplication, and tandem duplication play a vital role in the extension of a gene family (Yang et al., 2017). To investigate the expansion of the GmCPP family in soybean, we executed gene duplication analyses in the soybean genome (Supplementary Table 3). Out of all studied gene pairs, 10 gene pairs were attributed to segmental duplication. We also explored the non-synonymous divergence (Ka) versus synonymous (Ks) values for the GmCPP gene pairs. It was found that nine duplicated gene pairs showed Ka/Ks values <0.5, whereas, two duplicated gene pairs showed Ka/Ks values between 0.5 and 1.0 (Supplementary Table 3). Generally, Ka/Ks of the studied gene pairs were <1, showing that the GmCPP gene family experienced purifying selection pressure with restricted functional differences.

For GmCPP-6, a polymorphic site was identified in the coding region. KASP marker was developed at the SNP site. KASP assay results showed that soybean accession having HEX tail has “C” allele while accession having FAM tail has “T” type allele (Figure 6).

Two allelic variations of GmCPP-6, i.e., GmCPP-6-T and GmCPP-6-C were identified in the studied soybean germplasm. GmCPP-6-C was the most frequently occurring allelic variation available in 58.6% of studied soybean accessions. Marker trait association analyses exhibited that at unique field sites, GmCPP-6-T was associated with higher thousand seed weight in both environments (Figure 7).

In the present work, we identified 81 CPP genes in nine different organisms i.e., chlorophytes (C. reinhardtii), lycophytes (S. moellendorffii), Brassicaceae (A. thaliana and B. rapa), Fabaceae (G. max, G. soja, and C. cajan), and Poaceae (Z. mays and O. sativa). Phylostratum analyses of CPPs indicated that the primitive plant pedigree as CPPs were existing in chlorophytes, showing that CPP genes originated from early land plants and ortholog genes of CPP are existing across kingdom Plantae. Phylogenetic analyses were performed to establish the evolutionary relationship among the studied species. All CPP genes were divided into six different clades which showed that most of the G. max genes exhibited a close relationship with G. soja genes and indicated that both species share a common ancestry. G. soja genome has shown to have 0.9154 GB consensus sequences, covering ∼98% of G. max genome sequence (Schmutz et al., 2010). Gene structure analyses showed that GmCPPs have a higher number of intron which indicates that GmCPP belong to the primitive gene family group. Sequence logos for conserved amino acid residues were also highly conserved in G. max and G. soja. Both N and C terminals of GmCPPs and GsCPPs are conserved. These results indicate that the GmCPP and GsCPP genes are evolutionarily conserved which might be helpful to underpin the pattern of CPP protein sequence conservation in other members of kingdom Plantae.

The estimates of biophysical parameters of all GmCPP gene family members delivered helpful information. Biophysical properties predicted that 10 out of 12 GmCPPs were positioned in the nucleus. The values of pI and grand average of hydropathicity of all GmCPPs indicated that all CPP proteins were hydrophilic (<0) and alkaline (∼7) (Supplementary Table 4).

The structure of the gene is a vital component that might be contributed by deletion and/or insertion incidents (Lecharny et al., 2003). In past, genome-wide studies have demonstrated that the loss or gain of introns during eukaryotic divergence was widespread (Rogozin et al., 2003; Roy and Penny, 2007). Gene structure analyses showed that all GmCPP genes have varied intron lengths that might play crucial roles in the functional divergence of GmCPP genes. It has been well-documented that introns play an important role in the evolution of different species (William Roy and Gilbert, 2006). In the current study, we observed that the number of intron for GmCPPs ranged from seven to nine indicating that G. max has evolved a long time ago (>Million years). Roy and Gilbert also advocated that earlier evolved species have more introns as compared with the newly evolved species (William Roy and Gilbert, 2006). Ten motifs were identified which showed that CPP proteins might function in different biological pathways allied with other co-factors. The motif distribution pattern of CPP proteins indicated that the distribution was relatively conserved and minimal divergence among the proteins from different groups might be linked with the specific biological function associated with soybean development and stress tolerance.

Transcription is governed by the binding of TFs to promoter cis-acting regulatory elements. Various studies have reported the crucial role of cis-acting elements in the processes of plant growth and stress responses (Fankhauser and Chory, 1997). In this study, cis-acting elements related to plant development and stress responses were identified in the upstream region of GmCPPs.

The uneven distribution of GmCPP genes on chromosomes of G. max shows probable gene loss or addition through whole genome or segmental duplication incidents. It has been reported that gene duplication and divergence generally lead toward evolution (Chothia et al., 2003). Gene duplication creates functional differences, which is crucial for speciation and adaptableness in changing environmental conditions (Conant and Wolfe, 2008). Gene duplication indicates that the aligned sequences share > 70% similarity and coverage length > 80% of the entire length (Yang S. et al., 2008). The two duplicated genes present on the different chromosomes of the same sub-genome might be the consequence of segmental or whole genome duplication, whereas, their presence on the same chromosomes might be the consequence of tandem duplication (He et al., 2012). Tandemly duplicated genes tend to be positioned together on chromosomes whereas, in segmental or whole genome duplication, the duplicated genes are generally distributed throughout the genome (Schauser et al., 2005). Approximately 65 million years ago whole genome and segmental duplications in primitive plant species contributed to the expansion of a number of gene families (Barakat et al., 2009; Wang et al., 2013) and contributed genomic complexity to kingdom Plantae (Cannon et al., 2004).

In the present work, we characterized 12 GmCPP genes, three times the number of CPPs present in chlorophytes, which indicates that CPP experienced expansion during their evolution. As reported previously, expansion in the genome permitted many crop plant species to acclimatize to environmental conditions (Ramsey and Schemske, 1998). We noticed that segmental and whole genome duplication were the major reasons responsible for the expansion CPP gene family in G. max. Segmental type duplication is the main contributor during evolution and it has happened in numerous plant genomes which contain many duplicated chromosomal blocks (Cannon et al., 2004). For instance, many A. thaliana gene families experienced evolutionary dynamics that led toward gene family expansion (Baumberger et al., 2003; Wang et al., 2008). Our results showed that GmCPPs grouped into pairs (Segmental duplication) which shows an ancient expansion in the gene family in G. max. To estimate the selection and environmental pressure, non-synonymous (Ka) and synonymous (Ks) rates of substitution (Ka/Ks) were computed. We noticed that Ka/Ks values of GmCPP genes were <1 illustrating that GmCPP gene family experienced strong purifying selection pressure.

Soybean PAN-genome might be helpful for bridging the phenotype to genotype gap in soybean breeding. Recently, PAN-genome has been used to discover genes for flowering time in G. soja (Li et al., 2014). Marker-assisted selection of elite alleles in breeding programs is vital for ongoing soybean breeding. The utilization of elite allelic variations in cultivars can be enriched if effective molecular platforms are available (Rasheed et al., 2017; Majeed et al., 2018; Ur Rehman et al., 2021). In this study, we used the PAN-genome of G. max to explore sequence polymorphism for GmCPP-6. Although, sequence polymorphism was explored in all studied GmCPP genes the allelic variation was only identified in the CDS region of GmCPP-6. Hence, all other genes were excluded for marker-trait association analyses. The absence of polymorphism in all other GmCPP genes is possibly due to allele fixation during evolution or because of the lower number of G. max accessions available for PAN-genomics studies. More work is required for further confirmation or to investigate these two possibilities. Converting sequence polymorphism to gel-free (KASP) markers enable SNPs to be more efficiently applied in selecting desirable alleles in marker-assisted breeding. Moreover, the KASP assay procedure is cost-effective. In the current study, soybean accessions having GmCPP-6-T had higher thousand seed weight under both environmental conditions i.e., WW and WL. Moreover, RNA-Seq Atlas of G. max also reported higher expression of GmCPP-6 in seeds (14–25 days after fertilization) (Severin et al., 2010). Generally, yield-related parameters of crop plants are administered by several genes and are strongly influenced by external stimuli. Pyramiding of favorable alleles might be helpful for continued improvement in soybean. The developed molecular marker will be useful for marker-assisted breeding in soybean which can be used in combination with other molecular markers.