Sobhi F.Lmlom,Yong Zhng,Bohong Su,Hito Wu,Xi Zhng,Jinong Fu,Bo Zhng,Li-Jun Qiu,*

aNational Key Facility for Gene Resources and Genetic Improvement/Key Laboratory of Crop Germplasm Utilization,Ministry of Agriculture,Institute of Crop Sciences,Chinese Academy of Agricultural Science,Beijing 100081,China

bPlant Production Department,Faculty of Agriculture Saba Basha,Alexandria University,Alexandria 21531,Egypt

cKeshan Branch of Heilongjiang Academy of Agricultural Sciences,Keshan 161606,Heilongjiang,China

dCollege of Agronomy,Northeast Agricultural University,Harbin 150030,Heilongjiang,China

eSchool of Plant and Environmental Sciences,Virginia Polytechnic Institute and State University,Blacksburg,VA 24060,USA

ABSTRACT

1.Introduction

The yield of soybean[Glycine max(L.)Merr.]is of profound importance in the global market place for its heavy utilization as a food and feed source for humans and animals[1].Branch number(BN),a major component of plant overall structure,is one of the most important factors underlying soybean yield.The number and distribution of branches on the soybean plant influence both lodging resistance and seed yield by determining the plant architecture associated with light interception[2].BN is affected by numerous environmental factors,including plant density[3],planting pattern[4],planting date[5],photoperiod[6,7],and nutrient conditions[8].BN determines the optimum number of plants per square meter and photosynthetic efficiency,both of which directly influence yield[9].Identifying genes and QTL associated with regulation of soybean BN will advance plant architecture research and cultivar development.

The tremendous advancements in soybean genomic research have accelerated the precision of QTL mapping for key agronomic traits.The availability of the soybean reference genome sequence of cultivar Williams 82 has allowed more precise QTL mapping and gene mining[10–12].In soybean,several genetic linkage maps have been constructed,employing RFLP(restriction fragment length polymorphism),SSR(simple sequence repeat),and single-nucleotide polymorphism(SNP)markers[13–16].To date,a total of 21 QTL associated with BN in soybean are listed in Soybase(https://www.soybase.org/search/qtllist by symbol)[2,17–19].Most of those QTL were distributed on 10 chromosomes:4,5,6,10,11,14,15,17,18,and 19.Among these chromosomes,6 and 11 contain major BN QTL with PVE＞10%[2,17,18].These QTL/genes can be considered possible targets for marker-assisted selection for BN.

In the present study,high branching soybean genotype Kennong 24 and low branching genotype Kenfeng 19 were identified during the breeding process and used to develop segregating populations of F2,F2:7,and F2:8.We identified a novel QTL controlling soybean BN,qBN-1,using the F2population,which was subsequently fine-mapped using RIL populations.We then developed and analyzed two backcross populations,and narrowed down the region of qBN-1 to two putative candidate genes.Based on the qRT-PCR analysis,Glyma.06G208900 was selected as a novel candidate gene controlling BN.Results from this study can be valuable in future cloning and transformation efforts for the BN gene as well as genetic improvement of plant architecture of soybean.

2.Materials and methods

2.1.Population development and phenotyping of BN

The RIL population included 599 lines derived from a cross between Kennong 24(high branching number)and Kenfeng 19(low branching number).F1plants were self-pollinated in 2015 to produce an F2population consisting of 599 individuals in 2015,from which the individuals were.The population was advanced by single-seed descent to produce F2.7and F2:8lines.The parental genotypes along with the F2,F2:7,and F2:8populations were planted along with the parental genotypes in Keshan of Heilongjiang province in 2015,2018,and 2019,respectively.The row length,row spacing,and plant interval were 3.00,0.40,and 0.10 m,respectively.At maturity,edge plants were removed and data for BN per plant was from individuals in F2and for 10 randomly selected plants from each line in F2:7and F2:8.BN of each plant was manually recorded as the number of effective branches on the main stem that had two or more nodes with at least one mature seed pod at harvest[2,18].

To confirm the effect of branch-number QTL detected in the current study,two backcross populations were constructed.KN24 was used as a recurrent parent to cross with F1plants derived from KN24 and KF19 to produce KN24BC1F1plants,which were selfed to develop a 1305 BC2F2population.At the same time,KF19 was backcrossed to selected F1plants to obtain BC3F1plants.These were selfed to develop 1712 BC3F2plants.The BC2F2and BC3F2populations were planted in Keshan county in 2018.The cultivation conditions were the same as those used for the RILs in 2018.A two-tailed Student’s t-test was performed to determine the statistical significance between KF19BC3F2and KN24BC2F2.A total of 42 SSRs(http://soybase.org)were used to evaluate the relationship between two sets of selected lines and their parents.Compared to KN24,both KN24BC2F2-4-1 and KN24BC2F2-4-14 showed similarities of 95.42% and 94.10%,respectively,and compared to KF19,KF19BC3F2-14-88 and KF19BC3F2-4-92 showed similarities of 94.74%,and 96.04%,respectively.The four lines and their parents were used for RNA expression analyses.

2.2.Bulked segregant analysis(BSA)

To study the molecular genetic mechanism underlying BN in soybean,we performed QTL analysis of BN in F2population.Bulked segregant analysis(BSA)was performed to identify SSR markers linked to genes influencing BN.Two bulks with contrasting BN were constructed.One bulk was designed by combining the genomic DNA of 50 plants showing the extreme BN of eight.Likewise,the other bulk combined the genomic DNA of plants with BN zero.The two bulks,as well as their parents,were screened with 543 SSR markers covering all 20 chromosomes to identify polymorphic markers that are possibly associated with the BN.

2.3.Genotyping with SSR markers

Fresh leaves were used to extract DNA according to the cetyltrimethylammonium bromide(CTAB)[20]method with some modification.DNA sequences of 14 SSR markers were obtained from Soybase(http://soybase.org)(Table 1).The PCR program was as follows:94 °C for 5 min,34 cycles of 30 s at 94°C,30 s at 55°C,and 1 min at 72°C,with a final extension of 5 min at 72°C.PCR products were then separated on 6%(w/v)polyacrylamide gel followed by silver staining.

2.4.QTL mapping and verification

QTL were identified by inclusive composite interval mapping(ICIM)using IciMapping software version 4.1 with mapping function Kosambi[21–23].The logarithm of odds(LOD)score threshold was set to 3 with a recombination frequency of 0.3.The adjusted window size was 5 cM.The ICIM-ADD model for biparental populations(BIP)was adopted.The LOD threshold was computed by 1000 permutations at the 0.05 probability level,and the walking speed along the chromosome was set to 1.0 cM[21,23].

Table 1–Fourteen polymorphic markers on chromosome 6 used to localize branch-number QTL(qBN-1).

2.5.RNA extraction,reverse transcription PCR,and quantitative real-time PCR

For RNA extraction,fresh tissues of leaves,stem,axillary meristem,shoot apical meristem(SAM),and root were collected from KN24,KF19,and the two backcross populations(BC2F2with homozygous with KN24 and BC3F2homozygous with KF19)at the vegetative stage V4,when the fourth trifoliate leaves appear,and at reproductive stage R1,when one flower blooms at any node on the main stem.Samples were frozen in liquid nitrogen and kept at-80 °C.Using an RNA Prep Pure Plant kit(Tiangen Co.,Beijing,China),total RNA from each tissue was extracted and treated with DNaseI(Thermo Fisher Scientific Inc.,Grand Island,NY,USA).A cDNA synthesis was performed using a SuperScript II kit(TaKaRa Biotechnology,Dalian,Liaoning,China).Real-time PCR was performed using an SYBR Premix Ex Taqkit(TaKaRa)on an ABI 7300 Real-Time PCR System.Three replicates were run for each sample.The soybean Actin11 gene was used as an internal control[24].The relative expression level against the Actin11 gene was quantified using the 2-ΔΔCTmethod[25].

3.Results

3.1.Phenotypic variation in the populations and their parents

There was a significant difference in BN between the two parents.In 2015,the BN mean was 6.0 for KN24 and 0.6 for KF19,respectively.While,in 2018,the mean BNs were 7.9 for KN24 and 1.1 for KF19,respectively.A slight reduction in BN was observed in 2019,with an average BN 6.2 for KN24 and 0.8 for KF19,respectively(Table 2).The phenotypic stability of BN for KN24 and KF19 over three years indicated that they were suitable for use as parents to develop segregating populations and identify QTL for BN.Tables 2 and 3 show the descriptive statistics of the BN of all populations including F2in 2015,F2:7,BC2F2,and BC3F2in 2018 and F2:8in 2019.The BN that was investigated five populations showed continuous variation and followed the normal distribution(Fig.1).

3.2.Preliminary mapping of QTL for BN in an F2 population

Of 543 SSR markers screened for polymorphism,only 250 distributed across 20 chromosomes were polymorphic between KN24,KF19 and two DNA pools derived from the F2population.Of these,two polymorphic markers on chromosomes 6 and 18,corresponding to the genomic regions of qBN-1 and qBN-2,respectively,appeared to govern BN phenotype.We found out two SSR markers of BARCSOYSSR_06_0717 and BARCSOYSSR_06_1441 that associated with qBN-1.To map the region of qBN-1,we obtained 10 polymorphic SSR markers within this region from Soybase(https://www.soybase)and identified 599 individuals of F2population.Eventually,the qBN-1 locus was localized to a 2.2-Mb region on chromosome 6 flanked by BARCSOYSSR_06_0993 and BARCSOYSSR_06_1070(Fig.2).A total of 84 gene models were present in this region(Glyma.W82.a2.V1.1).

Table 2–Descriptive statistics of BN in parental lines of KN24 and KF19 and their segregating populations of F2,F2:7,and F2:8 in 2015,2018,and 2019.

3.3.Fine mapping and validation of qBN-1

To fine-map qBN-1,an F2:7RIL population of 599 lines was developed and screened with 10 SSR markers that showed polymorphism between two parents.Then the qBN-1 was located within a 115.67-kb genomic region on chromosome 6 between markers BARCSOYSSR_06_1048 and BARCSOYSSR_06_1053.Four additional SSR markers obtained from Soybase(https://www.soybase.org)were used to confirm the qBN-1 in F2:8lines.Results showed that qBN-1 was localized to the same region,with the highest R2value(22.69%)and the highest LOD score(32.07)as in the F2and F2:7populations(Table 4).Two SSR markers,BARCSOYSSR_06_0993 and BARCSOYSSR_06_1035,were also used for genotyping 1712 BC3F2and 1305 BC2F2plants.Then,870 plants recombinant between these two markers were genotyped with the 10 SSR markers.Further screening of 370 recombinants in the region between markers BARCSOYSSR_06_0993 and BARCSOYSSR_06_1070 was performed for fine-mapping qBN-1.This approach further confirmed the localization of qBN-1 between markers BARCSOYSSR_06_1048 and BARCSOYSSR_06_1053(Fig.3).This region harbored only two genes.One,Glyma.06G208800,encodes a calcium-binding cml15-related protein and the other,Glyma.06G208900,a phospholipidtranslocating ATPase.

Table 3–Descriptive statistics of BN in two backcross populations.

3.4.Expression patterns of two candidate genes

To confirm which of genes related to BN,qRT-PCRs were performed to analyze the expression pattern of the two genes in KN24 and KF19 as well as two sets of selected lines.qRTPCR showed that the two genes Glyma.06g208800 and Glyma.06g08900 showed different expression levels for lowand high-branching genotypes in multiple tissues and at the V4 and R1 stages(Fig.4).Glyma.06G208800 gene showed low expression in all tissues,but only leaves at the R1 stage showed a significant difference between KN24 and KF19 as well as between two sets of selected lines:KN24BC2F2-4-1/KN24BC2F2-4-14 and KF19BC3F2-14-88/KF19BC3F2-4-92.However,Glyma.06G208900 showed much higher expression than Glyma.06G208800.The expression of Glyma.06G208900 appeared to be low and showed no significant difference in root and stem at the V4 and R1 stages,whereas the expression of Glyma.06G208900 in leaves was higher but significantly different only at the R1 stage.Compared to tissues of roots,stems,and leaves,the expression of Glyma.06G208900 was higher level differed significantly between high-and lowbranching genotypes in both axillary and shoot apical meristems.However,the expression of Glyma.06G208900 was significantly lower in high-branching than in low branching genotypes in the axillary meristem and was significantly higher in high-branching than in low-branching genotypes in the shoot apical meristem.

4.Discussion

The BN of soybean plants is a major determinant of soybean architecture and seed yield.Given that it is controlled by complex spatial-temporal regulation of axillary bud outgrowth following axillary meristem initiation[26,27],detecting the genetic basis holds great potential to enhance the breeding selection and efficiency in order to develop high yielding soybean cultivars.In the present study,qBN-1 was identified in a large interval of about 24 Mb on chromosome 6 between the SSR markers BARCSOYSSR_06_0717 and BARCSOYSSR_06_1441 using BSA.In the F2generation,we used 10 polymorphic markers within the big region to narrow down qBN-1 to a physical interval of 2.2 Mb between markers BARCSOYSSR_06_0993 and BARCSOYSSR_06_1070.The qBN-1 region overlapped with those of QTL previously reported by Chen et al.[18]and Shim et al.[17].

To delimit the candidate region,F2:7and F2:8RILs were screened using an additional SSR markers.The qBN-1 was fine-mapped between SSR markers BARCSOYSSR_06_1048 and BARCSOYSSR_06_1053 with a candidate region of 115.67 kb.Using the backcross population,qBN-1 was confirmed to lie in the same candidate region between markers BARCSOYSSR_06_1048 and BARCSOYSSR_06_1053.Although a large number of BN-related QTL have been detected in the soybean genome,all have shown a considerably large confidence interval[2,17–19,28].For instance,three branching-number QTL were mapped on chromosomes 3,6,and 10 using an F2population with 154 plants derived from a cross between Charleston and Dongnong 594[18].A BN QTL was located on chromosome 10 using 126 F5recombinant inbred lines(RILs)[28].Five QTL,qBr1 to qBr5 on chromosomes 6,19,11,17,and 18 respectively,associated with BN have been identified using F9recombinant inbred lines(RILs)[2].Among these,a major QTL on chromosome 6(qBR6-1),with a LOD score of 10.3 and 14.5% of the phenotypic variation in BN,contained 13 genes[17].However,based on the genome sequence of Glycine max Wm82.a2.v1 in Phytozome 12(https://phytozome.jgi.doe.gov/pz/portal.html),in our finemapping region,there were only two predicted genes within this 115.7 kb interval(between BARCSOYSSR_06_1048 and BARCSOYSSR_06_1053)which may control BN.We further discovered that qBN-1 was located near the E1 gene.Sayama et al.[2]reported that the E1 gene had a putative pleiotropic effect on BN.Also,another study revealed that E1 gene significantly correlated with the BN in various F2populations and one F2:3population[29].It was reported that the E1 locus was linked to the Satt365 marker in linkage group(LG)C2(Gm06)[30],but in our study,the E1 locus located close to marker 06-1043 was about 191 kb from qBN-1.Resequencing analysis of two parents showed no SNP variation or InDel in the E1 locus.Therefore,E1 might be not correlated to BN in our study.

Fig.1–Distribution of BN in five populations.(a)599 F2 plants in 2015.(b)599 F2:7 lines in 2018.(c)599 F2:8 lines in 2019.(d)1305 KN24BC2F2 plants in 2018.(e)1712 KF19BC3F2 plants in 2018.

The outgrowth of axillary buds is inhibited by the active shoot apex,in a phenomenon referred to as apical dominance[31].Decapitation abolishes apical dominance and triggers the growth of one or more axillary meristems because auxin,which is synthesized in the shoot apex,is mobilized to the lower parts of plants and inhibits branch outgrowth[31].Considering the relevance of shoot apex to branch development,we compared the expression levels of two genes showing transcriptional activity in the SAM and axillary meristems.Glyma.06G208900 was significantly downregulated in axillary meristems in high branch-number genotypes and upregulated in SAM,strongly suggesting that this gene influences BN in soybean.

Bioinformatics and sequence homology analysis(https://soybase.org/)revealed that Glyma.06G208800 encodes a CML15-related calcium-binding protein.Numerous studies indicated that CMLs functions are associated with both biotic and abiotic stresses development,and considerable evidence has demonstrated that these proteins are not probable to have redundant roles or,in contrast,play vital and precise roles in coordinating ecological responses of plants.In addition,many CMLs are now known to recognize a specific target[32–38].CMLs also contribute to several aspects of plant development including CML42 involved in trichome branching[39],CML25 and CML7 in root hair elongation[40,41],CML39 in early seedling establishment[42],and CML23 and CML24 in flowering[43].We also identified a

phospholipid-translocating ATPase(Glyma.06G208900)which was upregulated in the SAM and downregulated in axillary meristems in KN24 and high-branching lines.The P-type-ATPase influences cell growth and auxin trigger the proton pump,resulting in loosening of the cell wall,which occurs when the acid-labile bonds are broken or via the initiation of lytic enzymes within the membrane activating lytic enzymes within the wall[44].Recent reports confirm that auxininduced cell elongation involves auxin-mediated regulation of H+-ATPase activity by phosphorylation[44,45].Regulation of the PM H+-ATPase is achieved by important factors that control plant physiology such as hormones,environmental stresses,phytohormones,phytotoxins,and light[46].In particular,several studies showed that auxin specifically increased the level of PM H+-ATPase by two to three times in elongating tissues such as maize(Zea mays)coleoptiles[47,48].Therefore,regulation of H+-ATPase activity might play a dual role during phototropism to modulate the proportion of protonated auxin and thereby auxin influx,and to promote cell wall acidification and thereby cell elongation.PM H+-ATPase activity was measured in PM-enriched fractions isolated from peach tree buds or their underlying tissues[49].In Arabidopsis thaliana,P4-ATPases are vital for establishing phospholipid asymmetry between the two leaflets in the lipid bilayer.In eukaryotic cells,this is an essential transport activity required for generating an initial membrane curvature preceding vesicle budding in both endocytosis and exocytosis[50–52].The influence of phospholipidtranslocating ATPase on BN for Glyma.06G208900 awaits further study.Our results may support the cloning of a gene influencing BN and the development of functional markers for marker-assisted selection in soybean breeding.Moreover,our findings provide valuable information for identifying the candidate gene of the qBN-1 for BN in the near future.

Fig.2–Identification of qBN-1 by QTL mapping.Marker names are shown to the left of plots.(a)QTL plot of the qBN-1 locus in the F2 generation between markers BARCSOYSSR_06_993 and BARCSOYSSR_06_1070,identified in 599 F2 plants.(b)qBN-1 between BARCSOYSSR_6_1048 and BARCSOYSSR_6_1053,identified in a 599-line F2:7 RIL population.(c)qBN-1 between BARCSOYSSR_6_1048 and BARCSOYSSR_6_1053 confirmed in a 599-line F2:8 RIL population.

Fig.3–Validation of qBN-1 locus.(a)1712 KFBC3F2 and 1305 KN24BC2F2 plants were screened with two SSR markers:BARCSOYSSR_06_0993 and BARCSOYSSR_06_1135,to identify recombinants.(b)The qBN-1 locus was localized between markers BARCSOYSSR_06_993 and BARCSOYSSR_06_1070 on chromosome 6.(c)The qBN-1 locus was further localized to a 115.67-kb region flanked by markers BARCSOYSSR_06_1048 and BARCSOYSSR_06_1053.Numbers below SSR marker names are numbers of recombinant plants.Letters to the right of“Phenotype”values indicate significant differences between recombinant lines.

5.Conclusions

In this study,qBN-1,a major QTL controlling BN on soybean chromosome 6,was fine mapped to 115.67-kb in RIL population of F2:7and F2:8and confirmed in backcrossing population of KF19BC3F2and KN24BC2F2.Within this interval,two putative candidate genes were analyzed.Findings of expression analysis showed that Glyma06G208900 may be the candidate gene within the region of qBN-1.The findings of this study elucidate the feasibility to implement qBN-1 into other soybean breeding lines using marker-assisted selection.

Declaration of competing interest

Authors declare that there are no conflicts of interest.

Acknowledgments

This research was supported by the National Key Research and Development Program of China(2016YFD0100201 and 2016YFD0100304),the Platform of National Crop Germplasm Resources of China(2016-004 and 2017-004)and the Agricultural Science and Technology Innovation Program(ASTIP)of the Chinese Academy of Agricultural Sciences.

Fig.4–Relative expression levels of Glyma.06G208800 and Glyma.06G208900 in KN24,KF19 and two sets of near-isogenic lines(NILs)at two stages:V4 and R1.**indicated expression level was significant different at P＜0.05.

Author contributions

Li-Juan Qiu conceived and supervised the research project.Sobhi F.Lamlom performed data analysis and QTL mapping and wrote the manuscript.Yong Zhang and Haitao Wu collected phenotyping data.Bohong Su,Xia Zhang,and Haitao Wu provided assistance and advice on experiment.Li-Juan Qiu,Jindong Fu,and Bo Zhang revised and improved the manuscript.