Six Arabidopsis IQM protein sequences were used as reference, and their homologous sequences in soybean were retrieved from Phytozome v13 database (http://phytozome-next.jgi.doe.gov ) through Protein Basic Local Alignment Search Tool (BLASTP) program. The Pfamscan (https://www.ebi.ac.uk/Tools/pfa/pfamscan/ ) and SMART (http://smart.embl-heidelberg.de/ ) databases were used to annotate the conserved domains of candidate sequences. Finally, proteins containing the complete IQ motif were identified as members of the soybean IQM family. Information on soybean IQM genes, including CDSs, genome sequences, location coordinates, lengths of open reading frames (ORF), number of amino acids, and molecular weight, were acquired from Phytozome v13. The physicochemical characteristics of the GmIQMs were generated using ExPASy (http://web.expasy.org/protparam/ ). Subcellular localization was predicted using the WoLF PSORT program (http://wolfpsort.org ).

To examine the intraspecific and interspecific evolutionary relationships of IQM genes, the predicted IQM protein sequences of multiple species were obtained from the corresponding databases. Arabidopsis and rice IQM protein sequences were downloaded from TAIR (http://www.arabidopsis.org) and RAP-DB (http://rapdb.dna.affrc.go.jp/ ), respectively, whereas soybean, alfalfa, tomato, maize, sorghum, and Brachypodium distachyon IQM protein sequences were retrieved from Glycine max Wm82.a2.v1, Medicago truncatula Mt4.0v1, Solanum lycopersicum ITAG4.0, Zea mays B84 v1.2, Sorghum bicolor v3.1.1, and Brachypodium distachyon v3.1 of Phytozome v13. ClustalW software was used for multisequence alignment. The parameters were set to the system default values [Pairwise Alignment—Gap Opening Penalty:10.00; Gap Extension Penalty: 0.10. Multiple Alignment—Gap Opening Penalty:10.00; Gap Extension Penalty: 0.20. Use Negative Matrix: off; Delay Divergent Cutoff(30)]. MEGA11 software was used to construct an unrooted phylogenetic tree using the neighbor-joining (NJ) method, and bootstrap analysis was conducted using 1,000 replicates (Wu et al., 2016). The generated phylogenetic tree was identified using iTOL (https://itol.embl.de/ ).

To compare the gene structures of GmIQMs, the distribution of exons and introns was analyzed using GSDS2.0 (http://gsds.cbi.pku.edu.cn) (Mei et al., 2021). The gff3 format files of soybean IQM genes downloaded from Phytozome v13 database were imported, then the system will automatically generate the structure chart contained CDS, UTR and Intro. To understand the similarities and differences in protein structures, the conserved motifs of encoded GmIQM proteins were identified using the MEME database (https://meme-suite.org/meme/ ), the maximum number of motifs was set to 15.

The distribution image of GmIQM genes on the soybean chromosome was generated by MG2C (http://mg2c.iask.in/mg2c_v2.1/ ) according to the gene information retrieved from the Phytozome v13 database. Gene duplication analysis was performed as previously described (Feng et al., 2014). The SoyBase browser (http://soybase.org/gb2/gbrowse/gmax1.01) was used to search for duplicate gene pairs. Coparalogs were considered to be duplicated in tandem if they were on the same chromosome and separated by five or fewer genes in a 100-kb region (Wang et al., 2010); otherwise, they were deemed to be fragment duplications. To further analyze the divergence of duplicated genes, the synonymous substitution rate (Ks) and non-synonymous substitution rate (Ka) were calculated using TBtools software. According to a rate of 6.1×10^-9 substitutions per site per year, the divergence time (T) was calculated using the Ks value and the formula: T = Ks/(2×6.1×10^-9)10^-6 Mya (Lynch and Conery, 2000).

The Arabidopsis culture conditions were obtained from Lv et al, 2019. Arabidopsis (Arabidopsis thaliana) Col-0 seeds were sown in nutrient soil after 3 days of vernalization at 4°C and then cultured at 22°C under 16 h light/8 h dark conditions for approximately 4 weeks. Rosettes were cut before flowering to prepare protoplasts for subcellular localization and BiFC assays.

Soybean (Glycine max) Williams82 seeds were placed on absorbent paper and placed in the dark at room temperature for three days to germinate. Seedlings that grew consistently in a sponge hole tray were cultured in 1/4 strength Hoagland’s solution for 14 days under 16 h light/8 h dark conditions at 25°C and 6000 lx with 80% relative humidity. Subsequently, 50 mM NaHCO₃, 150 mM NaCl, 20% (w/v) PEG6000, 100 µM ABA, 100 µM MeJA, and 2 mM SA were added to Hoagland’s solution to simulate various abiotic stress conditions. Equal amounts of leaves and roots were collected for RNA isolation at 0, 3, 6, 12, and 24 h after treatment.

Yeast two-hybrid assays were performed according to the manufacturer’s protocol. Briefly, the full-length CDS of GmCaM (Glyma.19G121900) and 15 GmIQMs were cloned separately into pGADT7 and pGBKT7. The plasmids, AD-GmCaM and BD-GmIQMs, were co-transformed into AH109 yeast cells using the PEG/LiAC method. Interactions in yeast were tested on SD/-Trp/-Leu/-His and SD/-His-Trp-Leu-Ade plates containing 20 μg/mL x-α-gal. Cotransforming with the AD empty vector and BD-empty vector was used as a negative control, while AD-AtCAT2 and BD-AtIQM1 was used as a positive control (Lv et al., 2019).

For construction BD-GmIQM^LQ/VQ yeast expression vectors, LQ or VQ motif in CDS of five representative GmIQMs (GmIQM1a, 112-113aa LQ; GmIQM2a, 135-136aa LQ; GmIQM3a, 63-64aa VQ; GmIQM1e, 127-128aa LQ; GmIQM6a, 66-67aa LQ) was deleted using the PCR product fragment splicing method in vitro. Mutant protein sequences were cloned into the pGBKT7 vector. Thereafter, GmCaM binding to the mutant GmIQMs was analyzed using the yeast two-hybrid assay described above.

The protoplasmic transformation for both experiments was performed as previously described (Lv et al., 2019). For subcellular localization, the full-length CDS of 15 IQM genes in soybean Williams82 was cloned into the pGreenII-35S-GFP plant expression vector, respectively. Thereafter, pGreenII-35S-GFP-GmIQMs plasmids were transformed into protoplasts isolated from the leaves of 4-week-old Arabidopsis wild-type Col-0. After 16 h of incubation at 22°C, the GFP signal was observed using a Zeiss LSM800 confocal microscope at 488 nm absorption and 507 nm emission wavelengths. The pGreenII-35S-GFP empty vector was used as a negative control.

For the BiFC assays, the CDSs of GmIQM1d/GmIQM2c and GmCaM were cloned into binary pSAT1-nEYFP and pSAT1-cEYFP vectors, respectively. Thereafter, the pSAT1-nEYFP-GmIQM1d/GmIQM2c and pSAT1-cEYFP-GmCaM plasmids were transformed into protoplasts isolated from the leaves of 4-week-old Col-0. After 16 h of incubation at 22°C, the EYFP signal was observed using a Zeiss LSM800 confocal microscope at 488 nm absorption and 530 nm emission wavelengths. Co-transformation with nEYFP empty vector and cEYFP-GmCaM served as negative controls.

RNA was isolated using the Eastep Super Total RNA Extraction Kit (Promega, LS1040) according to the manufacturer’s protocol. For the tissue-specific expression analysis, total RNA was extracted from roots (no nodule), cotyledons, stems, leaves, flowers, pods and nodules of Williams82. For the gene expression analysis under abiotic and biotic stress, total RNA was extracted from 14-days-old leaves and roots of Williams82, with or without treatment. RT-qPCR was performed as previously described (Bu et al., 2014; Lv et al., 2019), and first-strand complementary DNA (cDNA) synthesis was performed using the PrimeScript RT Reagent Kit (Takara, RR047A), according to the manufacturer’s instructions. cDNAs were used as templates for RT-qPCR with gene-specific primers. Primer sequences are listed in Supplemental Table S2 . RT-qPCR was performed using 384-well plates with SYBR Premix Ex Taq II (Takara, RR820A) and a Roche Light Cycler 480 Real-Time PCR system. The soybean tubulin gene (Glyma.05G207500) was used as an internal control for mRNA (Li et al., 2017) and the relative expression levels of the genes were calculated using the 2^−CT method. Each sample was collected from five independent plants, independent biological assays were repeated three times, and similar results were obtained. Data from one representative replicate are shown. Values are means average from three technical measurements.

As revealed in previous studies, IQM proteins are plant-specific calcium-independent calmodulin-binding proteins that contain one IQ motif (Zhou et al., 2010). To identify the IQM gene family in soybean, we performed BLASTP analysis of the Glycine max genome via alignment with reported IQM proteins in Arabidopsis. Finally, 15 candidate GmIQM genes containing the IQ calmodulin-binding motifs were confirmed using conserved domain analysis. According to the order on chromosomes of GmIQMs and their homology relationship with AtIQMs, these 15 GmIQM genes were named GmIQM1a-1f, GmIQM2a-2c, GmIQM3a-3c, GmIQM5, and GmIQM6a/6b. The characteristics of each IQM protein, including location on chromosome, location coordinates, open reading frame (ORF) length, number of amino acids, molecular weight, theoretical pI (isoelectric points), and predicted subcellular localization, are listed in Table 1 . According to these data, the protein sequences of the 15 soybean IQM genes did not significantly differ in length (467-661 aa) and molecular weight (52.29-74.58 kDa). Further, most proteins had high isoelectric points (pI 7.14-9.34), except GmIQM2a (pI 6.42) and GmIQM2b (pI 6.24). Almost all GmIQM proteins were predicted to localize in the nucleus, while only few were predicted to localize in the chloroplast or cytoplasm ( Table 1 ).

Based on our previous study, IQMs can interact with CaMs in Arabidopsis and rice (Zhou et al., 2012; Fan et al, 2021 suggesting that binding to CaM is a common characteristic of IQM proteins. To further verify IQM in soybean, yeast two-hybrid experiments were carried out to determine whether IQM could also be combined with CaM in soybean. The results showed that only two vector combinations of AD-GmCaM and BD-GmIQM1d/2c could not grow on SD/-Trp/-Leu/-His triple dropout and SD/-Trp/-Leu/-Ade/-His quadruple dropout culture media ( Figure 1A ), indicating that most GmIQMs could interact with GmCaM in yeast cells, except GmIQM1d and GmIQM2c.

To further assess the interaction between GmIQM1d/GmIQM2c and GmCaM in vivo, an enhanced yellow fluorescent protein (EYFP)-based bimolecular fluorescence complementation (BiFC) assay was performed. Plasmids carrying nEYFP-GmIQM1d/GmIQM2c and cEYFP-GmCaM fused expression vectors were co-transformed into protoplasts of Arabidopsis Col-0. The co-transformation with the nEYFP empty vector and cEYFP-GmCaM was used as negative controls. Fluorescence signals were observed in the test group using a laser confocal microscope ( Figure 1B ).These results revealed that GmIQM1d/GmIQM2c interacts with GmCaM in vivo and a particular component of plant cells, which is nonexistent in yeast cells, and may be indispensable for their interactions. Overall, the reliability of soybean IQM proteins was confirmed using protein interaction assays.

To understand the similarity and evolutionary relationship between GmIQM proteins, we constructed an unrooted phylogenetic tree of 15 soybean IQM protein sequences. Based on the results, all members of the soybean IQM gene family could be divided into three major subfamilies (I, II, and III; Figure 2A ). Similar to the analysis results for Arabidopsis and rice, the IQM members in soybeans mainly belonged to subfamilies I and II. Notably, subfamily I had the most members, with seven genes, while subfamily III had the fewest IQM members, with three genes. Further, the most closely related members within each subfamily formed sister-gene pairs.

The diversity of gene structure is well known as an important foundation for gene family classification (Zhao et al., 2017; Mengarelli and Zanor, 2021; Zhou et al., 2021). To further detect the features of each IQM subfamily in the soybean genome, the exon/intron organization of GmIQM genes was generated ( Figure 2B ). The structure chart revealed that most soybean IQM genes (eleven/fifteen) contained eight exons, except GmIQM1d, which contained ten exons, and three members of subfamily III (GmIQM3a, 3b, and 3c), which had nine exons. The relatively uniform exon numbers indicate the structural conservation of GmIQM.

The most closely related GmIQM gene pairs in the same subfamily shared similar gene structures, both in exon distribution and intron length. Nevertheless, one pair (GmIQM1d/GmIQM5) displayed an obvious discrepancy ( Figure 2B ); compared with GmIQM5, GmIQM1d had two extra short coding sequences in front of the first exon, which may be due to intron gain events during the long evolutionary period.

The conserved domains or motifs of 15 GmIQM genes were examined using the MEME database. Fifteen potentially conserved motifs were identified ( Figure 3A ; Figure S1 , and S2 ). Motif 5 was annotated to encode the IQ motif based on a data search of Pfamscan and SMART. The distribution of motifs in 15 GmIQMs is consistent with the gene clustering that

Most closely gene pairs in the same subfamilies had identical motifs, except GmIQM1d/GmIQM5 ( Figure 3A ), suggesting functional similarities or redundancies among these GmIQM proteins. Further, the 15 GmIQM proteins had common motifs (motifs 1, 2, 3, 5, 6, 7, 8, and 9); however, some motifs were recognized to be specific to certain subfamilies. For example, motifs 12 and 15 are unique to subfamilies II and I, respectively. Thus, variations in these specific motifs may cause functional differentiation of the IQM proteins in soybean.

The IQ motif is key to the interaction between CaM and CaMBP (Rhoads and Friedberg, 1997). To determine whether the combination of GmCaM and GmIQMs is dependent on the IQ motif, we constructed a BD vector with the soybean IQM protein sequence deleted 6Q (LQ/VQ) motif. Thereafter, AD-GmCaM and BD-^LQ/VQ were co-transformed into yeast cell AH109. As all transformants did not grow on the SD/-Trp/-Leu/-Ade/-His quadruple dropout culture medium ( Figure 3B ), the IQ motif was identified to be necessary for functional GmIQM and its conservation in different species.

To accurately understand the orientation of the 15 soybean IQM genes on each chromosome, a chromosomal map we constructed based on the location information retrieved from the soybean database. The distributions of these genes on chromosomes appeared to be wide but unbalanced. Fifteen GmIQM genes were mapped on 12 of the 20 soybean chromosomes. Nine chromosomes, including chromosomes 2, 5, 7, 10, 11, 13, 17, 18, and 19, contained only one IQM gene. Three chromosomes, including chromosomes 8, 9, and 15, had two IQM genes, while the remaining eight soybean chromosomes did not contain any IQM gene ( Figure 4 ). This biased distribution pattern of IQM genes has been observed in Arabidopsis and rice genomes (Zhou et al, 2010; Fan et al, 2021). Although two IQM genes were present on chromosomes 8, 9, and 15, however, there was no clustering from GmIQM genes on the same chromosomes, suggesting that the expansion of IQM family genes in soybean may not have been produced from tandem duplications.

Gene duplication is considered an important source of biological evolution. There are three main types of gene duplication: segmental duplication, tandem duplication, and transposition (Sankoff, 2001; Flagel and Wendel, 2009; Panchy et al., 2016). Among these duplication types, segmental duplication is a major contributor to the amplification of many gene families (Cannon et al., 2004; Zhu et al., 2020; Li et al., 2021; Wang M, et al., 2021; Zhao et al., 2021);. Herein, gene duplication events were examined to further understand the expansion mechanism of the IQM family in soybean. According to the data of synteny (Glycine recent duplication) from SoyBase browser, we found six pairs of

GmIQM genes located in duplicated blocks, indicating that these genes were generated by segmental duplication. The remaining three genes (GmIQM1f, 2c, and 3c) lacked duplicated pairs in their corresponding synteny blocks ( Figure 4 and Table 2 ), which aligns with the results of clustering analysis. We determined whether tandem duplication also played a role in adjacent GmIQM genes on the same chromosome. It has been reported that a pair of genes is separated by three or fewer genes within a 100-kb region on a chromosome, this may be due to tandem duplication (Wang et al., 2010; Feng et al., 2014; Wang X, et al., 2021). According to this criterion, no pair was generated by tandem duplication of the soybean IQM genes. Therefore, segmental duplication was identified to contribute significantly to the expansion of the IQM gene family in soybean.

To understand the evolutionary selection of duplicated GmIQM genes, we calculated the Ka/Ks ratio (substitution ratio of non-synonymous/synonymous) for each pair of duplicated GmIQM genes. In general, Ka/Ks = 1 indicates that both genes drift neutrally, Ka/Ks>1 indicates accelerated evolution with positive selection, and Ka/Ks<1 indicates a functional constraint with negative or purifying selection of the genes (Hurst, 2002; Nekrutenko et al., 2002; Feng et al., 2014; Ma et al., 2014). The Ka/Ks ratios from all duplicated gene pairs were found to be less than 0.4 ( Table 2 ), suggesting that the evolution of the soybean IQM gene family was mainly influenced by negative or purifying selection, thereby limiting the functional differentiation of duplicated genes. Based on the divergence rate of 6.1×10^-9 synonymous mutations per site per year (Lynch and Conery, 2000; Zhang et al., 2011; Chen et al., 2014; Hyun et al., 2014; Meng et al., 2015; Zhang et al., 2019), the duplication of these paralogous pairs was estimated to occur between 7.263 and 14.674 million years ago (Mya) ( Table 2 ).

Seventy-four putative IQM protein sequences from eight species, including four dicotyledons, Glycine max (15), Arabidopsis thaliana (6), Solanum lycopersicum (9), Medicago sativa (8); and four monocotyledons, Oryza sativa (8), Zea mays (11), Sorghum bicolor (10), and Brachypodium distachyon (7), were analyzed using the neighbor-

joining (NJ) method to further assess the phylogenetic relationship of IQMs in plants. All IQMs were found to cluster into three distinct subfamilies (I, II, and III) ( Figure 5 ). This result aligns with the classification of GmIQM genes; subfamilies I and II had the most members (77%) in the combined phylogenetic tree, while subfamily III had the fewest genes (23%). Each subfamily was further divided into two subgroups (a and b) of dicots and monocots ( Figure 5 ), indicating that significant differences in the structure and function of IQM genes between dicots and monocots.

The subcellular localization of 15 soybean IQM proteins was determined to further confirm their functional sites in cells. The full-length CDS sequence of each GmIQM was cloned and constructed into a pGreenII-35S-GFP vector, and different plasmids were transformed into protoplasts of Arabidopsis wild-type Col-0. After incubation, the GFP signal was observed under a confocal microscope ( Figure 6 ). The results showed that GmIQM1a and GmIQM1b were only located in the nucleus; GmIQM3a, 3b, and 3c were located in the nucleus and cytoplasm; and the remaining IQMs were mainly localized in the cytoplasm. These location patterns matched the clustering results, members of the closest relative clade had the same localization in the cell, indicating that these genes may have the same or similar functions. Further, the practical localization of these genes did not completely match the prediction by the Wolf PSORT programs.

Every stage of plant growth and development is closely regulated by a large number of genes. To understand the potential role of soybean IQM genes in the plant life cycle, quantitative polymerase chain reaction (qRT-PCR) was performed to determine the expression levels of IQM genes in soybean roots (no nodule), cotyledons, stems, leaves, flowers, seeds, pods, and nodules respectively. As depicted in the heat map, soybean IQM genes from the same subfamily or the most closely related paralogous gene pairs had similar gene expression patterns. For example,

Three members of subfamily III (GmIQM3a, 3b, and 3c) were significantly expressed in eight tissues. However, in subfamily I, only two genes, GmIQM1c and GmIQM1e, which are sister pairs in the closest clade, had universal expression, and the remaining genes (GmIQM1a, GmIQM1b, GmIQM5, GmIQM1d, and GmIQM1f) were not expressed in most tissues ( Figure 7 ). Notably, some GmIQM genes have tissue-specific expression characteristics, such as GmIQM5 and GmIQM1d, which had markedly high transcript abundance in nodules but were not expressed in other tissues, and IQM6a, which only had a slight expression in flowers ( Figure 7 ), suggesting that these genes might have vital functions at specific stages of soybean development.

The type and number of cis-acting elements in gene promoters are well known to largely determine gene function. To explore the potential role of soybean IQM genes, we analyzed the GmIQM gene promoter sequences within 2 kb upstream of the start codon using the PlantCARE database. In addition to the core elements, TATA-box and CAAT-box, various stress and hormone response elements were found in the promoter of each soybean IQM gene, such as MBS, GT1, ERE, MYB, MYC (drought, high salt, low temperature, etc.), ABRE (ABA), CGTCA-motif, TGACG-motif (MeJA), and TCA-element (SA) ( Figure 8 ). Therefore, GmIQMs could play a role in the response of soybean to biotic and abiotic stresses. However, the number and distribution range of abiotic response elements were significantly higher than those of MeJA and SA stress-related elements, especially MYB and MYC, which existed in almost all soybean IQM promoters. Thus, the soybean IQM gene family may be more important for abiotic stress response than biotic defense.

According to the analysis of cis-acting elements in promoters, IQM genes are most likely regulated by abiotic stresses and hormonal stimuli in soybean. To further confirm the function of the GmIQM genes, soybean wild-type Willims82 was treated with different stressors, including NaCl, NaHCO₃, PEG, SA, ABA, and MeJA, to simulate saline, alkaline, drought, and biological stress conditions. After treatment, the expression levels of eight representative IQM genes in soybean leaves and roots were determined using qRT-PCR. Gene expression profiling revealed that these IQM genes were strongly induced by NaCl and PEG, both in the leaves and roots ( Figure 9 ). Thus, these IQM genes play an important role in the response to high salinity and drought conditions in soybean. These genes were also upregulated by NaHCO₃, SA, ABA, and MeJA; however, compared with NaCl and PEG, the induction degree of their expression levels was relatively low under these stresses, and the transcriptional change mainly occurred in the roots, and not in the leaves ( Figure 9 ). These results indicated that IQM family genes also function in alkali and biological stress responses to some extent.