The variety Nannong 513 provided by Nanjing Agricultural University, was a temperature-sensitive variety. The seeds were grown in pots with peat soil and vermiculite (7:3), the germination was conducted in a plant incubator, and the temperature was set at 25°C. Four treatment groups were set in this experiment with three independent biological replicates. The plants were treated at the V1 stage of soybean (the first trefoil stage in the vegetative growth period of soybean), and the treated plants were pretreated with foliar application of exogenous melatonin (100 μmol/L) at night, while the other plants were sprayed with distilled water as control. After two days of spraying, half of the plants in both groups were put in another plant incubator at 4°C to conduct cold stress, while the other plants were still kept at 25°C as normal temperature control (Gai et al., 2020). Plant leaves were sampled at 24 h after cold treatment for transcriptome analysis and qRT-PCR detection, and 3 d after cold stress for physiological and biochemical determination.

The determination of the activity of peroxidase (POD) and catalase (CAT), and the glutathione redox homeostasis (GSH/GSSG), hydrogen peroxide (H₂O₂), and malondialdehyde (MDA) were conducted with POD, CAT, GSH, GSSG, H₂O₂, and MDA assay kits, respectively (Suzhou Grace Biotechnology Co., Ltd.). Evans Blue staining was performed according to the method of Xia et al. (2009).

The melatonin (MT) content in soybean was determined according to Yan et al., 2019, leaves were ground to a powder in liquid nitrogen and then extracted with 1.5 mL of chloroform at 4°C for 15 h. After centrifugation of the extraction mixture, the chloroform fraction was evaporated to dryness and dissolved in 100 µL of 42% methanol. Aliquots of 10 µL were subjected to HPLC using a fluorescence detector system. The contents of abscisic acid (ABA) and indole-3-acetic acid (IAA) were determined by ELISA reagent kits (Shanghai Enzyme-linked Biotechnology Co., Ltd. Shanghai, China) (Yu et al., 2022). The data were subjected to analysis of variance with SPSS, and the means were compared using Duncan’s t-test at the 5% level. meanwhile, the data was visualized using origin 8.0 (microcal Inc, Northampton, mA, USA).

Input 118 known A. thaliana B3 gene IDs into the Arabidopsis genome database (TAIR, http://www.arabidopsis.org/, 2020) to obtain protein sequences. Candidates from the soybean B3 family were investigated with BLASTP using Arabidopsis thaliana B3 protein sequence as probe. In the SMART database (http://smart.embl-heidelberg.de/, 2020), the domain of Pfam (PF02362) was identified and screened in the phytozome database. Finally, 145 GmB3 genes were identified and named GmB3-1–GmB3-145 according to their chromosomal positions. The theoretical isoelectric points (pI), the number of amino acids, and the grand average of hydropathicity (GRAVY) of all predicted B3s was then determined by ExPASy (Wilkins et al., 1999). Multisequence matching and Maximum-Likelihood analysis were performed using 1000 replicates as bootstrap values (Morgulis et al., 2008).

The MEME tool (http://meme.nbcr.net/meme/, 2020) was used to detect the B3s’ motifs (Bailey et al., 2009). The exon-intron sequence of B3 genes was determined by comparing the coding sequence and genome sequence of B3 genes using the Gene Structure Display Server (GSDSv2.0; http://gsds.cbi.pku.edu.cn/, 2020) (Hu et al., 2015).

Coordinate correspondence between DNA and protein sequences is determined using Gene-wise (Simmons et al., 2019). Whole-genome protein sequences and gene positions for soybean were retrieved from EnsemblPlants (http://plants.ensembl.org/index.html, 2020), MCScanX was utilized to analyze B3s gene duplication events (Wang et al., 2012). TBtools were adopted to visualize the results.

The upstream 1.5 kilobases (kb) genomic DNA sequences of the GmB3s were retrieved from the soybean genome and the putative cis-regulatory elements in promoter regions were identified using the PlantCare database (http://biinformatics.psb.ugent.be/tools/plantcare/, 2020) (Rombauts et al., 1999).

The B3 gene expression data in different tissues were obtained from the Phytozome database (https://phytozome-next.jgi.doe.gov/pz/, 2022).

Plant leaves were sampled at 24 h after cold treatment for transcriptome analysis, each treatment was repeated with three times (Yu et al., 2021). The total RNA was extracted, and the library was built, and the quality was tested, and the transcriptome was sequenced, analyzed and annotated by LC SCIENCES (Hangzhou, China) (Sheng et al., 2017). During the detection of DEGs, a fold-change ≥ 1.5 and a false discovery rate < 0.05 were used to screen different expression genes. Differential expression profiles of differentially expressed genes were presented. Amazing HeatMap software (Chen et al., 2020) was used to generate a heatmap.

Total RNA was extracted using TRIzol^® reagent (Invitrogen, Carlsbad, CA, USA). One microgram per RNA sample was used as the template for the synthesis of the first-strand cDNA, using the ReverTra Ace™ qPCR RT Master Mix with gDNA Remover (TOYOBO Co., Osaka, Japan). Subsequently, qRT-PCR was performed with SYBR^® Select Master Mix RT-PCR System (Takara) on an optical 96-well plate. Select actin as an internal reference. All the primers used for gene expression analysis were shown in Table S1 . Relative expression level of genes was calculated using formula 2^−ΔΔCT (Livak and Schmittgen, 2001; Cui et al., 2019). Three independent biological replicates were analyzed.

We utilized PlantRegMap (Internet) http://plantregmap.gao-lab.org/index- chinese.php (Tian et al., 2020) to predict transcription factors associated with GmB3. Simultaneously, the transcription factors (TFs) network was visualized using Origin 8.0.

To assess the cold-caused injury to soybeans, we analyzed its Chlorophyll and H2O2 content, MDA accumulation, and electrical conductivity. The cold treatment significantly increased the content of Chlorophyll, MDA, and electrical conductivity by 23.4%, 29.9%, and 263.6%, respectively, compared with the control. However, exogenous melatonin application substantially decreased the content of MDA and electrical conductivity by 11.4% and 46.7%, respectively, compared with the only cold treatment, and chlorophyll content was further increased. Interestingly, the content of H2O2 in soybean with foliar application of melatonin decreased by 28.8%, compared with cold treatment alone ( Figures 1A–D ). This indicated that exogenous melatonin treatment could alleviate the damage caused by cold stress in soybean.

We measured the activity of antioxidant enzymes including POD and CAT to unveil how melatonin alleviates cold-induced oxidative stress. The cold treatment significantly enhanced the activity of POD and CAT by 97.9% and 54.5%, respectively, compared with the control ( Figures 1E, F ). While compared with cold treatment alone, the activity of POD and CAT increased by 26.4% and 59.3%, respectively, with foliar application of melatonin. In addition, compared with the control, the GSH/GSSG ratio of soybean leaves after cold treatment decreased by 41.6%, compared with the control, but increased by 26.5% after melatonin application, compared with cold treatment alone ( Figure 1G ). The results of Evans blue staining also demonstrated that the damaged tissues in soybean leaves increased after cold treatment, while the application of melatonin alleviated the damage to soybean leaves ( Figure 1H ). Soybean seedling growth was seriously affected by the cold, compared with the control, the development of soybean seedlings was significantly inhibited and the leaf surface was damaged after cold stress, while melatonin treatment alleviated the leaf surface damage of soybean after cold stress ( Figure 1I ). This indicated that exogenous melatonin treatment could improve the antioxidant enzyme activity and GSH/GSSG ratio of soybean against cold stress.

To evaluate the effect of exogenous melatonin on endogenous melatonin and other hormones in soybean under cold stress, the content of melatonin, indoleacetic acid, and abscisic acid in the leaves of soybean seedlings were analyzed under different treatments ( Figure 2 ). Compared with the control, the content of endogenous melatonin in soybean seedlings after exogenous melatonin treatment increased in both normal temperature and cold treatment groups, especially for the cold treatment, the melatonin content in soybeans increased more significantly with the application of melatonin. Although no significant difference appeared in melatonin content in soybeans between the cold treatment and control. Interestingly, compared with the control, the cold treatment increased the content of IAA and ABA in soybeans, while both were further increased with the application of melatonin. The results suggested that exogenous melatonin might alleviate soybean response to cold stress by regulating the hormone content in soybeans under cold stress.

B3 transcription factors play an important role in plant tolerance to abiotic stress. To further study the effect of exogenous melatonin on cold stress, the expression of several representative soybean B3 genes was determined after treatment with melatonin under cold stress by qRT-PCR ( Figure 3 ). The results explored that the expression of B3-001, B3-003, B3-078, B3-083, B3-093, B3-095, B3-108, and B3-123 was induced by cold treatment and increased by 3.36-fold, 1.03-fold, 2.04-fold, 2.16-fold, 29.52-fold, 80.52-fold, 1.94-fold, and 1.38-fold, respectively, compared with the control. While the application of exogenous melatonin substantially further improved the expression of all the above genes, compared with only cold treatment, which up-regulated by 52.4%, 58.1%, 25.4%, 130.4%, 15.4%, 240.5%, 199.7%, and 141.2%, respectively. Similarly, the transcription of B3-049 gene was also significantly induced with the application of melatonin, compared with the cold treatment alone, although no significant difference was observed between cold and normal temperature treatment. There was no significant difference in the expression of the above genes between melatonin treatment and control under normal tmperature, except for B3-078 gene. The above results indicated that melatonin could alleviate the cold stress of soybean by inducing the expression of B3 genes.

We used Arabidopsis B3 protein sequence as reference, screened candidate B3 protein through blast comparison of soybean genome, and screened 145 B3 genes through conservative domain (PF02362) and redundancy elimination. The candidate genes were named GmB3-001 to GmB3-145 based on the location on chromosomes ( Figure 4 ). The soybean B3 genes were unevenly distributed on 20 chromosomes, among which the 20th chromosome was focused for the distribution of 12 soybean B3 genes. The B3 protein sequence consisted of 72-1136 amino acids, with an average length of 538 amino acids. The relative molecular weight and isoelectric point of B3 protein ranged from 8632.16 kDa (GmB3-140) to 127058.28 kDa (GmB3-117) and 4.86 (GmB3-100) to 11.02 (GmB3-066), respectively. Among them, 75 B3 proteins were acidic (with an isoelectric point < 7), while 70 B3 proteins were alkaline (isoelectric point > 7) ( Table S2 ).

The phylogenetic tree was constructed using the B3 protein sequences of soybean. The B3 family of soybean could be classified into four subgroups (namely B3-I, B3-II, B3-III, and B3-IV) based on the phylogenetic tree analysis. These B3 groups (I-IV) consisted of 57, 25, 11, and 52 members, respectively ( Figure 5A ). To elucidate the evolutionary relationship between the B3 gene family, we constructed a co-linear map of B3 genes in Arabidopsis and soybean, which revealed that 20 pairs of B3 genes presented collinearity between soybean and Arabidopsis ( Figure 5B ), and 26 pairs of B3 genes existed collinearity among 14 chromosomes in soybean ( Figure 5C ).

To identify the conserved structure of soybean B3 protein, 30 motifs were predicted by MEME motif analysis. Similar motif structures were present in the B3 members of the same subgroup, although different motifs appeared in various subgroups. For instance, most B3 proteins of subgroup contained 7 to 25 motifs, except for GmB3-35 and GmB3-85, containing one motif. Most GmB3 proteins of subgroup II contained 2 to 8 motifs, while Motif 30 was only found in subgroups II. Most B3 proteins of subgroup III only contained 2 motifs, although GmB3-60 and GmB3-77 contained one motif. And B3 proteins of subgroup IV contained 1 to 7 motifs ( Figures 6A, B ). The B3s structure of exon and intron was detected to obtain gene structure. Although different gene structures were exhibited in four subgroups, a relatively consistent gene structure of B3s was found in each subgroup ( Figure 6C ).

To investigate the regulatory mechanism of the B3 gene, we scanned the sequence of the promoter codon 2000bp upstream of ATG and obtained a number of cis-acting elements associated with plant hormones and stress response. Hormone response elements were mainly induced by auxin (AuxRR-core), gibberellin (GARE), and abscisic acid (ABRE). Among them, 84 B3 genes were found with ABA response elements, accounting for 57.9% of the total genes. In addition, 7 B3 genes with auxin response elements and 15 B3 genes with GA response elements were found. The predicted stress response elements mainly included LTR (low-temperature response), MBS (drought induction), and ARE (anaerobic induction), among which, 103 B3 genes containing ARE elements that were the most frequent, accounting for 71.0% of the total, while 34 B3 genes containing MBS and 26 B3 genes containing LTR were revealed ( Figure 7 ). These results suggested that B3 genes might be involved in the response of soybean to hormones and abiotic stresses.

The transcription level of B3 in different tissues (including roots, stems leaves, and flowers) was analyzed based on the Phytozome database. Remarkable difference appeared in B3s expression levels of various soybean tissues ( Figure 8 ). Thereinto, the expression of GmB3-116 was only high in flower, while the expression of GmB3-093 was only high in stem. However, GmB3-123 was only highly expressed in root, while GmB3-078 was only weakly expressed in root. These results suggested that B3 genes were expressed in a tissue-specific manner in the soybean.

Based on the transcriptome data, the response of B3 gene family members to cold and melatonin was analyzed. Compared to control, 39 B3 genes were induced and 59 B3 genes were inhibited by cold stress. The transcription of 40 B3 genes of soybean was induced, however, the expression of 57 B3 genes was reduced by melatonin treatment compared with the control. Similarly, the transcription of 42 B3 genes was enhanced, while the expression of 57 B3 genes was decreased by cold and melatonin treatment compared to the cold treatment alone ( Figure 9A ). It suggested that B3 genes responded to melatonin under cold treatment in soybean.

To further investigate the regulatory mechanism of the B3 gene of soybean, we predicted that a total of 34 transcription factor families interacted with the B3 gene. Among all B3 genes, 18 transcription factors were predicted to interact with GmB3-004 and GmB3-013, with the largest number ( Figure 9B ). Interestingly, transcriptome data revealed that 14 transcription factor families contained genes that could be induced by cold, while regulated by melatonin treatment ( Figure 9C ). The results indicated that the B3 gene regulated soybean stress response through interaction with transcription factors, suggesting that the B3 gene played an important role in plant resistance to abiotic stress through interaction with transcription factors.