text:
Genome-wide association studies provide genetic insights into natural variation of seed-size-related traits in mungbean

Although mungbean (Vigna radiata (L.) R. Wilczek) is an important legume crop, its seed yield is relatively low. To address this issue, here 196 accessions with 3,607,508 SNP markers were used to identify quantitative trait nucleotides (QTNs), QTN-by-environment interactions (QEIs), and their candidate genes for seed length (SL), seed width, and 100-seed weight (HSW) in two environments. As a result, 98 QTNs and 20 QEIs were identified using 3VmrMLM, while 95, >10,000, and 15 QTNs were identified using EMMAX, GEMMA, and CMLM, respectively. Among 809 genes around these QTNs, 12 were homologous to known seed-development genes in rice and Arabidopsis thaliana, in which 10, 2, 1, and 0 genes were found, respectively, by the above four methods to be associated with the three traits, such as VrEmp24/25 for SL and VrKIX8 for HSW. Eight of the 12 genes were significantly differentially expressed between two large-seed and two small-seed accessions, and VrKIX8, VrPAT14, VrEmp24/25, VrIAR1, VrBEE3, VrSUC4, and Vrflo2 were further verified by RT-qPCR. Among 65 genes around these QEIs, VrFATB, VrGSO1, VrLACS2, and VrPAT14 were homologous to known seed-development genes in A. thaliana, although new experiments are necessary to explore these novel GEI-trait associations. In addition, 54 genes were identified in comparative genomics analysis to be associated with seed development pathway, in which VrKIX8, VrABA2, VrABI5, VrSHB1, and VrIKU2 were also identified in genome-wide association studies. This result provided a reliable approach for identifying seed-size-related genes in mungbean and a solid foundation for further molecular biology research on seed-size-related genes.

Background
Mungbean (Vigna radiata (L.) R. Wilczek) is a basic source of protein and carbohydrate, as it contains approximately 20% protein and 75% carbohydrate, and is a traditional and important legume in Asia (Somta et al., 2007). Due to its short life cycle (60-75 days), relative drought tolerance, and the ability to restore atmospheric nitrogen in association with Rhizobium/Bradyrhizobium bacteria, mungbean plays a crucial role in cropping systems and soil improvement (Somta et al., 2007;Alam et al., 2014).
The crop is generally grown as a cash crop in cerealbased farming systems. However, the major constraint in mungbean production is low seed yield. The average seed yield of mungbean is only approximately 700 kg per ha (Islam et al., 2015). Therefore, improving seed yield is the main goal in mungbean breeding. Understanding the genetic basis underlying seed-size-related traits is critical for the genetic improvement of mungbeans. In mungbeans, the ideotype of high-yielding cultivars are generally characterized by a large seed size, a short and synchronous maturity, a low sensitivity or insensitivity to day length, and the resistances to insects and disease (Fernandez et al., 1988). However, the knowledge on genes related to seed size has been limited. Moreover, the genes involved in the pathway of seed developments are not yet fully known.
Seed weight is the most important yield component and directly proportional to seed yield per plant in mungbean. To date, there have been seven studies of QTLs for seed weight in mungbean. Most of these studies are based on bi-parental segregation populations derived from interspecific crosses between cultivated and wild (V. radiata var. sublobata) mungbeans, and only two studies have evaluated seed size in more than one environments. The number of QTLs identified in those studies ranged from 3 to 11. Humphry et al. (2010) reported 11 loci for seed weight using SSR-marks, and Mei et al. (2009) identified a major QTL associated with both bruchid resistance and seed mass. Nonetheless, no candidate gene was identified for this trait.
Knowledge regarding seed development pathway is also a valuable source for transgenic strategies to improve crop production. As reported, there are several signaling pathways that control seed size, including the G-protein signaling, ubiquitin proteasome pathways, mitogen-activated protein kinase (MAPK) signaling, auxin pathways, and some transcriptional regulators (Li et al., 2019). In Arabidopsis, GPA1, AGB, and AGG3 were involved in G-protein-signaling pathways. DA1, DA2, SOD2, UBP15, EOD1, and SAMBA were involved in ubiquitin proteasome pathways. In addition, ABA2, ABI5, SHB1, MINI3, IKU2, and CKX were involved in the HAIKU (IKU) pathway. Additional genes were found to be related to seed size developments, but their pathways are uncertain, such as KIX8, BES1, MES1, and KLU (Orozco-Arroyo et al., 2015;Li et al., 2019). However, some reports have been focused on genetic foundation and molecular mechanism of seed developments in mungbean.
Genome-wide association studies (GWASs), along with multi-omics analysis, have been frequently used to mine candidate genes for most important agronomic traits in crops. Integrating GWAS with comparative genomics, transcriptome analysis, and molecular experiments, genes have been identified to be associated with complex traits (Liu et al., 2020c). For example, Gong et al. (2022) conducted a GWAS with high-quality single nucleotide polymorphism (SNP) data and seed-size traits, and found that Cla97C05G104360 and Cla97C05G104380, which are involved in abscisic acid metabolism, played important role in regulating the seed size in watermelon. Duan et al. (2022) identified GmST05 to be associated with soybean seed size through the GWAS of 1800 soybean germplasm resources, and GmST05 differed significantly at the transcriptional level. Liu et al., 2022a,c used GWASs and biological experiments to identify a pleiotropic gene GmPDAT for seed size-and oil-related traits in soybean, and a salt-stress-tolerance gene VrFRO8 in mungbean. Nonetheless, the related genes responsible for seed-size-related traits remained unknown in mungbean.
To address the above issues, 196 mungbean accessions with 3,607,508 SNP markers were used to conduct GWAS for seed length (SL), seed width (SW), 100-seed weight (HSW) using 3VmrMLM (Li et al., 2022b), efficient mixed-model association expedited (EMMAX) (Kang et al., 2010), genome-wide efficient mixed-model association (GEMMA) (Zhou and Stephens, 2012), and compressed mixed linear model (CMLM) (Zhang et al., 2010) methods. Candidate genes around quantitative trait nucleotides (QTNs) and QTN-by-environment interactions (QEIs) for the three traits were predicted by transcriptomics and comparative genomics. Key candidate genes were verified by RT-PCR analysis. Moreover, genes in seed-development-regulation pathway were also mined by comparative genomics. It should be noted that VrEmp24/25 and VrKIX8 were found to be associated with SL and HSW, and a major gene VrPAT14 (LOD = 61.95, r 2 = 5.80%) was identified in QEI detection via 3VmrMLM.

Plant materials and treatments
A diverse set of 196 mungbean accessions including 20 wild and 176 cultivated accessions from 23 countries, were used in this study (Supplementary Data Set 1). All the accessions were planted in a randomized complete block design with two replicates in an experimental field of Kasetsart University, Kamphaeng Saen Campus, Nakhon Pathom, Thailand in 2018 and 2020. In each replicate, each accession was planted in a single row 2.5 m long with 12.5 cm intra-row spacing (ca. 20 plants/row) and 50 cm inter-row spacing. Cultural practices were performed according to Park (1978). SW (mm), SL (mm), and HSW (g) were measured. At maturity. The SL and SW traits for each accession were averaged based on 20 seeds and 100SW for each accession was averaged based on three replicates.

Whole-genome resequencing
The young leaves of the above 196 mungbean accessions were collected 1 week after planting. The DNA was extracted in 2018, using the CTAB method (Smith et al., 2005). Short reads sequenced by an Illumina HiSeq 4000 platform (Illumina, San Diego, CA, United States), and mapped to scaffolds using Burrows-Wheeler-Alignment Tool (BWA) (Version 0.7.15) 1 (Li and Durbin, 2009). Genome Analysis Toolkit (GATK) was used to select SNP and indel 2 (McKenna et al., 2010). Sulv 1 genome was selected as the reference genome in the GATK analysis (Yan et al., 2020). High-quality SNPs and Indel variations were obtained as the following steps. (a) Retaining concordant sites both identified by GATK and VCFtools were retained (Danecek et al., 2011). (b) Filtering out SNP with quality value below 30, removing SNPs with an average coverage depth < 8× and with minor allele frequency (MAF) less than 5%. (c) Deleting insertions and deletions (InDels) with length less than10 bp were deleted. A total of 3,607,508 SNPs were identified.
Genome-wide association study for seed width, seed length, and 100-seed weight Only the SNPs with MAF ≥ 0.05 and missing rate < 10% were used in GWAS (Pongpanich et al., 2010). The lines with more than 95% missing for trait were filtered out (Liaw and Wiener, 2002). SW, SL, and HSW, and the above SNP markers in 196 mungbean accessions were used to conduct GWAS using four different methods, including 3VmrMLM (Li et al., 2022b) via software IIIVmrMLM (Li et al., 2022a), EMMAX (Kang et al., 2010), GEMMA (Zhou and Stephens, 2012), and CMLM (Zhang et al., 2010). The probability threshold for significant QTNs was set at 1/m = 2.77e-07 (m = 3,607,508) for all the four methods (Xu et al., 2018;Zhang Y. M. et al., 2019;Zhang Y. M. et al., 2019), and the LOD score threshold for suggested QTNs was set at LOD ≥ 3.0 for 3VmrMLM (Li et al., 2022b). Heatmaps of the linkage disequilibrium was generated by LDheatmap package (Shin et al., 2006), haplotype analysis was conducted by LDheatmap package (Barrett et al., 2005). The averages for those traits measured in 2018 and 2020 were used in GWAS.

Candidate gene identification
Candidate genes for salt tolerance were mined in the follow steps. (a) All the genes between the 30 Kb around regions for each of the significantly QTN were mined, where the LD-value was about 20 Kb in mungbean, (b) mined the Arabidopsis, rice and soybean homologous genes of those candidate genes, which were reported related to seed developments, seed production, phytohormone signaling pathways and carbohydrate metabolism pathways, etc. (Li et al., 2019), as the candidate genes. (c) The selected genes showing different expression between two groups of mungbean accessions contrasting in seed size (large seed vs. small seed) (see below) were considered as candidate genes.
In the key candidate gene identification, the extracted RNA in two large-seeded accessions at 10 and 25 DAF were treated with RNase-free DNase I (Promega, Madison, WI, United States). After reverse transcription, the cDNA was used as a template for RT-qPCR using the Takara Bio TB Green Premix Ex Taq (Tli RNase H Plus). The detail progress was described by Liu et al. (2022b). Reactions were run on a Bio-Rad CFX96 system. EVM0007380 (homologous of At3g18780) was used as the CK in this experiment. Primers were designed by NCBI and tested by RCR of tubulin. The t-test was adopted in the hypothesis testing, P < 0.05, P < 0.01, and P < 0.001 indicated significant probability levels at 0.05, 0.01, and 0.001, respectively. Information of the primers used is presented in Supplementary Table 1.

Genome-wide association studies for seed-size-related traits in mungbean
Detection of main-effect quantitative trait nucleotides for seed-size-related traits in each environment After removing the SNPs with an average coverage depth < 8× and with a MAF less than 5%, we identified more than 3.6 million SNP markers. In the single-environment analysis, the phenotypic observations for each trait in 196 accessions measured in 2018 and 2020 were used to associate with 3,607,508 SNPs using 3VmrMLM, EMMAX, GEMMA, and CMLM under the situations of five subpopulations and polygenic background control (kinship matrix) (Supplementary Data Set 3). As more than 10,000 QTNs were identified by GEMMA for HSW in 2018, the relevant results were not used in the subsequent analysis. As a result, 208 significant QTNs were identified for the above traits. Thirteen significant QTNs were simultaneously identified in two environments by two GWAS methods (Supplementary Table 3; Supplementary Data Set 4), some significant QTNs are presented in Figure 2. For example, Chr10-25206533-25223155 The frequency distributions of seed-size-related traits. Frequency distributions of HSW (A) (g), SL (B) (mm), and SW (C) (mm) in 196 mungbean accessions, which were measured in 2018 (brown bar) and 2020 (black bar). SD, standard deviation. The associations of HSW with SW and SL, the average dates of those traits measured in 2018 and 2020 were used in the partial correlation analysis (D).

Detection of quantitative trait nucleotides for seed-size-related traits in multiple environments
To detect more stable QTNs, three seed-size-related traits of 196 mungbean accessions measured in 2018 and 2020 were used to associate with 3607508 SNP markers using twoenvironment 3VmrMLM joint analysis. As a result, 32, 33, and 18 significant QTNs were identified for HSW, SL, and SW, respectively (Supplementary Table 3), and had a 1.08% average proportion of total phenotypic variation explained by each QTN. Moreover, eight significant QTNs were identified (Supplementary Table 4). For example, Chr1-8161305-8347626 (LOD = 24.09∼36.33) and Chr10-25222572-25223133 loci (LOD = 29.75∼37.89) were detected to be associated with HSW and SL, respectively (Supplementary Tables 3, 4).
Based on all the above main-effect QTNs in singleand multiple-environment analysis, five stable QTNs across various methods and/or two environments were found (Supplementary Table 5 Detection of quantitative trait nucleotide-by-environment interactions for seed-size-related traits in multiple environments All the above datasets in GWAS were used to detect QEIs using 3VmrMLM. As a result, 5, 10, and 5 significant QEIs were The P1-values were calculated by CMLM, EMMA, and 3VmrMLM, The P2-values were calculated using paired t-test from the average FPKM values at three stages between two high seed weight (n1 = 2) and tow seed weight (n2 = 2) mungbeans, and their significances were marked by * (0.05 level); FC and NA represent fold change and no expression, respectively. found to be associated with HSW, SL, and SW, respectively (Supplementary Figure 2; Table 2). Among these QEIs, 5 had zero dominant-by-environment interaction effects, and 7 had zero additive-by-environment interaction effects. For example, the two loci Chr4-26262890 and Chr4-31677341 for HSW had only additive-by-environment interaction effects of 0.12 ( Supplementary Figures 2A-C In addition, among these QEIs, the QEI locus Chr9-24007163 for SW had large effect, and r 2 was 5.8% (Supplementary Figure 2B, LOD = 61.95). The additive and dominant effects in environment 1 were −0.14 and −0.098, respectively.
Using the same approach described above, among 65 genes around 20 QEIs, four were homologous to previously reported seed development related genes in rice and A. thaliana ( Table 2), although new experiments are necessary to explore these novel GEI-trait associations. The four genes were described as below. VrFATB was linked to the locus Chr4-30176682 (Supplementary Figure 2A). As described in Bonaventure et al. (2003) and Sun et al. (2014), FATB is  The P-values were calculated using paired t-test from the average RPKM values at three stages between two high seed weight (n1 = 2) and tow seed weight (n2 = 2) mungbeans, and their significances were marked by * (0.05 level); FC and NA represent fold change and no expression, respectively.

FIGURE 3
The expression of eight key candidate genes. The expression profiling of eight key candidate genes significantly associated with seed-size-related traits. The expression profiling of eight key candidate genes between two high-seed-weight and two low-seed-weight mungbeans (A). Real-time PCR analysis of the eight key candidate genes; the t-test was used to test the significant differences of genes expression between two high-seed-weight mungbeans at 10 DAF and 25 DAF (B). DAF, days after flowering.
a major determinant of saturated fatty-acid synthesis, and increases FATB activity at low temperature during seedling establishment caused high saturated fatty-acid content in plant. VrGSO1 was linked to the locus Chr4-42563100 (Supplementary Figure 2A). As observed in Creff et al. (2019), GSO1 was a stress signal-pathway-related gene, and stress-associated MPK6 protein acted downstream of GSO1 in developing embryo. VrPAT14 was linked to the locus Chr9-24007163 (Supplementary Figure 2B). In Zhao et al. (2016), PAT14 was involved with NPR1-dependent salicylicacid signaling. VrLACS2 was linked to the locus Chr8-32848165 (Supplementary Figure 2C), in which VrLACS2 was essential for normal cuticle development in Arabidopsis (Schnurr et al., 2004) and CrLACS2 suppression resulted in 50% less oil, yet with a higher amount of chloroplast lipids under N-deprivation (Bai et al., 2022).

Haplotype analysis of the main candidate genes
Two DEGs, VrEmp24/25 and VrKIX8, were detected in the single-and multi-environment analyses (Figures 4A,B), and verified by RT-qPCR. Their haplotypic analyses were described as below.
In the haplotype analysis of VrEmp24/25, five SNP markers were found to be within VrEmp24/25 and the promoter region (Supplementary Data Set 8), and the two SNP markers in VrEmp24/25 were used to consist of three haplotypes ( Figure 4D). Among the three haplotypes, hap 1 (5.17 g) had significantly higher HSW than hap 2 (1.58 g) and hap 3 (4.50 g; P = 2.11E-29) (Supplementary Table 7). Thus, hap 1 is elite haplotype. And the elite haplotypes TT made up more than 90.9% (160/176) in the cultivated mungbeans. VrEmp24/25 with elite haplotype frequencies less than 45% in wild mungbeans (Supplementary Table 7; Figure 4) can be exploited for the improvement of mungbean cultivars.
Based on these results, we deduced that these two SNP and six SNP cause the difference expression of the VrEmp24/25 and VrKIX8 gene, respectively. The discovery of VrEmp24/25 and VrKIX8 two domestication/improvement genes can accelerate breeding selections and facilitate ideal crop designs.

Expression patterns of seed development pathway genes in mungbean
As seed development pathway genes were largely unknown in mungbean, we mined seed development pathway genes by comparative genomics and transcriptomics analysis. As a result, 54 genes in seed-development pathway were identified in this study (Figure 6; Supplementary Data Set 9), such as two GPA1, one AGB, and one AGG3. In the ubiquitin proteasome pathways, two DA1, one DA2, one SOD2, one EOD1, and one UBP15 rather than SAMBA were identified. In the auxin pathways, two ABA2, one ABI5, three SHB1, five IKU2, and three CKX2 rather than IKU1 and MINI3 were identified (Figure 6A). Five transcription factors including three BES1, and two SOD7 were identified. Moreover, 16 genes for seed size developments were found to be with uncertain pathways, including three KIX8, five MES1, and one KLU (Figure 6A; Supplementary Data Set 9). Among the 54 genes, 13 genes were significantly differentially expressed (P-value < 0.05, t-test) between two low-seed-weight (nos. G169 and G171) and two high-seed-weight (no. G141 and G143) accessions in the 196 mungbean accessions using the transcriptome data at 10, 15, and 25 DAF (Figure 6B; Supplementary Data Set 8). Moreover, almost 90% of the 54 genes (48/54) had higher expressions in the early stage of seed development (10 and 15 DAF) than in the late maturation stage (25 DAF), including VrKIX8 (EVM0032114), which was commonly identified in the GWAS by 3VmrMLM for HSW and SL. And EVM0010067/VrABA2, EVM0033315/VrSHB1, EVM0028440/VrABI5, and EVM0030447/VrIKU2 were also identified in the GWAS by 3VmrMLM, within 100 Kb region of significant QTNs (Table 1).
We also did the PPI analysis among the seed development pathway genes, and found five pairs of PPIs were larger than the medium confidence value of 0.40 (Supplementary  Table 7), indicating the existence of significant PPIs, i.e., EVM0013794.1 (VrAGG3) and EVM0006667.1 (VrDA2) Genetic analysis of VrKIX8. Local Manhattan plots for HSW in multi-environments. LOD ≥ 3.0 for the 3VmrMLM as the significant QTN (A). LD heatmaps surrounding Chr1-8161305-8347626 loci (B). Genes around the significant QTN region, shown at the bottom (C). Haplotype analysis of VrKIX8 (D). Wil, the wild accessions. Cul, the cultivated accessions. symbol "**" means omit the same sequence part.

Discussion
The high-yield and efficiency breeding progress of mungbeans have been limited by the lack of ideal yieldrelated genes. At present, few QTNs or QTLs of yield-related traits in mungbeans have been reported (Kang et al., 2014). This study provided a genetic analysis of seed-size-related traits in mungbeans, to improve the accuracy of significant QTNs, we used multiple genome-wide M0017 association studies combined with multi-omics analysis to mine candidate genes associated with yield-related traits. Firstly, a total of 98 QTNs and 20 QEIs were identified using 3VmrMLM, while 95 and 15 QTNs were identified using EMMAX, and CMLM, respectively. Then, in the identification of candidate genes, 12 key candidate genes were mined, and seven of them including VrKIX8, VrEmp24/25, and VrPAT14 were evidenced by transcriptome analysis and RT-qPCR analysis. Lastly, through haplotype analysis, the thirtieth amino acid of VrEmp24/25 in the elite haplotype was changed from Ile to Phe. And there were six SNP in the promoter and 5 UTRs of VrKIX8, however, the amino acid sequence of VrKIX8 in the elite haplotype was not changed. The results provided the theoretical basis for both the functional identification of seed-size-related genes and for quality improvements in mungbean breeding.
Multiple genome-wide association studies methods combined with multi-omics analysis in mining candidate genes In the GWAS, how to identify candidate genes around significant QTNs has been a challenge. Liu et al. (2020c), Zhang et al. (2021), and Gong et al. (2022) selected the 100kb interval upstream and downstream of the significant QTN as the candidate interval in watermelon and soybeans. Usually, the interval has been chosen according to the LD decay values.
Notably, 3VmrMLM showed more powerful ability in the detection of significant QTN than GEMMA, EMMAX, and CMLM, as it found more differentially expressed key candidate genes than other methods. The combination of 3VmrMLM and multi-omics analysis in the genetic analysis of complex traits was helpful.
Besides, in maize, the loss function of EMP24 and Emp25 would impair embryo and endosperm development (Xiu et al., 2020). EMP24 was required for the splicing of nad4 (Ren et al., 2019), and the lack of either Nad4 or Nad5 blocked the assembly of complex I holoenzyme in Arabidopsis (Ligas et al., 2019). The loss of the steady-state level of mitochondrial nad5 mature mRNA blocked the assembly of complex I and caused an arrest in endosperm development . Lastly, the elite haplotypes of VrEmp24/25 (TT) made up the main proportion of more than 90.9% in cultivated mungbeans, 45% in wild mungbeans ( Figure 4E). The HSW in hap 1 haplotypes accessions was significantly higher than that in hap 2 and hap 3 (P = 2.11E-29). It was reported that a single amino acid completely prevented the appearance of the enzyme in the medium, and we inferred that the related variation could lead to the change in enzyme activity (East et al., 1990;Alfson et al., 2018).
There have four evidences to take VrKIX8 as another important seed-size trait gene. Firstly, VrKIX8 associated with Chr1-8161305-8347626 (LOD = 24.09∼36.33) for HSW and SL were detected in multi-environment by 3VmrMLM ( Figure 5A; Supplementary Table 5). Secondly, VrKIX8 (LOD = 24.09∼36.33) had significantly differentially expressed between high-and low-HSW accessions ( Figure 3A). Then, in Arabidopsis, the disruption of KIX8/9 and PPD1/2 could cause large seeds due to increased cell proliferation and cell elongation in the integuments (Liu et al., 2020a). In soybeans, the loss of the function GmKIX8-1 showed a significant increase in the size of seeds and leaves. In addition, the increase in organ size was due to the increased cell proliferation, rather than cell expansion. GmKIX8-1 showed negatively regulated cell proliferation in plants (Nguyen et al., 2021). Lastly, the elite haplotypes of VrKIX8 (ATCGAA) made up the main proportion of more than 73% in cultivated mungbeans, 40% in wild mungbeans. Moreover, there are four SNPs in the promoter and of VrKIX8, and two SNPs in the CDS region, however the amino acid sequence did not change between the elite haplotypes and the other haplotypes ( Figure 5C). The HSW in hap 1 haplotypes accessions was higher than that in hap 2 to hap 4 (P = 1.19E-21). We supposed that the mutations may have influenced the translation efficiency of VrKIX8 and caused low expression in cultivated accessions during mungbean domestication.
We also analyze the possible correlation between the main seed development pathways. Among the 54 genes, five genes (VrAGG, VrDA2, VrAGB, VrGPA1-1, and VrGPA1-2) consisted of five pairs of significant PPIs. Interestingly, four pairs PPIs were found to be in the G-protein-signaling pathway, and one pair of PPIs was found to be in the G-protein-signaling and the ubiquitin proteasome pathways (Figure 6; Supplementary Table 6). Ubiquitin proteasome pathway is an important pathway for the selective degradation of proteins and seed development (Smalle and Vierstra, 2004), and the G-protein-signaling pathway is a ubiquitous cell transmembrane signal transduction pathway in eukaryotes (Huang et al., 2006). Moreover, mutations in GPA1 or AGB1 could cause short flowers (Lease et al., 2001;Ullah et al., 2001). The overexpression of AGG3 promoted seed and organ growth by increasing cell proliferation, and loss-offunction mutations in AGG3 caused small seeds and organs (Chakravorty et al., 2011;Li et al., 2012). The ubiquitin receptor DA1 could control seed size by restricting cell proliferation in maternal integuments (Li et al., 2008). DA1 functioned synergistically with DA2 to restrict seed growth, and DA2 physically interacted with DA1 in vitro and in vivo (Song et al., 2007;Xia et al., 2013). This interaction could mediate the interactions between the G-protein-signaling pathway and the ubiquitin proteasome pathway, which might offer an important clue in the mechanism analysis of seed development.
In addition, 48 genes had higher expressions in the early stage of seed development than in the late maturation stage of seed development, indicating that seed-development-related genes function primarily in the early stages of seed development, which was consistent with the findings of Zuo et al. (2022) in soybean.

Conclusion
This study conducted GWAS for seed-size-related traits in mungbeans. 98 QTNs and 20 QEIs were identified using 3VmrMLM, while 95, >10,000, and 15 QTNs were identified using EMMAX, GEMMA, and CMLM, respectively. A total of 12 key candidate genes were mined, which were homologous to known seed-development genes in rice and A. thaliana. VrEmp24/25 and VrKIX8 were identified as main candidate genes around two stable QTNs, the two candidate genes were further confirmed by RT-qPCR and haplotype analysis, and prevalent haplotypes of VrEmp24/25 and VrKIX8 may be useful in mungbean breeding.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: The WGS sequencing data of 196 mungbean accessions was uploaded to NGDC, with subCRA011538, subSAM100395, and PRJCA010704 ID.


s2fieldsofstudy:
Agricultural And Food Sciences
Environmental Science

provenance: AgriculturalAndFoodSciences,2022-2022-Agricultural And Food Sciences-2022.jsonl.gz.gz

text:
Diversity of Microbial Communities of Pinus sylvestris var. mongolica at Spatial Scale

Soil microorganisms play an indispensable role in the forest ecosystem. It is necessary to study the soil microorganisms in Pinus sylvestris var. mongolica, which is one of the afforestation species widely planted in the northern sandy region of China. We collected soil samples of P. sylvestris at large spatial scales and analyzed bacterial and fungal community composition differences using high-throughput sequencing techniques. The results showed that: (1) the richness index of different sandy lands was significantly different. The α-diversity of bacteria was the highest in Mu Us Sandy Land, and the α-diversity of fungi was the highest in Horqin Sandy Land. (2) The dominant phyla of bacteria were Actinobacteria, Proteobacteria, Chloroflexi and Acidobacteria, while the dominant phyla of fungi were Ascomycota and Basidiomycota. The relative abundance of dominant phyla was different. (3) Temperature and precipitation were the main driving factors of bacterial and fungal community change at large spatial scale. In addition, bacteria were also affected by total nitrogen, soil organic carbon and pH content; fungal community was affected by pH. The microorganisms showed obvious differences in geographical distribution, which could provide ideas for promoting sustainable management of P. sylvestris stand.

Introduction
Soil microorganism: a general term for microscopic organisms living in soil that are indistinguishable from the naked eye. It is estimated that one gram of soil contains up to hundreds of millions of microorganisms [1], among which bacteria and fungi are the key organisms, given their great diversity and richness. It has been dubbed "microbial dark matter" due to its vast variety and unknown function [2]. The high heterogeneity of the soil environment, which is an important hub connecting various spheres and ecosystems, breeds a great diversity of microorganisms; it is also the engine of the soil element cycle and maintains the sustainable development of human beings and the Earth's ecosystem [3,4]. Soil microorganisms play an indispensable role in the forest ecosystem [5,6]. It is the bridge between aboveground vegetation and the underground soil system [7,8], participating in carbon and nitrogen cycles, regulating the decomposition of organic matter and energy flow [9][10][11][12][13][14], and affecting the succession of plant community structure [8]. The sensitive response of microbial diversity to the surrounding environment reflects the stability and function of an ecosystem to some extent and can be an important indicator of soil change [15].
The diversity and composition of microbial communities show great spatial variability and are associated with changes in some biological or abiotic factors [16][17][18]. Above-ground vegetation provides nutrients needed for the survival of soil microorganisms in the process of material exchange with the external environment and changes the composition of the microbial community by affecting soil properties through litter [19] and root exudates [20,21]. The complex environmental characteristics of soil and geographic factors play a prominent role in determining the soil microbial community structure [22]. There is no conclusion that soil pH has a significant impact on microbial diversity. Some previous studies have described a strong influence of pH on the composition of microbial communities over different scales. Bacterial community structure changes significantly even in a very small variation of pH values [23][24][25]. Fierer and Jackson [26] and Griffiths et al. [27] have demonstrated that pH is a primary factor in a continent-wide survey of microbial communities; other studies have shown that pH has nothing to do with bacterial community structure [28]. There is a coevolution between plants and microorganisms. Plants provide nutrients for microorganisms, and microorganisms transform litter into inorganic nutrients needed for plant growth through decomposition. Rhizosphere nutrients and bioavailability directly or indirectly effect the microbiome assembly and functions [29]. Banerjee et al. [30] found that soil moisture effected microbial community structure more than nutrients. Similarly, the influence of geographical factors and climatic factors, such as temperature, altitude and precipitation on soil microbial diversity, has not been determined [31][32][33].
Pinus sylvestris var. mongolica (P. sylvestris for short below), a geographical variety species of the Scots pine, is the widely planted arbor for constructing protective plantations in "Three North" (Northwest, North and Northeast of China) [34,35], on account of its fast growth and strong adaptability to poor environmental conditions. Up until now, most of the studies on P. sylvestris have focused on the diversity of fungi or ectomycorrhizal fungi at a particular location, the change of fungal community with a shift ingrowth stage or the microbial difference between plantation and natural forest and have ignored understanding the microbial community assembly across spatial scales.
In this study, we collected soil of P. sylvestris from different areas and studied the community structure of bacteria and fungi, in order to provide some basis for solving the decline phenomenon of P. sylvestris in introduced areas. High-throughput sequencing technology was applied to study the composition of the microbial community. To some extent, the composition of microbial communities in different introduced areas reflects the preference of their respective environments. We attempted to explore the response of microorganisms to environmental factors at a large spatial scale, which helped to clarify the reasons for the decline of P. sylvestris in the introduction area, providing some theoretical support for the introduction of the P. sylvestris plantation. We hypothesized that (1) microbial diversity was different in different locations and that the diversity of a healthy sand was higher than that of declining area; (2) the composition of microbial community was more sensitive to climatic factors than soil factors.

Site Information
We conducted this study in four typical areas of P. sylvestris forest (Figure 1). Specific site information is shown in the Table 1. Honghuaerji Forest Park (Inner Mongolia Autonomous Region), located in the southeast of Hulunbuir (HB) Sandy Land, with a sandy P. sylvestris system and a grassland wetland landscape, is at an elevation of 500-1000 m, with a temperate continental monsoon climate, characterized by rain and heat in the same season [36]. Beautiful mountains and rivers, rich flora and fauna make it a local tourist hotspot.
Saihanba Mechanical Forest Farm (Hebei Province), located in the southern edge of Otindag (OD) Sandy Land, is at an elevation of about 1510-1939 m, with a cold temperature continental monsoon climate, with long winters and low temperatures [37]. The forest coverage rate has been greatly improved by artificial afforestation and ecological management in the past 60 years. There are vast forests and grasslands, clear plateau lakes and historical relics of the Qing Dynasty. Manchu and Mongolian people have lived here for a long time, and their national cultures blend with each other, giving birth to rich folk customs. Its unique geographical location breeds abundant animal and plant resources, as well as natural hydrological resources and cultural attractions, which bring it unique charm.
Yulin Sandy Forest Park (Shaanxi Province), located in the southern edge of Mu Us (MU) Sandy Land, is at an elevation of about 1120 m, with a semi-arid continental monsoon climate. Its annual precipitation rate varies greatly; there are more rainstorms in the summer, adequate light and heat which are sufficient for plant growth, and large temperature differences between day and night [38]. The sandy land was seriously eroded by wind and sand at first but now, the forest coverage rate has been significantly improved by ecological management; the cultural landscape and tourism industry have been developed.
Zhanggutai Sandy Forest Park (Liaoning Province), located in the southern margin of Horqin (HQ) Sandy Land, is at an elevation of 220-240 m, with a temperate continental monsoon climate, characterized by simultaneous humidity and heat, and little precipitation, mainly concentrated in the summer [39]. It is an ecological landscape combining desert wonders and artificial forests. The reserves of silica sand and peat are abundant; the seedlings of P. sylvestris have been cultivated for several decades. Honghuaerji Forest Park (Inner Mongolia Autonomous Region), located in the southeast of Hulunbuir (HB) Sandy Land, with a sandy P. sylvestris system and a grassland wetland landscape, is at an elevation of 500-1000 m, with a temperate continental monsoon climate, characterized by rain and heat in the same season [36]. Beautiful mountains and rivers, rich flora and fauna make it a local tourist hotspot.
Saihanba Mechanical Forest Farm (Hebei Province), located in the southern edge of Otindag (OD) Sandy Land, is at an elevation of about 1510-1939 m, with a cold temperature continental monsoon climate, with long winters and low temperatures [37]. The forest coverage rate has been greatly improved by artificial afforestation and ecological management in the past 60 years. There are vast forests and grasslands, clear plateau lakes and historical relics of the Qing Dynasty. Manchu and Mongolian people have lived here for a long time, and their national cultures blend with each other, giving birth to rich folk customs. Its unique geographical location breeds abundant animal and plant resources, as well as natural hydrological resources and cultural attractions, which bring it unique charm.
Yulin Sandy Forest Park (Shaanxi Province), located in the southern edge of Mu Us (MU) Sandy Land, is at an elevation of about 1120 m, with a semi-arid continental monsoon climate. Its annual precipitation rate varies greatly; there are more rainstorms in the summer, adequate light and heat which are sufficient for plant growth, and large temperature differences between day and night [38]. The sandy land was seriously eroded by wind and sand at first but now, the forest coverage rate has been significantly improved by ecological management; the cultural landscape and tourism industry have been developed.
Zhanggutai Sandy Forest Park (Liaoning Province), located in the southern margin of Horqin (HQ) Sandy Land, is at an elevation of 220-240 m, with a temperate continental monsoon climate, characterized by simultaneous humidity and heat, and little precipitation, mainly concentrated in the summer [39]. It is an ecological landscape combining desert wonders and artificial forests. The reserves of silica sand and peat are abundant; the seedlings of P. sylvestris have been cultivated for several decades.

Soil Sampling
Experimental samples were collected in April-May 2021; plantations of two sampling sites in each test location were selected for study; three 20 × 20 m 2 plots were taken from each place as repeated groups. P. sylvestris soil samples of 0-20 cm layer collected with a soil drill were pooled into sealed bags; stones, herbs roots, and litters were excluded with Microorganisms 2022, 10, 371 4 of 14 2 mm mesh [40]. The plots were separated by more than 100 m. A five-point sampling method was adopted, with three repeated soil samples collected in each point. All the composite samples (4 sandy lands × 2 places × 3 × 5 = 120 composite samples) were stored in dry ice for transport to the lab for analysis.

Environmental Variables
The soil samples were divided into two parts: one was naturally air-dried for the determination of soil physical and chemical properties; the other was used for high-throughput sequencing of soil microorganisms.
For the determination of soil properties, we referred to the method of Bao [41]. The soil pH value was obtained by measuring the suspension with a soil-water ratio of 1:2.5 using a pH meter; soil organic carbon (SOC) was determined by titration of ferrous sulfate heated by potassium dichromate; total nitrogen (TN) was determined by using the Kelvin method, with Se, CuSO 4 , and K 2 SO 4 as catalysts, adding concentrated sulfuric acid digestion, with 1 g of soil. MADAC (Anti-molybdenum antimony colorimetry) was used to estimate the phosphorus content in samples from soil. The determination of cation exchange capacity (CEC) was determined by ammonium acetate exchange using the Kjeldahl method.
Climate data for each sampling location were obtained from the China Meteorological Data Service Center (CMDC, http://data.cma.cn/en; accessed on 15 July 2021), including the mean monthly temperature (Ta), mean monthly precipitation (Pa), monthly average maximum and minimum temperature (T+ and T−).

Soil DNA Extraction and Amplification
Total genomic DNA were extracted using the DNeasy Power Soil Pro Kit (QIAGEN, Frankfurt, Germany) and stored at −20 • C prior to further analysis. A NanoDrop NC2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) was used to quantify DNA; the quality of DNA was detected by 1.2% agarose gel electrophoresis. For PCR amplification of the bacterial 16S rRNA genes at the V3-V4 region, the primers were 338F (5 -ACTCCTACGGGAGGCAGCA-3 ) and 806R (5 -GGACTACHVGGGTWTCTAAT-3 ) [42]. For the ITS1 regions of fungi, the primers were ITS1F (5-GGAAGTAAAAGTCGTAACAAGG-3) and ITS2 (5-TCCTCCGCTTATTGATATGC-3) [43]. The PCR amplification was performed by Applied Biosystems 2720 Thermal Cycler (Thermo Fisher Scientific, Waltham, MA, USA); the Quant-iT PicoGreen dsDNA assay was performed by a Microplate Reader FLx800 (BioTek, Burlington, VT, USA). PCR amplicons were purified and recovered by adding Vazyme VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China) and quantified with the fluorescent reagent of Quant-iT PicoGreen dsDNA assay kit (Invitrogen, Carlsbad, CA, USA); then, amplicons were mixed in proportion to the sequencing amount, and pair-end 2 × 250 bp sequencing was performed using the Illlumina MiSeq platform with a NovaSeq 6000 SP reagent kit at Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China).

Sequence Analysis
Microbiome bioinformatics were processed using QIIME2 2020.2 with slight modification according to the official tutorials (https://docs.qiime2.org/2020.11/tutorials/; accessed on 27 June 2021). Briefly, a Demux plugin was used for the sample splitting of mixed sequences and a Cutadapt plugin was used to cut primers [44]. Then, a DADA2 plugin was used for quality filtering, denoising and chimera removal of the sequences. [45]. After quality screening, non-singleton amplicon sequence variants (ASVs) were aligned with mafft and used to construct a phylogeny with fasttree2 [46,47]. The taxonomic assignment of the sequences was determined based on the bacterial SILVA reference database for 16S rRNA and the fungal UNITE reference database for ITS. Raw sequencing data have been submitted to the NCBI database (accession number: PRJNA785088).

Statistical Analysis
Sequence data analyses were mainly performed using QIIME2 and R packages (version 4.0.3). Microbial alpha diversity indices at the ASV level, including the Shannon diversity index, the Chao1 richness estimator, and Pielou's evenness were calculated using the ASV number in QIIME2, and visualized as box plots. Under the condition that the test data meet the homogeneity of variance, one-way analysis of variance (ANOVA) was conducted to compare the soil properties, with p < 0.05 considered statistically significant, supplemented by post hoc Tukey's tests. SPSS version 23.0 software (SPSS Inc., Chicago, IL, USA) was used for the analysis of variance. Non-metric multidimensional scaling (NMDS) was performed to test for differences in microbial communities in different locations based on unweighted UniFrac distance. Redundancy analysis (RDA), canonical correlation analysis (CCA) and Mantel tests were used to explain the relationship between the environmental factors and microbial composition. NMDS, RDA, CCA, and Mantel tests were performed in R.

Soil Properties
The soil nutrient content in HB Sandy Land and OD Sandy Land was significantly higher than that in HQ Sandy Land and MU Sandy Land. The maximum value of TN, AP and CEC was found in HB Sandy Land; the maximum value of SOC was found in OD Sandy Land, the minimum value of TN, SOC was found in MU Sandy Land, and the minimum value of AP and CEC was found in HQ Sandy Land. The soils of MU Sandy Land were alkaline compared to the acidic soils of other sandy lands (Table 2).

Diversity of Bacterial and Fungal Communities
A total of 4,684,298 high-quality bacterial sequences were obtained from the samples; the sequences were grouped into 6730 ASVs. A total of 8,927,828 high-quality fungal sequences were grouped into 3687 ASVs.
For bacteria, diversity indices based on identified ASVs (Figure 2) showed that the diversity and richness indices at MU possessed the highest values, while the lowest values were found at OD. The Shannon index has a certain difference between MU and OD; the Pielou index showed a certain difference between HB and HQ. The Chao1 indices measurements differed significantly between MU and other locations (p < 0.05).
For fungi, the diversity and richness indices at HQ showed the highest values, while the lowest values were found at HB. The Shannon index possesses a certain difference between HB and HQ. The Chao1 indices measurements differed significantly between HQ and other locations (p < 0.05). There was no significant difference in the Pielou index of different sandy lands.
Pielou index showed a certain difference between HB and HQ. The Chao1 indices measurements differed significantly between MU and other locations (p < 0.05).
For fungi, the diversity and richness indices at HQ showed the highest values, while the lowest values were found at HB. The Shannon index possesses a certain difference between HB and HQ. The Chao1 indices measurements differed significantly between HQ and other locations (p < 0.05). There was no significant difference in the Pielou index of different sandy lands.

Bacterial and Fungal Community Composition
The relative abundances were calculated on the ASVs numbers ( Figure 3).

Bacterial and Fungal Community Composition
The relative abundances were calculated on the ASVs numbers ( Figure 3). In the composition of the bacterial community, a total of 25 phyla were detected. There were six dominant bacterial phyla with relative abundance that reached more than 1%, Actinobacteriota was the most abundant phylum (49.88%), followed by Proteobacteria (21.15%), and then Chloroflexi (10.23%), Acidobacteriota (7.18%), Firmicutes (4.75%), and Gemmatimonadota (2.31%). The bacterial community composition differed among test locations. Dadabacteria and Bdellovibrionota were only detected in the HQ Sandy Land, Armatimonadota, p_GAL15 and p_MBNT15 were only detected in the MU Sandy Land, Elusimicrobiota were absent in the OD Sandy Land and HQ Sandy Land.
(a)    In the composition of the bacterial community, a total of 25 phyla were detected. There were six dominant bacterial phyla with relative abundance that reached more than 1%, Actinobacteriota was the most abundant phylum (49.88%), followed by Proteobacteria (21.15%), and then Chloroflexi (10.23%), Acidobacteriota (7.18%), Firmicutes (4.75%), and Gemmatimonadota (2.31%). The bacterial community composition differed among test locations. Dadabacteria and Bdellovibrionota were only detected in the HQ Sandy Land, Armatimonadota, p_GAL15 and p_MBNT15 were only detected in the MU Sandy Land, Elusimicrobiota were absent in the OD Sandy Land and HQ Sandy Land.
Soil bacterial and fungal community structure with the different sandy land was analyzed using NMDS based on the unweighted UniFrac distance. NMDS analysis showed clear area separation in both bacteria and fungi ( Figure 4). ANOSIM agreed with the NMDS in that sampling location caused differences in bacterial and fungal community structure (r = 0.951, p = 0.001 and r = 0.963, p = 0.001, respectively).

The Response of Bacterial and Fungal Communities to Environmental Variation
Shannon diversity, Chao1 richness and Pielou evenness indices of bacteria significantly correlated with AP and pH ( Table 3). The most significant factors associated with the diversity index were soil physical and chemical properties and T-and Ta (p < 0.05). The Chao1 index was significantly correlated with all soil physical and chemical factors, except CEC. Only the Shannon index showed a certain correlation with meteorological data.

The Response of Bacterial and Fungal Communities to Environmental Variation
Shannon diversity, Chao1 richness and Pielou evenness indices of bacteria significantly correlated with AP and pH ( Table 3). The most significant factors associated with the diversity index were soil physical and chemical properties and T− and Ta (p < 0.05). The Chao1 index was significantly correlated with all soil physical and chemical factors, except CEC. Only the Shannon index showed a certain correlation with meteorological data. For fungi, the Shannon index was significantly correlated with the AP, pH, and CEC ( Table 4). The Pielou index was significantly correlated with all soil physical and chemical factors, except pH. The factors that were most significantly associated with the Shannon index included AP, pH and CEC. The only factor that affected the Chao1 index was pH. To evaluate the relationship of microbial community structure with environmental variables at spatial scales, RDA analysis of bacteria, CCA analysis of fungal community and the Mantel test were implemented, according to the results of the DCA analysis. The results showed that the first two axes explained 79.50% and 49.38% of the total variation in bacteria and fungi, respectively ( Figure 5).
Based on the Mantel test, climate factor was remarkably correlated with bacterial (Table 5) and fungal (Table 6) community structure. The main factors that affected the bacterial community of soil were TN, SOC and pH. For the fungal community, pH showed a relatively obvious effect compared to other soil factors (r = 0.081, p = 0.018). The difference between MU Sandy Land in bacterial and fungal community and the other locations was caused by the pH. Pa was the main factor that explained the differences in community composition between HQ Sandy Land and others.  To evaluate the relationship of microbial community structure with environmental variables at spatial scales, RDA analysis of bacteria, CCA analysis of fungal community and the Mantel test were implemented, according to the results of the DCA analysis. The results showed that the first two axes explained 79.50% and 49.38% of the total variation in bacteria and fungi, respectively ( Figure 5). Based on the Mantel test, climate factor was remarkably correlated with bacterial (Table 5) and fungal (Table 6) community structure. The main factors that affected the bacterial community of soil were TN, SOC and pH. For the fungal community, pH showed a relatively obvious effect compared to other soil factors (r = 0.081, p = 0.018). The difference between MU Sandy Land in bacterial and fungal community and the other locations was caused by the pH. Pa was the main factor that explained the differences in community composition between HQ Sandy Land and others.

Variation in Bacterial and Fungal Composition
The differences in microbial diversity are mainly in the richness index. The bacterial Chao1 richness index in MU Sandy Land was significantly higher than other locations, and the fungal Chao1 richness index in HQ Sandy Land was significantly higher than other areas. This was not consistent with our hypothesis that the diversity of healthy sand was higher than degraded sand. Guo et al. [43] studied the ectomycorrhizal community structure of P. sylvestris at different scales and showed the same results, i.e., that the level of diversity was not correlated with the health status of the stand. Microbial diversity is mainly influenced by tree species [48,49]; no definitive study has shown that the diversity of healthy vegetation is higher. There was no significant difference in the evenness index except for the bacterial community in HB Sandy Land and HQ Sandy Land, indicating that the distribution of bacterial community was relatively uniform in all locations.

Variation in Bacterial and Fungal Community Structure
We found that Actinobacteria, Proteobacteria, Chloroflexi and Acidobacteria were the dominant taxonomic phyla in the bacterial community composition-a finding in accordance with those of Zhao [50] and Kang et al. [51]. Actinobacteria, involved in the complex of soil organic matter decomposition [52,53], such as cellulose and lignin, which explains why the relative abundance was highest in all locations. In addition, it also possesses a strong tolerance to drought and can produce substances that stimulate and promote plant growth [54]. Proteobacteria are considered eutrophic, with higher relative abundance in soils with high available carbon content [55], which is consistent with the result that the relative abundance of Proteobacteria were HB (25.83%) > OD (20.83%) > HQ (20.20%) > MU (17.33%). However, in other studies on the soil bacteria of P. sylvestris, Chloroflexi does not belong to the dominant phylum [42,56]. Previous studies showed that Chloroflexi also participates in the cycle of carbon and nitrogen except, for fixing CO 2 [57]. Acidobacteria is a group of acidophilic bacteria, which also effect the degradation of cellulose [58]. Pankratov et al. [59] found that, although the degradation function of Acidobacteria was not as good as other known cellulose-degrading bacteria, it possessed strong stress resistance and could survive in cold northern acidic wetlands, playing an important role in cellulose degradation under such conditions. As far as fungi are concerned, Ascomycota and Basidiomycota were the main dominant groups, which has been confirmed in many studies [60][61][62]. The relative abundance of Basidiomycota was higher than that of Ascomycota in all samples collected. Ascomycota are dominated by saprophytic fungi, which are involved in the decomposition of organic matter and nutrient cycling [63,64]. Basidiomycota can decompose lignocellulose and cause symbiosis with plants to form mycorrhiza [9,61,65]. Although the Ascomycota and Basidiomycota are the dominant phyla in the MU Sandy Land, the relative abundance of the Ascomycota (56.07%) in the MU Sandy Land is greater than that of the Basidiomycota (38.46%). Ascomycota was able to decompose complex organic matter and held a faster evolution rate, stronger drought tolerance and radiation tolerance, and, thus, is more suitable for survival in harsh environments [62].
To sum up, the relative abundance of bacterial and fungal taxa in soil depends on the spatial variation of soil environment and specific characteristics of sampling sites [66].

Response of Bacterial and Fungal Communities to Environmental Variation
The significant differences in the microbial community structure of P. sylvestris from the original locality to the introduced area resulted from the environmental heterogeneity of climate and soil factors. According to the Mantel test, meteorological factors have significant effects on bacterial and fungal community structures at large scale. Prior studies noted that soil properties, precipitation, and temperature constitute the key drives of soil microbial composition and metabolic potentials at large ecological scales [67][68][69]. In an Alaskan boreal forest, the response of microbial communities to warming was decreasing richness [70]. Zhou et al. [71] investigated the effect of temperature on the formation of microbial diversity at various scales and found that temperature showed a regulatory effect on microbial community diversity and microbial activities. In terrestrial environments, metabolic rates and ecosystem productivity generally increase as temperatures rise [72][73][74][75]. Ten years of irrigation in a water-limited pine forest increased tree growth and canopy cover, as well as litter fall and root biomass, and microbial communities that thrived as a result of the increased availability of resources [76]. Bacteria and fungi communities showed different adaptability to precipitation changes. Fungi showed stronger tolerance to water stress due to the existence of mycelia [77,78]. The response of microorganisms to temperature and precipitation at different ecological scales still needs further comprehensive research. The high heterogeneity of soil granular structure, porosity, water holding capacity and other aspects provides living space, sufficient nutrient conditions and a suitable microenvironment for various microorganisms, which is the basis of rich soil microbial diversity [79].
In this study, Mantel test results showed that the pH was closely related to both the soil bacterial and fungal community compositions. It has been suggested that pH is a major regulator of microbial diversity and community composition at the regional scale [80]. Soil pH regulates microbial growth associated with plant nutrients availability by controlling the morphology of compounds in the soil. Zhang et al. [81] proved that soil pH could change the microbial utilization efficiency of nutrients, physiological metabolic activity and competition between populations, directly or indirectly affecting microbial diversity. In addition, Mantel tests revealed that the total nitrogen and soil organic carbon content was a major factor influencing the bacterial community structure. Decomposition of litter, roots and dead biological cells is the main source of forest soil organic carbon, which is an important part of the carbon cycle [82]. Sul et al. [83] found that soil microbial community structure and specific group distribution were most affected by soil organic carbon. SOC was also shown to mediate the geographic distribution of the microbial communities in complicated and volatile environmental situations [84][85][86]. Total nitrogen content was the main factor affecting the EM fungal community composition [43]. The diversity and composition of microbial communities vary greatly at the spatial scale; the complex interaction mechanism among soil environment, meteorological factors and plants still needs to be further explored. In addition, we only analyzed the differences of community composition in a single period; thus, further long-term studies on the seasonal dynamics of microorganisms are needed.

Conclusions
In this study, the response of microorganisms to environmental changes on a large spatial scale was studied. The results showed that there were obvious regional differences in bacterial and fungal community composition and different sensitivities to environmental factors. Microbial community structure showed a close correlation with meteorological factors; pH had a significant effect on the microbial community. In addition, bacterial community was sensitive to the content of N and C. Contrary to our expectations, the microbial α diversity in healthy stands was not higher than that in declining stands; further study is needed to explore this phenomenon.
The degree of microbial community response to the change of environmental factors varies in different sandy lands. Understanding this correlational analysis could provide theoretical support to reduce the decline of P. sylvestris forest in the introduced areas.


s2fieldsofstudy:
Agricultural And Food Sciences
Biology

text:
Comprehensive analysis of cucumber C-repeat/dehydration-responsive element binding factor family genes and their potential roles in cold tolerance of cucumber

Background Cold stress is one of the main abiotic stresses limiting cucumber (Cucumis sativus L.) growth and production. C-repeat binding factor/Dehydration responsive element-binding 1 protein (CBF/DREB1), containing conserved APETALA2 (AP2) DNA binding domains and two characteristic sequences, are key signaling genes that can be rapidly induced and play vital roles in plant response to low temperature. However, the CBF family has not been systematically elucidated in cucumber, and the expression pattern of this family genes under cold stress remains unclear. Results In this study, three CsCBF family genes were identified in cucumber genome and their protein conserved domain, protein physicochemical properties, gene structure and phylogenetic analysis were further comprehensively analyzed. Subcellular localization showed that all three CsCBFs were localized in the nucleus. Cis-element analysis of the promoters indicated that CsCBFs might be involved in plant hormone response and abiotic stress response. Expression analysis showed that the three CsCBFs could be significantly induced by cold stress, salt and ABA. The overexpression of CsCBFs in cucumber seedlings enhanced the tolerance to cold stress, and importantly, the transcript levels of CsCOR genes were significantly upregulated in 35S:CsCBFs transgenic plants after cold stress treatment. Biochemical analyses ascertained that CsCBFs directly activated CsCOR genes expression by binding to its promoter, thereby enhancing plant resistance to cold stress. Conclusion This study provided a foundation for further research on the function of CsCBF genes in cold stress resistance and elucidating its mechanism. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03664-z.

Introduction
Low temperature is a major environmental factor affecting the growth and seasonal distribution of field crops and horticultural crops. When plants are exposed to freezing temperatures, their cold tolerance and freezing resistance are increased, that is, this adaptation process is called cold adaptation [1]. Many plants, such as Arabidopsis and oilseed rape, subjected to low temperature stress, have evolved a complex set of cold-adapted mechanisms involving gene transcriptional regulation and a wide range of physiological, biochemical and metabolic changes [2][3][4]. In recent years, the molecular mechanism of cold adaptation has been studied extensively and the Open Access *Correspondence: 374141645@qq.com; 163.hwx@163.com 1 School of Biological Science and Technology, University of Jinan, Jinan 250022, China 2 Shandong Institute of Pomology, Tai'an, Shandong 271000, China key regulatory factors of this complex network have been explored. In the Arabidopsis genome, 4-20% of the genes are regulated by cold signaling. Many transcription factors including MYB, NAC, MYC, AP2/ERF (APETALA2/ ethylene response factor) and bZIP can induce the expression of stress-related genes to protect cells from injury at low temperature [5][6][7]. Studies have shown that the cold signal regulation pathway dependent on CBF (CRT/DRE binding factor), namely ICE1 (Inducer of CBF/DREB1 expression 1)-CBF-COR (cold-responsive) pathway, is the most clearly studied cold response pathway in plants [4,8,9]. In this pathway, ICE1 is rapidly expressed after low temperature induction and promotes the expression of CBF genes and CBFs bind to the CRT/ DRE (C-repeat/Dehydration responsive element) cis-element on the COR gene promoters to activate COR gene expression. Thus, plant resistance to low temperature can be improved [10][11][12].
CBF transcription factors belong to the DREB subfamily of APETALA2/Ethylene-Responsive Factor (AP2/ERF) family, so all CBF proteins have a highly conserved AP2 DNA binding domain, and there are also two conserved characteristic sequences (PKK/RPAGRxKFxETRHP and DSAWR) on both sides of the AP2 domain [13]. The former sequence is directly located in the upstream of AP2 DNA binding domain and may be related to protein transport [13]. Arabidopsis thaliana contains six CBF genes, namely AtCBF1 (DREB1B), AtCBF2 (DREB1C) and AtCBF3 (DREB1A), AtCBF4, AtDDF1 and AtDDF2, among which AtCBF1, AtCBF2 and AtCBF3 are extensively induced by low temperature and can all improve the cold tolerance of Arabidopsis thaliana. These three AtCBF proteins have high homology and are arranged in tandem array as AtCBF1-AtCBF3-AtCBF2 on the short arm of chromosome 4 [14][15][16]. The expression of AtCBF4 is induced by drought and ABA, but not by low temperature, while overexpression of AtCBF4 gene can enhance drought tolerance and cold tolerance of Arabidopsis thaliana [17]. AtDDF1 and AtDDF2 are related to salt tolerance of Arabidopsis thaliana [18].
Currently, CBF genes have been isolated from many plant species, such as apple, soybean, rice, tomato, wheat, barley and maize, and have been shown to play important roles in plant cold response [19]. Overexpression of AtCBFs in other species can enhance the cold resistance of each species, and the heterogeneous expression of CBFs from other plants in Arabidopsis can also enhance the cold tolerance of Arabidopsis [20][21][22]. When TaCBF14 and TaCBF15 genes from wheat were transferred into barley, the frost tolerance of transgenic plants was higher than that of wild-type barley [23]. Overexpression of OsDREB1A could enhance cold resistance of rice [24]. Transgenic potato with AtCBF1 gene enhanced cold tolerance and induced physiological changes related to adaptation to cold environment [20]. Overexpression of CsGG3.2 could enhance cold resistance via positively regulated the expression of CBF genes in cucumber [25]. Moreover, exogenous melatonin could increase the expression of CBF1 and enhanced the cold tolerance of cucumber seedlings [26]. CBFs are important regulators of plant growth and low temperature response, and the biological function of CBFs in regulating cold tolerance is highly conserved among plants, but also species-specific [27].
Low temperature rapidly activates the expression of CBF in plants, which then binds specifically to the DRE/ CRT cis-element of the COR gene promoters to induce its expression [28,29]. Many COR genes have been isolated from plants and are known as KIN (cold-inducible), ERD (early engend-inducible), LTI (low temperature induced) and RD (response to dehydration). These genes include LTI78, COR78 (RD29A), COR47, COR15A and KIN1 [30][31][32][33]. In addition to Arabidopsis, many COR homologues have been cloned from other plants [34]. The expression of CsCOR1 gene in tea was significantly induced by low temperature and drought, and the heterologous expression of CsCOR1 in tobacco enhanced its salt tolerance and dehydration tolerance [35]. Under cold stress, two cucumber CBF-inducible COR genes, CsCOR15b and CsKIN1, were higher expressed in CsGG3.2 overexpression plants [25].
The discovery of CBF transcriptional activators provide a new way to improve plant cold resistance, and lay a theoretical foundation for further discovery of key genes in plant cold resistance mechanism, which has a wide application prospect and important application value in crop and vegetable quality improvement. Cucumber (Cucumis sativus L.) is an economically important crop cultivated worldwide [36]. The functions of CsCBFs have not been systematically identified in cucumber. In this study, three CsCBF genes were identified in cucumber and the comprehensive analyses including the gene structures, conserved domains, phylogenetic analysis and cis-elements in promoters were further performed. In addition, the expression patterns of CsCBFs under different abiotic stresses were analyzed. Furthermore, overexpressed CsCBFs cucumber seedlings increased their cold tolerance by activating CsCOR genes expression. These results provide a basis for further research on cold tolerance mechanism of cucumber.

Identification and analysis of CBF genes in Cucumber
To identify CsCBF family genes in cucumber genome, the six Arabidopsis CBF protein sequences and the AP2/ ERF conserved domain and characteristic motif of CBF protein were employed as queries to search against the cucumber genome database using BlastP programme, respectively. Finally, three cucumber CBF family genes were identified, and named CsCBF1 to CsCBF3 according to their sequence similarity and phylogenies with individual AtCBF proteins. Multi-sequence alignment (MSA) analysis of the three cucumber CBF protein sequences and six Arabidopsis CBF protein sequences showed that the AP2/ERF domain and the characteristic motifs of CBF family were highly conserved, and the C-terminal and N-terminal of CsCBF proteins were significantly different (Fig. 1). The information of the CsCBFs, including the gene ID, gene name, chromosomal locations, isoelectric points (pI), and amino acid length was shown in Table 1. These three CsCBF genes were mapped on chromosome 3, 5 and 5 of cucumber, respectively, and the isoelectric points of CsCBF1, CsCBF2, CsCBF3 protein were 5.10, 5.16 and 4.85, respectively ( Table 1). The isoelectric points of all three CsCBF proteins were less than 7, indicating that the three proteins were acidic.

Phylogenetic analysis, gene structure and conserved motif analysis of CBF genes in cucumber, tomato and Arabidopsis
To analyze orthologous or paralogous relationships of the CBF genes from nine plant species, a phylogenetic tree of these genes was constructed from amino acid sequences (Fig. 2). As shown in Fig. 2, the 70 CBF proteins could be roughly divided into 18 groups, of which CsCBF1 and CsCBF2 belonged to group 5 and CsCBF3 belonged to group 8 ( Fig. 2; Table S1). To better evaluate the evolutionary relationships of the CsCBF proteins, we further analyzed the phylogenetic tree, gene structures and conserved motifs of six AtCBFs, seven SlCBFs and three CsCBFs (Fig. 3). As shown in Fig. 3A, the resulting tree categorized these CBF proteins into three clades, designated CladeI, CladeII, and Clad-eIII, and cucumber CBF members were classified into Clade I and III. Phylogenetic analysis also revealed that CsCBF1 and CsCBF2 had the highest homology with SlCBF1, SlCBF2 and SlCBF3, and all clustered in Clade I (Fig. 3A). CsCBF3 was classified into CladeIII, which included four SlCBF and two AtCBF proteins, and Cla-deII only composed of four AtCBF proteins.
The gene structure of CBFs from cucumber, Arabidopsis and tomato were also analyzed, which was consistent with the results of phylogenetic analysis (Fig. 3B). The number of exons in AtCBF, SlCBF and CsCBF genes was conserved, ranging from one to two exons. We found that the two clades, CladeI and Cla-deII, had same gene structures, which all contained only one exon and no intron (Fig. 3B). The six proteins, AtDDF1, AtDDF2 and SlCBF4-SlCBF7, in CladeIII  contained one exon, while CsCBF3 in CladeIII had two exons and one intron (Fig. 3B).
To further analyze the structural diversity and predict the function of the CBF proteins, the motif analysis of them was carried out by MEME ( Fig. 3C; Fig. S1). A total of ten distinct motifs, named Motif1-Motif10, were identified (Fig. S1). Motif1, which was representative AP2/ERF domain, and motif 2 and 4, which were the characteristic domains of CBF family, were identified in all CBFs. Some of the specific motifs were absent in certain clades. For example, motif 6, 7 and 9 existed only in Clade II, but were absent in all the members of the CladeI Collectively, CBF proteins with close evolutionary relationships in the phylogenetic tree generally had similar gene structures and conserved motifs, indicating that evolution of each subfamily in the three different species was relatively conserved.

Cis-elements identification of CsCBF gene promoters in Cucumber
Previous studies have indicated that most of CBF proteins regulate plant growth and tolerance to various abiotic stresses [37]. To investigate the biological function of CsCBF genes in cucumber, the potential cis-elements were identified on the 2-kb promoter regions of the CsCBF genes by PlantCARE (Table S2). As shown in Fig. 4, cis-elements responding to hormones such as abscisic acid (ABA), salicylic acid (SA), jasmonate acid (MeJA), auxin and gibberellin (GA) were presented on the promoters of CsCBF genes. Moreover, the CsCBF promoters also contained cis-elements in response to abiotic stresses, such as low temperature, defense and stress. In addition, the three CsCBF gene promoters all had light response signal elements (Fig. 4), suggesting that CsCBF genes might also be involved in the regulation of cucumber growth and development by light signal. In conclusion, the three CsCBF genes in cucumber may be involved in the response to multiple plant hormones and abiotic stresses.

Expression patterns of CsCBFs under different abiotic stresses
To investigate whether CsCBF genes respond to abiotic stress in cucumber, qRT-PCR was used to detect the expression levels of CsCBF genes under different stress conditions including low temperature (4°C), salt (100 mM NaCl) and ABA (100 μΜ ABA) (Fig. 5). Under low temperature treatment, the expression patterns of CsCBF1 and CsCBF2 showed a similar trend, which was firstly rapidly increased, and reached maximum values at 3 h, decreasing thereafter (Fig. 5A). However, the expression levels of CsCBF3 reached the maximum at 9 h, suggesting that CsCBF3 was less sensitive to low temperature than CsCBF1 and CsCBF2. Similar to the expression pattern after low temperature treatment, the expression levels of CsCBF1 and CsCBF2 achieved the maximum after 3 h of salt stress, while the expression level of CsCBF3 was highest after just 1 h of NaCl treatment, then rapidly decreased to the initial level and remained at a low level all the time (Fig. 5B). Different from the above treatments, CsCBF genes responded to ABA more rapidly, and all significantly increased at 0.5 h. The expression levels of CsCBF2 and CsCBF3 increased rapidly at 0.5h and reached the maximum at 1 h, then dropped to the initial level at 3h and increased again at 6h, decreasing  B 100 mM NaCl treatment. C 100 μΜ ABA treatment. The extracted RNA samples were collected at 0.5h, 1h, 3h, 6h, 9h, 12h and 24h after the corresponding treatment, respectively. The cucumber β-actin gene was used as an internal control, and three biological replicates were used for gene expression analyses. Error bars were the standard errors (SE). Different lowercase letters represented significant differences (P < 0.05) thereafter (Fig. 5C). Unlike CsCBF2 and CsCBF3, the expression level of CsCBF1 reached its highest value at 6h of ABA treatment (Fig. 5C). These results suggested that CsCBFs were involved in response to low temperature, salt and ABA.

Subcellular Localization of CsCBFs
PredictProtein software was used to predict the subcellular localization of CsCBF1, CsCBF2 and CsCBF3 proteins [38], and the prediction result showed that CsCBF1, CsCBF2 and CsCBF3 proteins were all in the nucleus (Fig.  S2). To verify the predicted results of CsCBF proteins, the fusion protein vectors 35S:CsCBF1-GFP, 35S:CsCBF2-GFP and 35S:CsCBF3-GFP were constructed, respectively (Fig. 6A). The GFP emited a green fluorescent signal under a laser-scanning confocal microscopy to determine where the gene is expressed. Microscopically, epidermal cells from tobacco leaves expressing the different CsCBF fusion proteins all only showed a fluorescence signal in the nucleus. As a control, the 35S:GFP fluorescence was observed throughout the whole cell (Fig. 6B). These data indicated that CsCBF proteins were all nuclear localization proteins, which were consistent with the previous prediction. The three CsCBF transcription factors may play roles in transcriptional regulation.

Overexpression of CsCBFs enhanced the tolerance of transgenic cucumber seedlings to cold stress
The expression levels of CsCBFs were significantly induced by low temperature (Fig. 5). In order to further investigate the response of CsCBF1, CsCBF2 and CsCBF3 to cold stress and their biological functions, the agrobacterium-mediated transient transformation experiments were conducted in cucumber cotyledons to clarify CsCBFs tolerance to cold stress. Cucumber seedlings overexpressing CsCBFs with GFP fluorescence signal in cotyledons were selected for subsequent experiments (Fig. S3A). The qRT-PCR analysis showed that the expression of CsCBF genes in their transgenic cucumbers were significantly higher than those in WT  (Fig. S3B), and different stress treatments all promoted the expression of CsCBF genes in transgenic cucumber (Fig. S3C-E). To test whether CsCBFs can enhance cold resistance of cucumber cotyledons, transgenic cucumber seedlings and WT were treated at 0°C, respectively. Before treatment, the 35S:CsCBFs transgenic plants and WT all grew well, while after 3 h of cold treatment, slight wilting appeared in cotyledons of WT compared with transgenic seedling overexpressing CsCBF1, CsCBF2 and CsCBF3, and serious wilting in WT showed more obvious difference from all transgenic seedlings after 24 h (Fig. 7A). After 48 h, the survival rate of WT was only 27%, while the survival rates of the 35S:CsCBF1, 35S:CsCBF2 and 35S:CsCBF3 transgenic plants were 60%, 67% and 56% ( Fig. 7A-B).
The contents of proline and MDA are important physiological indexes to measure cold resistance of plant [39]. Compared with WT plants, the 35S:CsCBFs transgenic plants all showed a significant increase in proline and great decrease in MDA content (Fig. 7C-D). In addition, to verify that overexpression of CsCBFs enhanced cold tolerance of cucumber seedlings by regulating COR genes, the expression levels of two CsCOR genes (CsCOR15A and CsKIN1) in CsCBFs overexpressed plants and WT were detected under cold treatment, respectively. As shown in Fig. 7E-F, the transcriptional levels of the two genes in transgenic plants were significantly higher than those in WT plants after cold treatment. Based on survival rate, physiological indexes and CsCOR genes expression, CsCBF1 and CsCBF2 played stronger roles in cold stress response than CsCBF3 (Fig. 7). Taken all together, these results suggested that overexpression of CsCBF genes in cucumber could significantly enhance cold tolerance of cucumber.

CsCBFs directly activate CsCORs expression by binding to their promoters
Overexpression of CsCBFs greatly induced the expression of CsCORs (Fig. 7E-F), and previous reports have shown that CBFs can directly bind to the promoter regions of COR genes to activate their expression [12]. To verify whether CsCBFs can bind to the promoters of CsCORs, the 2-kb promoter fragments of CsCOR15A and CsKIN1 were selected and inserted into the pHIS2 plasmid, respectively. The CDSs of CsCBFs were separately cloned into the pGADT7 vector. The yeast one-hybrid (Y1H) assays were carried out and the results showed that CsCBF1, CsCBF2, and CsCBF3 proteins all could specifically directly bind to the promoters of CsCOR15A and CsKIN1, but not empty pHIS2 vector (Fig. 8A). The transient GUS activity assays were carried out in tobacco leaves to verify the above results. The above DNA fragments were separately inserted into pCAMBIA1300-GUS vector containing GUS reporter gene, and the CDSs of CsCBFs were cloned into the pCAMBIA1300 plasmid to obtain 35S:CsCBFs recombinant plasmids. The results showed that all three CsCBFs could activate the expression of CsCOR15A and CsKIN1 in vivo (Fig. 8B). These data revealed that CsCBFs could directly activate CsCORs expression by binding to their promoters.

Discussion
Cucumber (Cucumis sativus L.) is an economically important crop cultivated worldwide and one of the main vegetables grown in the facility [36]. Low temperature is a major environmental factor affecting the growth and seasonal distribution of cucumber. The transcription factor CBFs (C-repeat Binding Factor) are key "molecular switch" for plant to sense low temperature signals and regulate their adaptive responses. These genes can activate the expression of several downstream cold-tolerance related functional genes to improve plant resistance to low temperature [40]. Therefore, bioinformatics analysis of CBF gene family in cucumber was conducted to provide theoretical basis for analyzing the mechanism of CBFs regulating cucumber cold stress response.
Studying the differences in gene structure and conserved motifs is an important reference for analyzing the evolutionary relationships of gene families [33]. In our study, multi-sequence alignment analysis of CsCBFs showed that the three CsCBF proteins all contained conserved AP2 domains and their flanks (Fig. 1). By analyzing the gene structures of CsCBFs, it was found that the three CsCBF genes in cucumber did not contain UTR (Fig. 3). Phylogenetic tree analysis showed that 70 CBF proteins from cucumber and other species were divided into 18 subgroups. Among them, the three CsCBF proteins are closely related to CBF proteins in apple and pomegranate (Fig. 2). Previous studies have shown that CBF genes play important roles in plant growth and development. Overexpression of AtCBFs results in plant growth retardation and delayed flowering [16,41].
Here, we showed that light response signal elements were found on the promoters of the three CsCBF genes (Fig. 4), suggesting that CsCBFs may be involved in the regulation of light signal in the growth and development of cucumber.
Studies have shown that light is also an important environment, and photochrome-interacting factors (PIFs) play a key role in regulating plant development [42]. Part PIFs can bind to G-box and E-box cis-elements in AtCBF promoters to regulate their transcription, and the PIF3/4/7 are negatively involved in the low temperature response pathways of plants and can negatively regulate CBF expression [43]. In our study, cis-elements that respond to abiotic stresses and hormones such as low temperature, defense and stress, abscisic acid, salicylic acid, jasmonate acid, auxin and gibberellin were simultaneously screened on promoters of CsCBF genes (Fig. 4). The expression analysis showed that CsCBFs were indeed regulated by low temperature (Fig. 5A), which will help further elucidate the molecular regulation mechanism of cold signaling.
The functions of CBF genes in apple, soybean, rice, tomato, wheat, barley and maize have been widely reported [19,[44][45][46]. Overexpression of DREB can enhance plant resistance to stress, and 35S:PpDBF1 transgenic tobacco has higher salt tolerance, drought tolerance and cold tolerance [47]. The heterologous expression of ZmCBF3 and CsCBF3 (tea) in Arabidopsis can significantly enhance the frost resistance of Arabidopsis thaliana [48,49]. Overexpression of CsCBFs in the cotyledons of cucumber seedlings could markedly enhance its cold resistance and the changes of Pro and MDA in 35S:CsCBFs transgenic cucumber seedlings indicated that CsCBFs could enhance cold tolerance of cucumber (Fig. 7). Low temperature induces the expression of CBFs, thereby activating the expression of downstream target genes. Heterologous expression of VvCBF1 in Arabidopsis can enhance the expression of AtCOR15A, AtRD29B and AtRD29A [50]. The expression of AtCOR15A, AtCOR47, AtKIN1 and AtRD29A in 35S:SmCBFs transgenic Arabidopsis thaliana were up-regulated, thus enhancing its cold resistance [51]. In this study, the transcript level of CsCOR genes was significantly upregulated in 35S:CsCBFs transgenic plants after cold stress treatment (Fig. 7). Y1H and GUS experiments demonstrated that CsCBFs, as transcription factors, could directly bind to promoters of COR genes and activate their expression (Fig. 8). These results provided a foundation for further research on the function of CsCBFs gene in cold stress resistance and elucidating its mechanism.

Conclusions
In this study, we comprehensively analyzed the cucumber CBF family genes. The expression patterns of CsCBF genes under different stress treatments were also investigated, and the roles of CsCBFs in cucumber cold tolerance were analyzed in detail by transient transgenic method. This study provided a foundation for further research on the function of CsCBF genes in cold stress resistance and elucidating its mechanism.

Genome-wide identification of CsCBFs in cucumber
To identify the CsCBF genes from cucumber genome database (http:// cucur bitge nomics. org/ organ ism/ 20), six Arabidopsis AtCBF proteins were used as query sequences and Blastp was used to search for the predicted proteins. All candidate genes were further confirmed by the existence of conserved characteristic sequences (PKK/RPAGRxKFxETRHP and DSAWR) and AP2 (PF00847.20) domains using the Pfam (available online: http:// pfam. janel ia. org) and Simple Modular Architecture Research Tool (SMART) datebase (http:// smart. embl-heide lberg. de).

Physicochemical properties of CsCBF proteins
ExPASy software (http:// web. expasy. org/ protp aram/) was used to analyze protein sequences of CsCBFs to predict amino acid length and isoelectric point (pI). The location of CsCBF genes on cucumber chromosome was determined according to the physical location information in cucumber genome database.

Phylogenetic analysis, gene structure and conserved motif analysis
Multiple sequence alignments of these proteins were performed using ClustalW with default parameters. An un-rooted phylogenetic tree was constructed with the full-length amino acid sequences of the 70 CBFs using MEGA 7.0, and the neighbour-joining (NJ) method was used with the following parameters: Poisson correction, pairwise deletion, and bootstrap (1000 replicates; random seed) [52]. The corresponding DNA and cDNA sequences of each predicted gene were downloaded from genomes, and the gene structures were analyzed as described by [53]. The conserved motifs in CsCBFs were identified using Multiple Expectation Maximization for Motif Elicitation (MEME) online program (http:// memesuite. org/ index. html) [53].

Identification of cis-elements on CsCBFs promoter in cucumber
The entire cucumber genome data were downloaded from the cucumber genome database (Chinese Long 9930: http:// cucur bitge nomics. org/), and the 2 kb sequences upstream of the transcription start site of of CsCBFs were extracted by TBtools. The cis-elements on the promoter regions of CsCBF genes were analysed by PlantCARE website (http:// bioin forma tics. psb. ugent. be/ webto ols/ plant care/ html/) [52].

Plasmid construction and transient transformation of cucumber cotyledons
To generate 35S:CsCBFs-GFP, the full-length coding sequence of CsCBF1, CsCBF2 and CsCBF3 were amplified and cloned into pCAMBIA1300 vector with a GFP tag, respectively. The recombinant plasmids were transformed into Agrobacterium tumefaciens LBA4404, and then transferred into 8-d-old cucumber cotyledons for subsequent cold tolerance tests [54]. The primers used are listed in Table S3.

Subcellular localization of CsCBFs
To determine the subcellular localization of CsCBFs, the empty GFP vector and the recombinant plasmid of 35S:CsCBF1-GFP, 35S:CsCBF2-GFP and 35S:CsCBF3-GFP were injected into tobacco leaf epidermal cells, respectively. The injected tobaccos were grown under normal conditions for about 48 hours. The fluorescent signal was observed by a fluorescence microscope.

Expression pattern of CsCBF genes under different abiotic stresses
The cucumber inbred line Xintaimici was used for transient genetic transformation and stress treatments, and all plants were cultured in a light incubator under 28 °C with 16 h light /20 °C with 8 h dark cycle conditions. The two-week-old cucumber seedlings were placed in a 4 °C incubator for low temperature treatment until the leaves were collected, and the leaves of two-week-old cucumber seedlings with consistent growth were sprayed with 100 μM ABA and 100 mM NaCl, respectively. The leaves were selected at 0, 0.5, 1, 3, 6, 9, 12 and 24h for subsequent quantitative analysis. The mixed leaves of five plants were one replicate, and each sample contained three biological replicates. Each treatment was repeated at least three times.

Yeast one-hybrid assays
The 2-kb promoter fragments of CsCOR15A and CsKIN1 were selected and inserted into the pHIS2 plasmid, respectively. The coding sequences (CDSs) of CsCBFs were separately cloned into the pGADT7 vector to obtain the constructs AD-CsCBFs. The optimal 3-AT concentration which could inhibit the growth of background histidine of pHIS2 vector was screened. Then the recombinant pHIS2 vector and AD-CBFs were co-transferred to yeast strain Y187 and grown on medium SD/-Trp-Leu-His with optimal 3-AT concentration. Empty vector pGADT7 was used as the control. Y1H was conducted according to [55]. The primers used are listed in Table S3.

Transient GUS activity assays
The ORFs of CsCBFs were separately inserted into the pCAMBIA1300 vector. The 2kb promoter fragments of CsCOR15A and CsKIN1 were separately inserted into pCAMBIA1300-GUS plasmids to obtain proCsCOR15A:GUS and proCsKIN1:GUS recombinant plasmids. The tobacco leaves were used to conduct GUS activity assays. The different combinations were injected into 5-week-old tobacco leaves by agrobacterium-mediated method. Under normal conditions, the injected tobacco grew 2-3 days for subsequent experiments. The transient activity assays were measured as described previously [54].