Zhaobo Li ,Qun Chng ,Zhuoran Gan ,Zhihong Hou ,Yuhang Zhang ,Yongli Li ,Haiyang Li ,Haiyang Nan ,Cn Yang ,Linnan Chn ,Sijia Lu ,Wnqian Shi ,Liyu Chn ,Yanping Wang ,Chao Fang ,Liping Kong ,Tong Su ,Shihn Li ,Kun Kou ,Lingshuang Wang ,Fanjiang Kong a,,*,Baohui Liu a,,*,Lidong Dong ,*

a The Innovative Academy of Seed Design,Key Laboratory of Soybean Molecular Design Breeding,Northeast Institute of Geography and Agroecology,Chinese Academy of Sciences,Harbin 130102,Heilongjiang,China

b University of Chinese Academy of Sciences,Beijing 100081,China

c Innovation Center of Molecular Genetics and Evolution,School of Life Sciences,Guangzhou University,Guangzhou 510006,Guangdong,China

d Agronomy College of Heilongjiang Bayi Agricultural University,Daqing 163319,Heilongjiang,China

e National Center for Soybean Improvement,National Key Laboratory of Crop Genetics and Germplasm Enhancement,Jiangsu Collaborative Innovation Center for Modern Crop Production,Nanjing Agricultural University,Nanjing 210095,Jiangsu,China

f Mudanjiang Branch,Heilongjiang Academy of Agricultural Sciences,Mudanjiang 157000,Heilongjiang,China

ABSTRACT Flowering time is an important agronomic trait for soybean yield and adaptation.However,the genetic basis of soybean adaptation to diverse latitudes is still not clear.Four NIGHT LIGHT-INDUCIBLE AND CLOCK-REGULATED 2 (LNK2) homeologs of Arabidopsis thaliana LNK2 were identified in soybean.Three single-guide RNAs were designed for editing the four LNK2 genes.A transgene-free homozygous quadruple mutant of the LNK2 genes was developed using the CRISPR (clustered regularly interspaced short palindromic repeats)/Cas9(CRISPR-associated protein 9).Under long-day(LD)conditions,the quadruple mutant flowered significantly earlier than the wild-type (WT).Quantitative real-time PCR (qRT-PCR)revealed that transcript levels of LNK2 were significantly lower in the quadruple mutant than in the WT under LD conditions. LNK2 promoted the expression of the legume-specific E1 gene and repressed the expression of FT2a.Genetic markers were developed to identify LNK2 mutants for soybean breeding.These results indicate that CRISPR/Cas9-mediated targeted mutagenesis of four LNK2 genes shortens flowering time in soybean.Our findings identify novel components in flowering-time control in soybean and may be beneficial for further soybean breeding in high-latitude environments.

Keywords:Soybean LNK2 CRISPR/Cas9 Genome editing Flowering time

1.Introduction

Soybean (Glycinemax (L.) Merr.) is one of the most important legume crops and provides the major plant protein for human food,animal feed,and industrial uses [1].It is also a typical short-day(SD)plant and highly sensitive to photoperiod[2,3].Soybean flowering is a key production trait that is sensitive to seasonal changes in day length.Twelve maturity loci (E1toE11andJ) controlling flowering and maturity time have been identified and characterized at the phenotypic and genetic levels in soybean [4–13].Among these loci,E1has a large influence on flowering time,and the E1 protein contains a B3-like domain [14].E2is an ortholog of theArabidopsis GIGANTEA(GI) gene [15],andE3andE4are homologs ofArabidopsis PHYTOCHROME A(PHYA2) andPHYA3genes [16,17].TwoFLOWERING LOCUS T(FT) homologs,FT2aandFT5a,have been identified as encoding components of florigen and controlling flowering time in soybean [18,19].The expression levels of these two genes are down-regulated byE1,E2,E3,andE4under long-day (LD) conditions,suggesting thatFT2aandFT5aare the major soybean flowering regulation targets [15,18].Subsequent studies have shown thatE9isFT2a,with a Ty1/copia-like retrotransposon SORE-1 inserted in its first intron in its recessive allele [20].Recessivee9allele causes delayed flowering [20].Several flowering-related genes have been cloned by reverse genetics in soybean.For instance,RAV belongs to the AP2/ERBP transcription factor family,and it is involved in plant development and flowering time in soybean [21].col1bmutants flower earlier than wild-type (WT) plants under LD conditions,and overexpression ofCOL1ain soybean plants causes later flowering than in WT plants under LD or natural conditions [22].

Recently,two novel circadian components were identified inArabidopsis thaliana:NIGHT LIGHT-INDUCIBLE AND CLOCKREGULATED 1(LNK1) andLNK2[23].LNK1 and LNK2 function as coactivators of dawn-phased MYB-like transcription factors,such as RVE4 and RVE8,to regulate target gene expression;the target genes includeTIMING OF CAB EXPRESSION 1(TOC1) andPSEUDORESPONSE REGULATOR 5(PRR5) [24–27].Other functions of LNKs in growth control and domestication have also been reported.InArabidopsis,the LNK family redundantly controls petiole length and leaf growth,which in turn influence rosette size and biomass accumulation in adult plants [28].In tomato,lnk2mutations led to a lengthened circadian period,likely adapting the cultivated species to the long summer days it encountered as it was moved from its equatorial origin [29].However,the potential functions of the LNK family members in soybean are still unknown.

The CRISPR (clustered regularly interspaced short palindromic repeats)/Cas9 (CRISPR-associated protein 9) system is a simple,efficient,and highly specific genome editing tool for soybean[30–32].Recently,it has been successfully applied to identify gene function and develop germplasm in soybean via editing of multiple genes[19,33–36].In this study,we identified the functional homeologs ofAtLNK2in soybean and edited soybeanLNK2with the CRISPR/Cas9 system to generate a homozygous quadruple mutant of theLNK2 genes.The quadruple mutant showed accelerated flowering under long photoperiods.Our results indicate that novel genes in theLNK2family control flowering time in soybean and may help drive improvements in soybean breeding in high latitude regions.

Fig.1.Phylogenetic tree of LNK2 proteins in plants.Dark blue and light blue shading indicates LNK2 proteins in legumes.Yellow shading indicates LNK2 proteins in other plants.Red font indicates four LNK2 proteins in soybean.

Fig.2.Schematic diagram of simultaneous site-directed mutagenesis of LNK2a, LNK2b, LNK2c,and LNK2d using the CRISPR/Cas9 technology.(A) Exon–intron structure of LNK2a,LNK2b, LNK2c,and LNK2d.Black boxes indicate exons and red stripes are editing sites.Continuous lines indicate introns and untranslated regions(UTRs).Nucleotide sequences indicate regions targeted by the sgRNAs designed in this study;nucleotides in red indicate the proto-spacer adjacent motif (PAM).(B) Schematic of the T-DNA region of the targeting vector designed for mutagenesis of the LNK2 genes using the CRISPR/Cas9 system. U3b indicates the Arabidopsis U3b promoter. U3d indicates the Arabidopsis U3d promoter. U6-1 indicates the Arabidopsis U6-1 promoter.sg1,sg2,and sg3 indicate the sgRNAs. Cas9 and Bar gene were expressed by 35S promoter. Bar indicates the phosphinothricin-resistance gene.LB,left border;RB,right border.

2.Materials and methods

2.1.Plasmid construction

gRNAs of theLNK2genes were designed using the online tool CRISPR-P (http://cbi.hzau.edu.cn/crispr/).Pairs of DNA oligonucleotides of the two sgRNAs were synthesized by Tsingke Biological Technology(Tsingke,Guangzhou,China)and annealed to generate dimers,which were then integrated into the pYLCRISPR/Cas9P35SB vector following Ma et al.[37].Cas9was expressed under the 35S promoter.Positive plasmids were introduced intoAgrobacterium tumefaciensstrain EHA101 for stable transformation.

2.2.Stable soybean transformation

The soybean cultivar Williams 82(W82)was transformed using theAgrobacterium-mediated method of Zeng et al.[38].T0generation transgenic soybean plants were screened by wiping 100 mg L-1glufosinate-ammonium solution onto the upper leaf surface at three vegetative stages (V3,V4,and V5).Transgenic plants (T0,T1)were identified by PCR amplification and sequencing.Different types of gene editing could be identified by sequence peaks.The heterozygous mutations showed superimposed peaks from the target sites to the end.Wild-type and homozygous mutations showed no superimposed peaks at the target sites.They were advanced to the T2generation for further analysis.PCR analysis was performed to confirm the presence of theCas9gene,usingCas9gene-specific primers (Table S1).

2.3.Quantitative real-time PCR (qRT-PCR) analysis

To assess the expression of theLNK2,E1,andFT2agenes in the quadruple mutant and WT plants,all of the soybean plants were grown in a growth chamber (Conviron,Winnipeg,Canada) at a constant temperature of 25 °C with an average light fluence rate of 200–300 μmol m-2S-1.Daylength schedules were 16 h light/8 h dark for LD conditions.The WT plants and quadruple mutant ofLNK2soybean plants under LD conditions were sampled at 20 days after emergence (DAE) and immediately frozen in liquid nitrogen.RNA was extracted using the RNA Prep Pure Plant Kit(CWBIO,Beijing,China).cDNA was synthesized using Oligo(dT)18 primer and the First cDNA transcriptase kit (Takara,Dalian,China).qRT-PCR was performed using a LightCycler 480 SYBR Green I Master(Roche,Mannheim,Germany) in Roche LightCycler480 system(Roche,Mannheim,Germany).A soybean house-keeping gene,Tubulin(Glyma.05G157300) was used as an internal control.The relative transcript levels of the target gene were analyzed following the relative quantification method (2-ΔΔCT) [39].

2.4.Phenotypic measurement

All plants were grown in the growth chamber under either LD conditions or short-day (SD) conditions (12 h light/12 h dark) at a constant temperature of 25 °C with an average light fluence rate of 200–300 μmol m-2S-1.The flowering and pod production times of each quadruple mutant and WT plant were recorded as days from emergence to the R1 stage (when the first flower appears at any node in the main stem) [40] and the pod setting stage.Plant height and node number were recorded every 4 days from 20 to 40 DAE under LD conditions.For analyses of flowering time,the pod production time and pod number of at least 10 plants were recorded.

2.5.Statistical analyses

For phenotypic evaluation,10 individual plants per genotype were used.For expression analyses using qRT-PCR,at least three plants were pooled per tissue sample,and at least three qRT-PCR reactions (technical replicates) were performed.The mean values for each measured parameter were compared using one-tailed,two-sample Student’st-tests in Microsoft Excel (Microsoft,Redmond,USA),whichever was appropriate.

2.6.Genetic marker development

TheLNK2a,2b,2c,and2dsequences of the W82 and mutant genomes were obtained by sequencing.Primers were designed with a product size of <210 bp.Two InDels (insertion-deletions)were developed based on the variations in theLNK2aandLNK2bgenes.One derived cleaved amplified polymorphic sequence(dCAPS) marker was developed based on the variations inLNK2c.One cleaved amplified polymorphic sequence (CAPS) marker was developed on the basis of the variations in theLNK2dgene.

Fig.3.Homozygous of soybean lnk2a2b2c2d generated by the CRISPR/Cas9 system in the T1 generation.(A)PCR-based genotype of the Cas9 gene of the transgenic soybean plants.‘‘+”indicates the Cas9 gene was detected.‘‘-”indicates the Cas9 gene was not detected.M indicates a marker.T1-1-1 to T1-1-8 indicate T1 plants from T0-1.T1-2-1 to T1-2-9 indicate T1 plants from T0-2.P indicates positive control.N indicates negative control.(B) Sequences of WT and mutant plants at target sites LNK2a-sg1, LNK2b-sg1,LNK2c-sg2, LNK2d-sg3,and LNK2d-sg2.Dashes indicate deleted nucleotides.Nucleotides in red indicate PAM.(C) The sequence peaks of the wild-type (WT) and mutants at target sites LNK2a-sg1, LNK2b-sg1, LNK2c-sg2, LNK2d-sg3,and LNK2d-sg2.Red arrowheads indicate mutation locations.

3.Results

3.1.Phylogenetic analysis of the LNK2 protein in plants

Protein sequence alignment was used to identify LNK2 protein homologs in legumes and other plants.Four soybean LNK2 homologs were named LNK2a (Glyma.04G141400.1),LNK2b (Glyma.11G154700.1),LNK2c (Glyma.13G199300.1),and LNK2d(Glyma.15G237600.1).Phylogenetic analysis showed that the four soybean LNK2 proteins were closer to LNK2 proteins in some legumes than to those in other plants(Fig.1).In legumes,the four soybean LNK2 proteins were closer to the LNK2 proteins inCajanuscajan(KK1022707),Phaseolus vulgaris(Phvul.011G209100.1,Phvul.005G029800.1),Trifolium pratense(Tp57577),andGlycine soja(XP 028228790.1) than to those inVigna unguiculata(QCD85557.1),V.radiata(XP014508800.1),V.angularis(XP017433373.1),Mucuna pruriens(RDX69660.1),Lupinus angustifolius(XP019413555.1),Arachis duranensis(XP015946873.1),Trifolium subterraneum(GAU39308.1),andMedicago truncatula(ABD28297.1) (Fig.1).

3.2.CRISPR/Cas9 system for targeted mutagenesis of soybean LNK2 genes

Three synthetic single-guide RNAs were designed for simultaneous knockout ofLNK2a,LNK2b,LNK2candLNK2din soybean.Three sgRNAs:sgRNA1(sg1),sgRNA2(sg2)and sgRNA3(sg3),were used to target the second exons ofLNK2a,LNK2b,LNK2c,andLNK2d(Fig.2A).Sg1 precisely targeted theLNK2aandLNK2bgenes,sg2 and sg3 targeted theLNK2candLNK2dgene (Fig.2A).The CRISPR vector encodesCas9was driven by the CaMV35S promoter and the three gRNAs were driven by theArabidopsis U3b,U3d,andU6-1promoters (Fig.2B).

Fig.4.Analysis of soybean LNK2 genes of the quadruple mutant and WT.(A)Alignments of the LNK2 amino acid sequences of the quadruple mutant and WT.(B–E)qRT-PCR analysis of 24-h expression dynamics of LNK2 in the leaves of the quadruple mutant and WT at 20 days after emergence (DAE) under LD (16 h light/8 h dark) conditions.Values shown are relative to the control gene Tubulin (TUB) and represent means ± standard error of the mean (s.e.m.) of three biological replications.Bars indicate s.e.m.Black and white boxes represent dark and light periods,respectively.

3.3.Transgene-free homozygous quadruple mutant of soybean LNK2

Transformation of 200 cotyledonary nodes of W82 yielded 13 T0lines.Twelve transgenic lines were identified by PCR usingCas9gene-specific primers (Fig.S1A).Sequencing of theLNK2genes showed superimposed peaks inLNK2a-sg1,LNK2b-sg1,LNK2csg2,LNK2d-sg2,andLNK2d-sg3,but not inLNK2c-sg3 (Fig.S1B).Thus,the transgene-encoded Cas9 and sgRNAs induced doublestrand breaks in theLNK2genes.Plants T0-1 and T0-2 were heterozygous quadruple mutants (Table S2).

The T1plants derived from plants T0-1 and T0-2 were screened to identify soybean homozygous quadruple mutants ofLNK2without the transgene.PCR amplification ofCas9showed that two T1plants derived from plant T0-1 and five T1plants from plant T0-2 lackedCas9(Fig.3A;Table S3).Sequencing of theLNK2 genes of Cas9-negative T1plants showed that one plant (T1-2-3) was a homozygous quadruple mutant ofLNK2(Table S3).It carried a 35-bp deletion atLNK2a-sg1,a 10-bp deletion atLNK2b-sg1,a 5-bp deletion atLNK2c-sg2,and a 1-bp deletion at bothLNK2d-sg2 andLNK2d-sg3 (Fig.3B,C),resulting in frameshift mutations in the LNK2a,LNK2b,LNK2c,and LNK2d proteins (Fig.4A).

3.4.Expression levels of LNK2a,2b,2c,2d in the quadruple mutant and WT plants

The expression levels ofLNK2a,LNK2b,LNK2c,andLNK2din the quadruple mutant and WT plants showed similar patterns,reaching their maximum values at 4 h and thereafter decreasing(Fig.4B–E).However,the transcript levels ofLNK2in the quadruple mutant were significantly lower than those in WT plants(Fig.4B–E).Thus,LNK2 may indirectly influence its own transcription level.

Fig.5.The quadruple mutant showed early flowering and pod setting under LD(16 h light/8 h dark)conditions.(A)Flowering phenotypes of WT and the quadruple mutant(T1-2-3)under LD conditions.(B)Statistics of the flowering time of WT and the quadruple mutant(T1-2-3)under LD conditions.DAE,days after emergen ce.(C)Statistics of the flowering time of WT and the quadruple mutant(T1-2-3)under SD(12 h light/12 h dark)conditions.(D)Pod setting phenotypes of WT and the quadruple m utant(T1-2-3)under LD conditions.(E)Statistics of the pod setting stage of WT and the quadruple mutant(T1-2-3)under LD conditions.(F)Statistics of the pod sett ing stage of WT and the quadruple mutant (T1-2-3) under SD conditions.Bars indicate the standard error of the mean (s.e.m.).All values are presented as mean ± s.e.m.(n=10 plants).Significant differences were identified by Student’s t-test (**, P <0.01).

Fig.6.LNK2 regulates the expression of E1 and FT2a. (A) Transcriptional levels of E1 in WT and the quadruple mutant of soybean T1-2-3 under LD (16 h light/8 h dark)conditions.(B)Transcriptional levels of FT2a in WT and the quadruple mutant of soybean T1-2-3 under LD conditions.Values shown are relative to the control gene TUB and represent means ± standard error of the mean (s.e.m.) of three biological replications.Significant differences were identified by Student’s t-test (**, P <0.01).

3.5.Analysis of the phenotype of the quadruple mutant of LNK2

Loss of function ofLNK2in the quadruple mutant significantly accelerated flowering (Fig.5A,B) and shortened pod production time (Fig.5D,E) under LD conditions.Under SD conditions,the flowering and pod production times were unchanged (Fig.5C,F),and were shorter than those under LD conditions.By qRT-PCR,thelnk2quadruple mutation suppressedE1expression and thus releasedFT2aexpression to facilitate flowering (Fig.6).These results indicate thatLNK2is involved in flowering time control in soybean under LD conditions.

Plant height and node number did not differ between the quadruple mutant and WT plants under LD conditions (Fig.S2A,B).Thus,soybeanLNK2plays a key role in flowering time but has no effect on plant height.

3.6.Development of genetic markers and inheritance of quadruple mutant alleles

The two InDels,one dCAPS,and one CAPS marker developed to identify the T1-2-3 mutant alleles are shown in Fig.7A.For the genotyping of the T1-2-3 mutants,PCR amplifications were performed usingLNK2-specific and InDel,dCAPS,and CAPS-specific primer pairs.The amplified products ofLNK2aandLNK2bon the WT genomic DNA templates were larger than those on the mutant genomic DNA templates(Fig.7B).The amplified products ofLNK2con the WT genomic DNA templates were cleaved by the restriction endonucleaseMboI.The sizes of products after cutting withMboI were 18 and 162 bp,but this was not the case for the mutant genomic DNA templates.Because the 18-bp fragment was too small to detect,only the 162-bp bands for WT appear in Fig.7B.Likewise,the amplified products ofLNK2don the WT genomic DNA templates were cleaved by the restriction endonucleaseMspI.The sizes of products after cutting withMspI were 20 and 186 bp,but this was not the case for the mutant genomic DNA templates,for which only the 186-bp fragment of the WT appears in Fig.7B.

Fig.7.Development of functional genetic markers for the mutant genes of soybean LNK2.(A).The functional genetic markers distinguishing the quadruple mutant and WT were designed from the edited sites of soybean LNK2 genes.Horizontal arrows indicate the molecular marker primes.Continuous lines indicate the nucleotide sequences of LNK2 genes.Dashes indicate deleted nucleotides.Nucleotide sequences indicate mutant regions;nucle otides in red indicate deleted bases,and underlined nucleotides indicate restriction endonuclease Mbo I and Msp I recognition sites.Green bases in the primer dCAPS2c-F indicate mismatched bases.(B) Phenotypic differences of the functional genetic markers of soybean LNK2 alleles between the quadruple mutant and WT.Boxes indicate the LNK2 genotypes of the quadruple mutant and WT tested by the developed markers.

4.Discussion

CRISPR/Cas9 technology is a powerful tool for rapid genome editing for crop improvement [30–32].It has been successfully used for gene editing in plant species includingArabidopsis,tomato,rice,wheat,and soybean [19,35,41–47].For instance,Cai et al.[19] obtained homozygous soybeanft2amutants using the CRISPR/Cas9 system and stable soybean transformation and found that the homozygous soybeanft2amutants had delayed flowering under both LD and SD conditions.CRISPR-Cas9-mediated targeted disruption ofFAD2-2caused mutations in the target sequence ofFAD2-2,and thefad2-2mutants showed both a high oleic acid content and a low linoleic acid content in soybean oil[48].In a recent study,CRISPR/Cas9-inducedGmPRR37-ZGDD mutants in soybean showed early flowering under natural LD conditions [49].In the present study,we designed three synthetic single-guide RNAs(sg1,sg2,and sg3) to edit fourLNK2genes.sg1 precisely targeted theLNK2aandLNK2bgenes,sg2 targeted theLNK2candLNK2dgenes,and sg3 targetedLNK2dbut did not work forLNK2c.We speculate that the reason why sg3 does not work forLNK2cis that we had few transgenic soybean lines,and sg3 has low editing efficiency forLNK2c.

In soybean,several QTL for agronomic traits,including plant height,flowering time,and salinity tolerance,have been cloned[11,14,50,51].CRISPR/Cas9-mediated multiplex gene editing provides a means of simultaneously manipulating several agronomic traits.Recently,multiplex gene editing has been reported for various crops[35,43].Miao et al.[43]used CRISPR/Cas9 technology to edit rice group I (PYL1–PYL6andPYL12) and group II (PYL7–PYL11andPYL13)PYLgenes and found that group I played more important roles than group II in stomatal movement,seed dormancy,and growth regulation.However,because the efficiency of stable transformation in soybean is lower than that in other crops,it is difficult to conduct multiplex gene editing.In the present study,we identified four LNK2 homeologs in soybean by sequence alignment with legume and other plant species.We then performed stable soybean transformation and obtained 12 T0transgenic lines.Although there is no scientific evidence that transgenic plants are harmful for human health,the debate on the safety of transgenic organisms continues.Accordingly,we generated two T1transgene-free homozygous quadruple mutants ofLNK2by selfing.The chance of producing a homozygous quadruple mutant ofLNK2was 1/256 (0.4%) from a heterozygous mutant ofLNK2,but the actual rate at which mutants appeared was 11.8%(2/17).We speculate that transgene-encoded Cas9 and sgRNAs again induced double-strand breaks inLNK2genes as the transgenic plants produced seeds.Our findings establish that the CRISPR/Cas9 technology offers great potential for soybean breeding.

The importance of flowering time in the adaptation of soybean is well established,and several major genes have been identified[11,14,52].For example,J,the ortholog ofArabidopsis thaliana EARLY FLOWERING3(ELF3),which is a major classical locus conferring the long-juvenile (LJ)trait,promoted flowering and improved adaptation of soybean to the tropics under short day lengths [11].Recently,we identified two homeologous pseudo-responseregulator genes,Tof11andTof12,and showed that they act viaLHYhomologs to promote the expression of theE1gene and to delay flowering under LD conditions [53].In the present study,the flowering and pod production times were significantly shorter in the quadruple mutants than in WT plants under LD conditions.However,plant height and node number did not differ between the quadruple mutant and WT plants under LD conditions.LNK2 promoted the expression ofE1and repressed that ofFT2ain soybean.Thus,LNK2 plays a key role in flowering time but does not affect plant height.Our results revealed an important new component in the flowering-time control of soybean.

5.Conclusions

Application of a CRISPR/Cas9 system to edit four soybeanLNK2genes yielded a transgene-free homozygous quadruple mutant ofLNK2.The mutant showed accelerated flowering under long photoperiod.These findings provide insight into the mechanisms underlying plant flowering time regulatory networks in soybean.

CRediT authorship contribution statement

Lidong Dong:wrote the manuscript,designed the experiments and managed the projects.Fanjiang Kong:designed the experiments and managed the projects.Baohui Liu:designed the experiments and managed the projects.Zhaobo Li:wrote the manuscript,performed experiments.Qun Cheng:wrote the manuscript,performed experiments.Zhuoran Gan:performed experiments.Zhihong Hou:performed experiments.Yuhang Zhang:performed experiments.Yongli Li:modified and revised the manuscript.Haiyang Li:performed data analysis.Haiyang Nan:performed data analysis.Cen Yang:performed data analysis.Linnan Chen:performed data analysis.Sijia Lu:performed data analysis.Wenqian Shi:performed data analysis.Liyu Chen:performed data analysis.Yanping Wang:performed data analysis.Chao Fang:performed data analysis.Liping Kong:performed data analysis.Tong Su:performed data analysis.Shichen Li:performed data analysis.Kun Kou:performed data analysis.Lingshuang Wang:performed data analysis.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work was supported by National Key Research and Development Program of China (2017YFD0101305),and the National Natural Science Foundation of China (31930083,31901568,31801384,31725021,and 31771815).

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2020.09.005.