Jing Yng, Zhenhu Guo,b, Lixin Luo, Qioli Go, Wuming Xio, Jifeng Wng, Hui Wng,Zhiqing Chen,, To Guo,
a National Engineering Research Center of Plant Space Breeding, South China Agricultural University, Guangzhou 510642, Guangdong, China
b Rice Research Institute of Heilongjiang Academy of Agricultural Sciences, Jiamusi 154026, Heilongjiang, China
Keywords:
ABSTRACT Early seedling vigor (ESV) is a major breeding target in rice, especially under direct seeding.To identify quantitative trait locus (QTL) affecting ESV, a recombinant inbred line population derived from a cross between 02428 and YZX, two cultivars differing in vigor during early seedling growth, was used for QTL analysis.Nine traits associated with ESV were examined using a high-density map.Of 16 additive loci identified,three were detected in two generations and thus considered stable.Four epistatic interactions were detected, one of which was repeated in two generations.Further analysis of the pyramiding effect of the three stable QTL showed that the phenotypic value could be effectively improved with an increasing number of QTL.These results were combined with results from our previous QTL analysis of the germination index.The lines G58 and G182 combined all the favourable alleles of all three stable QTL for ESV and three QTL for germination speed.These two lines showed rapid germination and strong ESV.A total of 37 candidate differentially expressed genes were obtained from the regions of the three stable QTL by analysis of the dynamic transcriptomic expression profile during the seedling growth period of the two parents.The QTL are targets for ESV breeding and the candidate genes await functional validation.This study provides a theoretical basis and a genetic resource for the breeding of directseeded rice.
Rice (Oryza sativaL.) serves as a staple food source for more than half of the global population [1].With the development of the social economy, rice cultivation and production systems have changed gradually.In recent years, direct seeding of rice has become increasingly popular because of its low cost and simplicity[2].Direct-seeded rice (DSR) is classified mainly as wet, dry, or water DSR [3].However, in practice, many factors restrict the development of DSR.For example, in water DSR, rice seeds are completely submerged, leading to poor or no germination and a low emergence rate [4].In wet or dry DSR, there is no standing layer of water to inhibit the growth of weeds,leading to high weed infestation and,consequently,increased yield loss[5–7].Seedlings with rapid uniform emergence and high vigor can access nutrients effectively and suppress the growth of weeds[8–10].In particular,high early seedling vigor (ESV) is crucial for the establishment of seedlings and their eventual biomass production or yield [11].It is thus desirable to breed rice cultivars with high ESV.
Quantitative trait locus (QTL) mapping is a main tool used to identify genetic factors affecting agronomic traits.Numerous QTL associated with seedling vigor have been detected [11–19].For example, Zhou et al.[13]used recombinant inbred lines (RILs)and composite interval mapping to identify nine QTL for ESV traits that correlated positively with one another.The individual QTL explained 4%–14%of the total trait variation.Xie et al.[17]mapped eight QTL for ESV using a RIL population and narrowed two major QTL,qSV-1andqSV-5c, to intervals of 1.13 Mb and 400 kb, respectively.A major QTL for seedling height,qPHS3-2,was fine-mapped.Within the mapped region, the geneOsGA20ox1, involved in gibberellin (GA) biosynthesis, was considered the most likely candidate gene [15].Zhang et al.[18]used a RIL population derived from 93–11 and PA64s and identified 57 QTL for ESV.With a BC3F2population derived from 93–11 and a chromosome segment substitution line (CSSL) harbouring segments from PA64s in the 93–11 background, a major QTL for seedling shoot length,qSSL1b,was fine-mapped within 80.5 kb.Recently, He et al.[20]reported that disruption of an isopropylmalate synthase gene,OsIPMS1,resulted in low germination speed and seedling growth under various conditions.Further analysis confirmed that regulation of seed vigor byOsIPMS1is associated with starch hydrolysis, glycolytic activity, and energy levels in germinating seeds.Although these studies provide us with some valuable information, the data remain insufficient.
High-density single-nucleotide polymorphisms (SNPs) are becoming more readily available owing to the development of sequencing technology [21].Genetic maps based on high-density SNPs show increased QTL mapping resolution, thus helping improve the efficiency of mapping [22].
In our previous study, we investigated the ESV of 200 conventional rice varieties, among which 02428 has high ESV, and YZX has medium ESV.This indicates that 02428 may carry favourable ESV genes.Therefore, in this study, we used a high-generation RIL population consisting of 192 RILs derived by single-seed descent from theindicacultivar YZX andjaponicacultivar 02428.The genetic map of this population consisted of 2711 bin markers, which were obtained via genotyping-by-sequencing (GBS)[23].The objectives were, using QTL mapping, to identify further QTL for ESV and, using transcriptome expression profiling, to identify differentially expressed genes as candidates for the QTL.Overall, these results broaden our understanding of the genetic and molecular bases of ESV in rice.In particular, some QTL are strong candidates as promising targets for markerassisted breeding of DSR varieties.The identified DEGs also provide valuable information for candidate gene verification and the dissection of gene regulatory networks affecting early seedling growth in rice.
The mapping population contained 192 RILs derived by singleseed descent from a cross between theindicaYZX andjaponica02428 cultivars.RILs from 6th- and 7th-generation populations were planted in respectively the wet and dry seasons of 2016,and the seeds were harvested as experimental materials for each season, and related characteristics were investigated after treatment.
The planting, management, harvesting and storage of the materials followed methods previously reported [24].The RILs and parents were planted in the experimental field of South China Agricultural University, Guangzhou, China (23.16°N,113.36°E).Twenty-five-day-old seedlings were transplanted into the experimental field with plant and row spacings both of 20 cm.Plants with the same heading date were selected, and the days to heading of each RIL was estimated; heading date was considered the timepoint when the leading panicle emerged from the leaf sheath to 1 cm.To synchronize seed maturation,six individual plants in the middle of each block were independently harvested on the 35th day after heading in the wet season and on the 40th day after heading in the dry season.The harvested seeds were dried at 42 °C for 5 days and then stored at -20 °C.
Three of the six harvested individual plants were randomly selected, and their seeds were oven treated at 50 °C for 7 days to break dormancy.All seeds were sterilized by soaking in 20%diluted bleach(6%–7%NaOCl)for 30 min and then rinsed six times with sterile water.A total of 100–200 healthy and plump seeds were selected from each plant and placed in a 9-cm-diameter Petri dish onto two layers of filter paper.Then, 10 mL of sterilized distilled water was added, and the dishes were placed in a chamber that had an 8-h light (200 μmol m-2s-1)/16-h dark photoperiod at 30 °C.Approximately 2 days later, 10 seeds (whose radicle or germ length had reached approximately 1 mm) per plant were selected, placed in a germination box (length 19 cm, width 13 cm, height 12 cm) and covered with two layers of filter paper,after which 20 mL of sterilized distilled water was added.The culture conditions were as described above.After six days, a WinRHIZO root image analysis system (Regent Instruments Inc.,Québec,Canada)was used to measure the seedling stem diameter(SSD), root length (RL), root surface area (RSA), root volume (RV),and root diameter (RD).Root fresh weight (RFW) and shoot fresh weight (SFW) were measured on a balance, total fresh weight(TFW) was calculated as the sum of RFW and SFW, and seedling height (SH) was measured with a ruler.
We used 192 RILs derived from the 02428/YZX cross to map the QTL involved in early rice seeding growth.The RIL population was genotyped by GBS, the Nipponbare genome (IRGSP 1.0) was used as the reference genome for SNP detection, and 85,743 highquality, biallelic, homozygous SNP markers were obtained.Using a sliding window approach,a high-density linkage map containing a total of 2,711 bin markers was constructed.All bin markers were evenly distributed across 12 chromosomes, the average physical length between the markers was 137.68 kb,and the average genetic distance was 0.86 cM.Among these 12 chromosomes,chromosome 10 had the fewest markers,at 162,and chromosome 3 had the most markers,at 311.The total genetic distance of the genetic map was 2343.68 cM, where Chr.10 was the shortest chromosome, at 118.78 cM,and Chr.1 was the longest,at 232.15 cM(Fig.S1)[23].QTL for ESV were mapped using the BIP function in QTL IciMapping v4.1 software [25].The inclusive composite interval mapping(ICIM-ADD) method was used to more precisely identify QTL parameters as follows: a 1 cM window size, 500 permutations and a 0.05 type I error rate were employed to call QTL.The logarithm of odds (LOD) score threshold was set to 2.5 [26].Epistatic QTL were identified via the inclusive composite interval mapping of epistatic QTL (ICIM-EPI) method, with the LOD score threshold set to 5.0 [27].Genes were annotated and analysed using the Rice Annotation Project (RAP) (http://rapdb.dna.affrc.go.jp) and Ensembl(http://plants.ensembl.org/Oryza_sativa)databases.
All methods followed our previous report [28].RNA samples were reverse-transcribed to cDNA using the High-Capacity cDNA Archive Kit (Applied Biosystems, Foster City, California, USA).qRT-PCR was conducted using the AceQ qPCR SYBR Green Master Mix Kit (Vazyme Biotech) according to the standard protocol, and the expression levels of the genes were determined on a StepOne-Plus system(Applied Biosystems,USA).Three replicates were used for each treatment.As an endogenous control, actin was used to normalize the obtained Ct values, and the relative expression values were calculated by the ΔΔCt method.Gene-specific primers were designed using NCBI primer BLAST (http://www.ncbi.nlm.nih.gov/tools/primer-blast/).The primer sequences of the four candidate genes are listed in Table S1.
Analysis of variance, correlation analyses,t-tests and Duncan’s multiple comparisons were performed using Statistical Analysis System (SAS, version 8.01) and Microsoft Excel.
In this study, we investigated the dynamic phenotypic changes in the 9 traits in two parents(Fig.1).By day 6,the differences in all the traits between the two parents were significant or extremely significant.In particular, SH showed extremely significant differences at all time points; RL and RSA were significant or extremely significant on days 3, 4, 5, and 6; RV was extremely significantly different on days 4, 5, and 6; RFW and TFW were significant or extremely significant on days 2, 4, 5, and 6; and SFW was significant or extremely significant on days 4, 5, and 6.However, RD and SSD became significant or extremely significant only on day 6.In addition, due to lateral root growth, the RD of both materials showed an increase followed by a decrease, while SSD showed a slow increase, and the remaining seven traits showed rapid increases.Notably, on day 6, the SH and RD of YZX were significantly higher than those of 02428, while the remaining 7 traits were higher for 02428 than for YZX.In general, the growth trends between 02428 and YZX were quite different, indicating substantial genetic differences between the two parents, which is beneficial for the identification of QTL.
All nine traits varied widely in the RIL population (Table 1, S2;Fig.S2).In the 6th generation,the coefficient of variation(CV)ranged from 10.25% to 26.44%, and in the 7th generation, it ranged from 10.48% to 28.04%.In particular, the traits showed differing amounts of transgressive segregation.Continuous single-peak distributions were observed for all traits, and the absolute skewness and kurtosis values of each trait were very close to 0, indicating that all of these traits were nearly normally distributed and suggesting that they were controlled by multiple minor-effect genes.
Fig.2 shows the pairwise phenotypic correlations among the nine traits.The results for the 6th and 7th generations were very similar.There were no strong correlations between aboveground and belowground traits, and the correlation coefficient between SFW and RFW was the largest: 0.62 in the 6th generation and 0.63 in the 7th generation.There was almost no correlation between SH and SSD among the three aboveground traits.The correlation coefficients SFW and SSD and between SFW and SH were similar at approximately 0.60.This result indicated that the contribution rates of SSD and SH to the SFW of early seedling shoots were not very different.Most of the five underground traits showed high pairwise correlations.The correlation between RSA and RL was the highest at 0.93 in both generations.Interestingly,RD showed a negative correlation with RL, with - 0.48 in the 6th and- 0.47 in the 7th generation.This finding was consistent with the change in the root system between 1 and 6 days observed for both parents;that is,with the growth of lateral roots,the diameter gradually decreased (Fig.1b).
QTL analysis revealed 29 QTL (Table 2; Fig.3).Of the 29, 14 were identified in the 6th and 15 in the 7th generation.The QTL were distributed on all chromosomes except 5,7,and 12.The phenotypic variation explained (PVE) ranged from 0.74% to 16.64%.QTL were detected for all nine traits, with the maximum number of QTL (five) associated with SFW and RL and only two QTL each with RD, RFW, SH, and SSD.QTL with the same physical position were considered identical.Of the final total of 16 loci, three were repeatedly detected in the two generations: locus 1 (qRL-1, 6th and 7th), locus 3 (qRSA-3, 6th and 7th;qRV-3, 6th and 7th;qRFW-3,7th;qTFW-3, 7th), and locus 13 (qSSD-9, 6th and 7th;qSFW-9, 6th;qRFW-9, 6th;qTFW-9, 6th).Locus 9, 12, and 16 were associated with two traits in a single generation.Of the identified loci,four have been reported previously [12,29–31](Table 2).This finding also confirms that our results are reliable.
Four epistatic interactions were detected across two generations for three traits(one for SSD,two for SH,and one for RD),distributed on chromosomes 2, 3, 4,6, 9, and 10 (Table 3).These QTL did not overlap additive QTL.The PVE of these epistatic interactions ranged from 9.25% to 20.00%.One pair of epistatic QTL for SH was detected in both the 6th and 7th generations,and the epistatic effect was negative, indicating that the recombinant phenotypic value was larger than the parental one.In summary,additive×additive epistatic effect is also an important genetic mechanism for ESV in terms of SH, SSD, and RD.
Among the 16 loci, three, namely locus 1, 3, and 13, were detected in both generations and thus considered stable.To further validate the effects of the three QTL,we first divided the individuals of the RIL population into the 02428 type (marker type 2) and YZX type (marker type 0) for each identified locus and then assessed the differences in the corresponding traits between the two genotypes (Table 4).The mean phenotypic values for all of the corresponding traits of RILs harboring favourable alleles were significantly (P<0.01) greater than those of RILs harboring unfavorable alleles.Some traits showed QTL in only a single generation,but extremely significant differences were observed between 0-type RILs and 2- type RILs in both generations.This finding suggested that these three stable QTL are credible and could increase the values of the corresponding traits.
To further confirm the three stable QTL with effects on ESV,without considering the effects of interactions among or environmental influences on these QTL.We performed pyramiding analysis of these three QTL, dividing into eight classes the recombinant types of the three QTL.We investigated the distribution of QTL in individual RILs based on the relationship between phenotypic values and the number of favorable QTL alleles present (Fig.4;Table S3).Although these three QTL did not affect all traits, most of the traits still showed an influence of other additive QTL.All phenotypic values but those of RD increased with increasing aggregation of favorable alleles.In summary, this result indicated that the pyramiding of favourable alleles could improve ESV in the varieties.
Both rapid, uniform germination and vigorous seedling growth are important characteristics of DSR cultivars.In our previous studies,we used the same population to identify QTL for seed germination ability at a normal temperature(30°C).Three QTL,qGI-6,qGI-7,andqGI-9,associated with germination index(GI)with large PVE and stable expression in the 6th and 7th generations were obtained[28].To identify lines with rapid and uniform seed germination and high ESV,we analyzed the genotypes of ESV and GI-associated QTL in this population.According to the source of the favorable allele,pyramiding of the favourable alleles of three QTL for ESV(locus 1,3, and 13) and three QTL for seed germination speed (qGI-6,qGI-7, andqGI-9) was detected in two RILs (G58 and G182)(Fig.5a).Taking the phenotypes of the 7th generation as a reference, we found that most values of the 10 traits of G58 and G182,including GI,were higher than the low values of the parents and the mean values of the population,and even transgressive segregation occurred, with a value higher than the high value of the parents.The phenotypes were improved to some extent (Fig.5b).Thus,pyramiding of the ESV-and GI-associated QTL could increase the rapid germination of seeds and ESV.Overall, G58 and G182 could serve as favorable allele donors in breeding and increase the efficiency of pyramid breeding of DSR.
Table 1 Phenotypic analysis of the YZX, 02428, and YZX× 02428 RIL population.
Fig.2.Correlations between nine traits associated with ESV in the RIL population.(a)6th generation.(b)7th generation.In the upper diagonal,the size of the circle and the intensity of shading indicate the magnitude of the correlation.The negative correlations are colored green and the positive correlations red.The lower diagonal shows the correlation coefficients, and * and ** represent significance at the 0.05 and 0.01 probability levels, respectively.
Markers flanking the QTL were used to anchor them to the reference genome IRGSP-1.0.A total of 259 genes were identified in the QTL regions(Table S4).The list of candidate genes was shortened by analysis of previously completed transcriptomic expression profiles[28].RNA-seq was performed using RNA extracted from all tissues(seed+seedling)of germinated YZX and 02428 seeds on days 0,2,3,and 4.Considering that the early seedlings of 02428 and YZX all grew rapidly within 1 to 6 days,the difference was obvious.Differential expression profiling (with false discovery rate ≤0.05 and absolute value fold change ≥2.0) of the two parents was accordingly performed on days 2, 3, and 4, and respectively 5128, 5421,and 5117 DEGs were detected (Fig.S3).The combination of QTL mapping and expression profiling yielded 34 candidate DEGs for the three stable QTL (Fig.6; Table S5).Of these DEGs, eight were from locus 1,11 from locus 3,and 15 from locus 13.
Table 2 QTL associated with ESV in two generations.
Fig.3.All QTL positions in the high-density map.The red symbols represent QTL detected in the 6th generation and the blue patterns those detected in the 7th generation.The words associated with different shapes are abbreviations of different phenotypes.The circles represent QTL detected in both the 6th and 7th generations.
To verify the accuracy of our expression profiling results, we further selected four genes for qRT-PCR according to gene expression level (Fig.7).They included anexpressed protein(Os01g0655500), anOsWAK26-OsWAK receptor-like protein kinase(Os03g0643250),aflavin-containing monooxygenase family protein(Os09g0549300), and anHMGd1 protein(Os09g0551600).Transcriptome sequencing of the two materials was not performed on days 5 and 6.Accordingly, qRT-PCR of the four genes was performed on days 4, 5, and 6.On these days, the four genes differed significantly in expression between the two parents, suggesting that they accounted for parental differences in early seedling growth.
Table 3 Digenic epistatic QTL in the RIL population identified by interval mapping.
Table 4 Summary of the phenotypic effects of three stable QTL.
A rice cultivar with strong ESV trait is desirable for direct seeding [3].For example, in the water DSR system, seeds are broadcasted into standing water, and rapidly growing buds help make contact with oxygen as soon as possible.A healthy and vigorous root system can help anchor seedlings in soil, thus reducing the population imbalance caused by dead or floating seedlings.In dry and wet DSR systems, one of the major threats is weedy rice[32,33].Robust seedlings can have a competitive advantage over weeds, and rapid growth of seedlings can allow fields to be irrigated as soon as possible and create a hypoxic environment,inhibiting the growth of weeds.However, the traits previously examined by researchers were not comprehensive enough, especially the indexes of the root system, which were mainly root weight and maximum root length [11,12,18].The period over which seedling vigor was measured also varied.In the present study, we collected phenotypic data six days after material treatment,which is a suitable irrigation period in production practices.Use of a root scanner to acquire RL,RD,RSA,RV,SSD,and SH values and a balance to measure SFW, RFW, and TFW helped us evaluate ESV comprehensively and improved the efficiency of QTL mapping.
In general,two parents with large differences in target traits are used to construct populations, facilitating the detection of QTL for the traits[34].In this study,the mapping population derived from 02428 and YZX showed extreme variation in ESV, and thus was suitable for QTL identification.
Fig.4.Pyramiding analysis of three stable QTL.(a)Summary of the traits associated with the three QTL.A yellow background indicates that the locus is associated with the trait.(b) Pyramiding effects for different numbers of favorable alleles of the QTL.Letters from a to d indicate significantly different values according to statistical analysis using Duncan’s multiple range test (α = 0.05).
In conventional QTL mapping, molecular-marker genotyping is time-consuming and labor-intensive [35].The markers are mostly distributed at low density and are not able to provide precise and complete information about the numbers and locations of QTL for the corresponding traits[36].In the present study,we used a highresolution genetic map containing 2711 bin markers.Compared with previous QTL intervals associated with ESV, the intervals in this study were much narrower and the physical intervals containing the three stable QTL were only approximately 450 kb on average.Another advantage of our study is the integration of QTL mapping and profiling analysis.Numerous studies[16,37–40]have suggested that combining QTL mapping and profiling analysis is a powerful method for successfully reducing the number of candidate genes.We identified 259 genes in the intervals and flanking sequences of the three loci (Table S4), and further combination of these genes with expression profiling data successfully reduced the number of candidate genes to 34.
Finally, we obtained the four most promising candidates.The first candidate isexpressed proteingene (Os01g0655500), whose expression increased with the growth of seedlings in both parents.Moreover,Os01g0655500showed significant or extremely significant differences in expression between the two parents at all time points(Table S4;Fig.7).The second candidate isOsWAK26-OsWAK receptor-like protein kinasegene (Os03g0643250), which belongs to the wall-associated kinase(WAK)gene family,and plays important roles in cell division or growth inArabidopsis thaliana[41].The third candidate is aflavin-containing monooxygenase family proteingene(Os09g0549300).YUCCAencodes a flavin monooxygenase–like enzyme inArabidopsis thaliana, which catalyses the rate-limiting step in auxin biosynthesis from tryptamine toN-hydroxyl tryptamine, and plays an important role in auxin biosynthesis [42].Therefore, we speculate thatOs03g0643250andOs09g0549300may play a role in the early growth of rice seedlings similar to that observed inArabidopsis thaliana.The last candidate is anHMGd1protein gene (Os09g0551600).The expression level ofOs09g0551600in 02428 increased with seedlings growth, which was consistent with the gene expression atlas of Minghui 63 rice covering the entire life cycle[43];that is,the gene expression level was high in the seedling stage.However, the expression level of this gene in YZX was almost unchanged (Table S4).These genes provide a foundation for understanding the mechanism of early seedling growth in rice.
Fig.5.Pyramiding of loci for ESV and GI in lines G58 and G182 with rapid germination and strong ESV.(a)Genetic background of G58 and G182.(b)Phenotypes of G58,G182,the population mean, and the two parents in the 7th generation.
Rice cultivars have been bred specifically for transplanting systems, and ESV has not been selected in conventional breeding for crop improvement programmes [3].In this study, we identified 16 loci, of which three showed repeatability.Some showed high PVE, for example locus 9 (qTFW-6,qSFW-6–1, 7th) and 10 (qSFW-6–2,7th).These loci warrant further examination by fine mapping and candidate gene analysis for ultimate application in markerassisted selection.Many studies [44–48]have shown that marker-assisted pyramiding of QTL is an effective method for improving traits,particularly quantitative or complex traits.In this study, we tested the pyramiding effect of the three stable QTL in the population by combining favorable alleles.
As expected,aggregating favorable alleles improved the phenotypic traits(Fig.4).However,aggregation is not the same as simple superposition.We must consider the influence of the environment on genetic locus effects,and epistatic interactions between QTL are also an important factor.In particular, in our study, four epistatic interactions were detected with high PVE (Table 3).Interestingly,our pyramiding analysis, which included two RILs (G58 and G182) with pyramiding of three rapid-germination QTL alleles,showed phenotypes consistent with our expectations.One of the main reasons may be that the digenic interactions detected in the present study did not involve additive QTL.The actual QTL pyramiding effect awaits further evaluation by marker-assisted selection of QTL in rice accessions with the same genetic background.
Fig.6.Candidate gene identification by QTL mapping and expression profiling.(a) Three stable QTL located on chromosomes 1, 3, and 9 identified by QTL analysis.(b)Detection of DEGs by expression profiling on days 2, 3, and 4.Red and green dots represent respectively up- and downregulated genes.Gray dots represent genes with no change in expression.
Fig.7.Relative expression levels of four DEGs.Error bar indicates the standard deviation of mean.* and ** represent significance at the 0.05 and 0.01 probability levels,respectively, according to a t-test.
CRediT authorship contribution statement
Tao Guo,Zhiqiang Chen,and Jing Yang designed the project,and Jing Yang performed all the experiments andwrote the manuscript.Zhenhua Guo, Lixin Luo, Qiaoli Gao, Wuming Xiao, Jiafeng Wang,and Hui Wangassisted in conducting the experiments and analyzing the data.Zhiqiang Chen and Hui Wang provided thedirection for the study and corrections to the manuscript.All authors read and approved the final manuscript.
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 research was supported by the Breeding New Varieties of Rice Suitable for Light and Simple Cultivation and Mechanized Production Project (2017YFD0100104), the Research and Development Plan for Key Areas in Guangdong Province(2018B020206002), and the China Agriculture Research System(CARS-01-17).Special thanks are due to the South China Agricultural University Doctoral Innovative Talents (Domestic Training)Cultivation Program (CX2019N044).The funding bodies had no role in the study design, data collection, analysis and interpretation, decision to publish, or writing of the manuscript.
Appendix A.Supplementary data
Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2020.08.010.