Update of Genetic Linkage Map and QTL Analysis for Growth Traits in Eucommia ulmoides Oliver

: Eucommia ulmoides (Tu-chung) is an economically and ecologically important tree species which has attracted worldwide attention due to its application in pharmacology, landscaping, wind sheltering and sand ﬁxation. Molecular marker technologies can elucidate the genetic mechanism and substantially improve the breeding e ﬃ ciency of E. ulmoides . The current research updated the original linkage map, and quantitative trait loci (QTL) analysis was performed on tree growth traits measured over 10 consecutive years in an E. ulmoides F1 population (“Xiaoye” × “Qinzhong No.1”). In total, 452 polymorphic markers were scored from 365 simple sequence repeat (SSR) primers, with an average of 1.24 polymorphic markers per primer combination. The integrated map was 1913.29 cM (centimorgan) long, covering 94.10% of the estimated genome and with an average marker density of 2.20 cM. A total of 869 markers were mapped into 19 major independent linkage groups. Growth-related traits measured over 10 consecutive years showed a signiﬁcant correlation, and 89 hypothetical QTLs were forecasted and divided into 27 distinct loci. Three traits for tree height, ground diameter and crown diameter detected 25 QTLs (13 loci), 32 QTLs (17 loci) and 15 QTLs (10 loci), respectively. Based on BLASTX search results in the NCBI database, six candidate genes were obtained. It is important to explore the growth-related genetic mechanism and lay the foundation for the genetic improvement of E. ulmoides at the molecular level. have and the


Introduction
As a relic plant that experienced the third glacial period, Eucommia ulmoides is a dioecious perennial deciduous tree (2n = 34) and the only surviving species of the genus Eucommia (Eucommiaceae) [1,2]. Historical documents show that E. ulmoides is naturally distributed in central China, protected by complex terrains [1], and it has been considered highly adaptable to various site conditions. Since its successful introduction into France in 1896, E. ulmoides has been introduced into Korea, Japan, Germany, Russia and the United States [3]. E. ulmoides bark has been used as traditional Chinese medicine for thousands of years [4,5]. Modern studies have revealed that E. ulmoides has medicinal and health effects, such as lowering blood glucose and blood lipid levels and regulating blood pressure [6,7]. The gutta-percha in E. ulmoides leaves, fruits and phloem tissues can be a substitute for rubber, exhibit good resistance to strong temperature and acid-base changes, and is also an excellent insulator [8]. The species is also used in landscaping and forestry because of its straight trunk, flourishing leaves, wind sheltering, sand fixation and positive impacts on water and soil conservation. The two parents were chosen because of the differences in leaf morphology characteristics, secondary metabolite content and climatic and geographic conditions of the origin areas. A wild phenotype "Xiaoye" was selected as the female parent, originating from the forest at Yantuo, Lingbao, Henan, with the following characteristics: small leaves, smooth bark, low contents of secondary metabolites and late budding and flowering times (convenient for controlled pollination). "Qinzhong No. 1" was selected as the male parent with a superior breeding variety, produced through selective asexual breeding and planted in the museum garden of Northwest A&F University, Yangling, Shaanxi (34 • 16 1" N, 108 • 4 21" E, 460 m above sea level), with the following characteristics: large leaves, rough bark, high levels of secondary metabolites and early budding and flowering times. "Qinzhong No. 1" also displayed high resistance to drought and cold, with a high growth rate [11].

DNA Extraction
We isolated DNA from fresh, young, fully expanded leaves of 152 F1 full sibs and two parental trees ("Xiaoye" and "Qinzhong No. 1"), according to a modified cetyltrimethylammonium bromide (CTAB) method [24]. The integrity of extracted DNA was checked via 1.2% agarose gel electrophoresis stained by ethidium bromide (EB), and the concentration was estimated with a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies Inc, Wilmington, Delaware, USA).

Simple Sequence Repeat (SSR) Analysis
Based on transcriptome data, 2200 pairs of SSR primers were randomly selected and synthesized by GENEWIZ Inc. (Suzhou, China) [22]. Eight samples (six F1 individuals and two parental trees) were used to assess the polymorphism and practicability of the SSR markers. A total of 365 SSR primers were validated as polymorphisms in the DZ0901 population [22,25] (Table S1) and were used for PCR amplification with 152 F1 individuals and two parents. The PCR amplification was carried out in 25 µL of reaction mix containing DNA (approximately 50 ng), dNTPs (0.2 mM), MgCl 2 (2.5 mM), each primer (0.4 mM), PCR (polymerase chain reaction) buffer and Taq DNA polymerase (1.5 U). The amplification protocol of SSR-PCR was as follows: 3 min at 94 • C; 30 s at 94 • C, 30 s at the locus-specific annealing temperature, 60 s at 72 • C for 35 cycles; 10 min at 72 • C (final extension step). The results of the PCR amplification were detected by polyacrylamide gel electrophoresis on 8% of non-denaturing polyacrylamide gel and then visualised through AgNO 3 staining. Polymorphic markers were recorded (repeat three times) as the presence or absence on the mapping population and then converted to "ll × l m" (heterozygous in female parent), "nn × np" (heterozygous in male parent) and "hk × hk" (heterozygous in both parent) to fit the map construction software JoinMap 4.0 (Plant Research International B.V. and Kyazma B.V., Wageningen, Gelderland, The Netherlands) [26].

Segregation Analysis and Linkage Map Construction
The E. ulmoides linkage map was updated with the pseudo-testcross mapping strategy using JoinMap 4.0 [26] combined with previous data [21]. A chi-squared test was performed to detect the segregation distortion with recorded SSR markers. Markers fitted to the expected Mendelian segregation ratio were applied for linkage analysis with the Kosambi mapping function and grouped following these parameters: 4.0 for the logarithm of odds (LOD) threshold, 0.45 for the recombination fraction threshold, 1.0 for the ripple value and 5.0 for the jump threshold. Linkage groups were drawn using the MapChart 2.2 software (Plant Research International B.V. and Kyazma B.V., Wageningen, Gelderland, The Netherlands) [27]. The arithmetic for estimating genome length (G 2 ) was proposed by Chakravarti et al. [28] in 1991: G 2 = G 1 × (m + 1)/(m − 1); G 1 refers to the observed length of one linkage group and m refers to the number of markers mapped on the linkage group. The computing method of the constructed linkage map coverage was assessed by G 1 sum/G 2 sum.

Growth Traits Assessment
Three vital growth traits (H, tree height; GD, ground diameter at 20 cm above ground level; C, crown diameter) were measured to assess the growth condition of the F1 progenies in October from 2010-2019 (at the end of the growing season). The 4-year (2010-2013, 1a-4a) phenotypic data for H and GD have been described in a previous study [21]. Here, 1 or 1a refers to 2010 traits, 2 or 2a refers to 2011 traits, and so on. Analyses of descriptive statistics and Pearson correlations of traits were carried out applying the statistical software SPSS 18.0 (Inc. SPSS, Prentice Hall, Upper Saddle River, New Jersey, USA, 2010) for Windows. Violin Plots (Omicshare Tools, Guangzhou Gene Denovo Biotechnology Co., Ltd., Guangzhou, Guangdong, China) were used to show the distribution and probability density for tree height (cm, H), ground diameter (mm, GD) and crown diameter (cm, C) of the F1 population DZ0901 over 10 years. The following formula was used to calculate the average growth rate (GR) of three traits over a 10-year period. Here, GV refers to the growth value, a and b refer to the measuring time (year), and a was larger than b.

Quantitative Trait Loci (QTL) Analysis
The QTL analyses were performed using the software MapQTL 5.0 (Plant Research International B.V. and Kyazma B.V., Wageningen, Gelderland, The Netherlands) [29]. The K-W (Kruskal-Wallis) analyses (nonparametric test) were performed to estimate the marker significance level and the marker-phenotype association. The interval mapping model (IM) was used to estimate the location of a supposed QTL. The multiple-QTL mapping (MQM) model was employed with the cofactors detected by the IM model. The significant logarithm of odds (LOD) threshold (p < 0.05) was determined via a permutation test (PT) with 1000 iterations. However, those QTLs with an LOD score ≥3 and lower than the LOD threshold from PT were also recorded. The phenotypic variations explained by the identified QTLs (R 2 ) were calculated by variance analysis. The software MapChart 2.2 was used for QTL representation [27]. The original sequence information of the SSR markers flanking the target trait QTLs were obtained [22], a BLASTX search was conducted in the NCBI database, functional annotation was performed and relevant information of candidate genes was summarised.

SSR Analysis
In total, 452 polymorphic SSR markers were obtained from 365 SSR primers, with an average of 1.24 polymorphic markers per primer combination. Among the 452 polymorphic SSR markers, 135 markers were identified as l m × ll markers (1:1 segregation ratio), 173 markers were identified as nn × np markers (1:1 segregation ratio), 144 markers were identified as hk × hk markers (3:1 segregation ratio), and the remaining 22 markers (4.87%) showed segregation distortion (p < 0.05) and were therefore excluded from mapping (Table S2). Combined with the previous data, a total of 1949 markers compliance with Mendelian segregation ratios put into use for the construction of the genetic linkage maps.

Genetics Linkage Map
The genetic linkage map of E. ulmoides consisted of 19 linkage groups with 869 segregating markers covering 1913.29 cM (Table 1, Figure 1, Figure S1). The maximum number of makers per linkage group was 151 (G1), the minimum number was 15 (G17), and the average number was 46. The total number of added SSR makers was 224, with a mean number of 11.79 per linkage group. The shortest length of each linkage group was 59.13 cM (G16), the longest was 172.21 cM (G3), and the average length was 100.70 cM. The mean genetic distance between adjacent markers was 2.20 cM. The estimated length of the E. ulmoides genome was 2033.151 cM. In this study, the constructed linkage map covered 94.10% of the estimated genome length. In addition, 286 markers were not linked to any of the existing linkage groups and were grouped in 88 small groups. One and three markers linked to G6 and G7, respectively, but could not be ordered. The remaining 790 markers did not link to any linkage group. cM. The estimated length of the E. ulmoides genome was 2033.151 cM. In this study, the constructed linkage map covered 94.10% of the estimated genome length. In addition, 286 markers were not linked to any of the existing linkage groups and were grouped in 88 small groups. One and three markers linked to G6 and G7, respectively, but could not be ordered. The remaining 790 markers did not link to any linkage group.

Growth Traits
The genetic variation of F1 population phenotypic traits was significant. The variation coefficients of tree height, ground diameter and crown diameter over consecutive 10 (8) years were 6.30-42.81%, 18.09-34.12% and 13.37-33.77%, respectively ( Table 2). These growth-related traits were obviously separated in the mapping population and normally distributed, which was suitable for QTL analysis ( Figure 2).

QTL Analysis
A total of 89 hypothetical QTLs were predicted, of which 30 were related to tree height, 38 were related to ground diameter and 21 were related to crown diameter. The above-mentioned QTLs included the results of the QTL analysis for the growth rates of three traits, among which 5, 6 and 6 QTLs were identified for growth rate of tree height, ground diameter and crown diameter, The white dots in the middle represents the median value. The black rectangle in the middle represents 25th to 75th percentiles range. The middle vertical line extended from the black rectangle represents the 95% confidence interval. The vertical axis of the shape outside the black rectangle represents the data dispersion, and the horizontal axis represents the distribution frequency. 1a represents for 2010, 2a represents for 2011, and so on.  (Table S3).

QTL Analysis
A total of 89 hypothetical QTLs were predicted, of which 30 were related to tree height, 38 were related to ground diameter and 21 were related to crown diameter. The above-mentioned QTLs included the results of the QTL analysis for the growth rates of three traits, among which 5, 6 and 6 QTLs were identified for growth rate of tree height, ground diameter and crown diameter, respectively (Table 3, Figure 3).  3 Genome-wide permutation test, * logarithm of the odds (LOD) value was significant at P < 0.05 based on 1000 genome-wide permutation tests. 4 Linkage Group. 5 logarithm of the odds (LOD) value at peak position. 6 Marker name nearest to the QTL position. 7 The proportion of phenotypic variance explained by the QTL. 8 Kruskal-Wallis significance level, given by the P values (* 0.1; ** 0.05; *** 0.01; **** 0.005, ***** 0.001, ****** 0.0005, ******* 0.0001). The hollow bars represents QTLs for growth rates. Thick bars and thin lines represent 1-LOD and 2-LOD confidence intervals of each QTL, respectively. QTLs name were described in Table 3. . QTL mapping on genetic linkage groups for growth traits in Eucommia ulmoides. Red, green and blue bars represent QTLs for height, ground diameter and crown diameter, respectively. The hollow bars represents QTLs for growth rates. Thick bars and thin lines represent 1-LOD and 2-LOD confidence intervals of each QTL, respectively. QTLs name were described in Table 3.   3 Genome-wide permutation test, * logarithm of the odds (LOD) value was significant at P < 0.05 based on 1000 genome-wide permutation tests. 4 Linkage Group. 5 logarithm of the odds (LOD) value at peak position. 6 Marker name nearest to the QTL position. 7 The proportion of phenotypic variance explained by the QTL. 8 Kruskal-Wallis significance level, given by the P values (* 0.1; ** 0.05; *** 0.01; **** 0.005, ***** 0.001, ****** 0.0005, ******* 0.0001).
We found 3 loci linked to 3 traits (Locus 2, 3 and 11), while 5 loci were detected for tree height and ground diameter (Locus 7, 15, 19, 24 and 27), 2 loci were detected for ground diameter and crown diameter (Locus 6 and 17), and no loci were detected for tree height and crown diameter ( Table 4). Based on the results of BLASTX, a total of six candidate genes were obtained. There was one QTL containing the KEGG (Kyoto Encyclopedia of Genes and Genomes) and GO (Gene Ontology) annotation information, respectively, and two failed to be mapped to the database (Table 5).

SSR Marker and Segregation Distortion
Interspecific SSR primers were used to construct genetic maps for numerous species, but E. ulmoides has no species of the same genus (Eucommia) as candidate SSR sources. Therefore, the development of polymorphic SSR markers by next-generation sequencing is economical, efficient and necessary. The markers used in this study were also obtained from the transcriptome database, and 4.87% segregation distortion (SD) were detected, which was lower than for Castanea sativa (sweet chesnut) [18], Populus spp. (poplar) [30], Juglans regia (walnut) [31], Ziziphus jujuba (jujube) [32] and Citrus clementina (clement pomelo) [33]. Segregation distortion is a common problem in population mapping, and a number of variables could have led to SD [34][35][36][37][38]. In the current study, there were also significant differences in parental origin and growth-related traits between "Xiaoye" and "Qinzhong No.1" [39], and the genetic difference between parents was also an important factor affecting SD. There were no three or four alleles in the SSR marker, which may be related to the conservatism of CDS (coding sequence) identified via transcriptome sequencing.

Genetic Linkage Map
The E. ulmoides genetic linkage map contained more information with added SSR markers. Marker density and genome coverage have improved. The number of linkage groups (the mapped marker number was over 15) in the E. ulmoides linkage map was updated from 13 to 19 [21], which was larger than the chromosome haploid number (17). In addition, we obtained some small linkage groups, the same as in Juglans regia [31]. These phenomena have been found in many species [16,18,31,33,40,41] and indicate that some gaps prevent connection between linkage groups belonging to the same chromosome, which needs to be optimised in the future by expanding mapping population, increasing parental combinations and adopting new markers. By contrast with the variety of mapping populations in crops [42,43], the growth cycle of trees was long, and the genetic offspring was difficult to be obtained. The F1 mapping population is the most common one, and abundant variation can be gained through interspecific hybridisation in many trees [18,33]. In E. ulmoides, parents with significant differences should be valued.

Growth Traits
The growth rates of traits in 5a (2014) and 6a (2015) changed from slow growth to rapid growth, which may be related to transplanting (2014, 5a). After transplanting, plant spacing increased, leaving more growth space for the trees. Hence, the F1 progenies showed a rapid increase in tree height, ground diameter and crown diameter after adapting to the new environment. However, the slow growth of 9a (2018) and 10a (2019), especially the stagnant growth of crown diameter, may have been caused by the limited growth space, indicating that the phenotype was decided by the interaction of genotype and environment [31]. The research objects in this experiment were decennial high trees with a massive root system. The human, financial and land resources needed for transplanting were the primary issues. The design of the nursery and seed orchard and the use of fields or woodlands are important issues in tree cultivation. As a long-term task, tree cultivation needs to consider the actual condition of tree species and to better plan the use of fields or woodlands to obtain reliable research results.

QTL Analysis
Growth-related characteristics are important economic indices of E. ulmoides. The QTLs identified for each trait were distributed in different linkage groups, suggesting that the trait may be regulated by multiple genes with large effect. Three traits (tree height, ground diameter and crown diameter) had a complex genetic background, which was also confirmed for Pistacia vera (pistachio) [40], J. regia [31], P. tremula (European aspen) [44], H. brasiliensis [16] and Camellia sinensis (tea tree) [45]. The phenotypic variation explained by each QTL had a wide range, and the sum of the phenotypic variation for one trait was more than 100%. Beavis [46] has proposed the "Beavis effect" in 1994, which stated that the mapping population size could affect the phenotypic variation accounted for QTLs. In addition, when the distance between one QTL and its co-marker is close, the MQM model will absorb the effect of another QTL, thus affecting the phenotypic variation result. This problem can be solved by increasing the population size and optimising the QTL analysis model.
The tree height-related QTLs of different years were identified in the same region, forming a total of 13 loci. The same phenomenon also appeared in traits of ground diameter (17 loci) and crown diameter (10 loci). Growth rate-related QTLs also appeared in the aggregation regions of QTLs, indicating that the QTL loci contained the major genes which have a significant regulatory effect on the growth traits. In this study, the frequent QTLs showed various phenotypic variations among different years, and no QTLs were detected for 10 consecutive years. Expression instability of growth-related QTLs among years has been observed in many species [47,48]. Conson et al. [16] measured tree height, number of whorls, and stem diameter at seven different ages in H. brasiliensis, while Yang et al. [49] measured height and diameter at breast height over 2 years in Pinus elliottii (pelliottii), and the QTLs detected at different times did not cluster in the same region. These findings indicate that variations could be possible because different genes are involved in genotype-environment interactions or QTL regulation according to the developmental stage [50].
The growth traits of trees are easily influenced by the surrounding environment [31], and with the growth and maturity of trees, the genetic regulation mechanisms may adjust accordingly [51].
In this study, the QTL loci identified in G1, G5, G8, G13, G18 and G19 were associated with two or three traits. The results of the phenotypic correlation analysis revealed plain positive correlations among tree height, ground diameter and crown diameter in E. ulmoides. The QTLs for different traits with high correlation tend to appear in clusters in H. brasiliensis [16], Populus (poplar) [15], P. elliottii [49], Prunus dulcis (almond) [19] and J. regia [31]. It was generally inferred that the gene may regulate and control the expression of multiple traits or that, QTLs/genes were linked closely [52,53].
The QTL detection should formally apply on genotypic value, where environmental variations have been estimated and subtracted from the phenotypic observations. However, no replicates were set in this study, and therefore, estimates of environmental variance were not allowed and phenotypic values cannot be corrected. In the further study, we will replicate the mapping population (via grafting or tissue culture) to grow in different fields with various site conditions, with the aim to better evaluate the environmental and genetic effects.

Candidate Genes for QTLs
This study relied on previous transcriptome data to obtain annotation information and did not expect to identify all candidate genes for QTLs affecting growth traits in E. ulmoides.
Gigantea protein (GI) was predicted to be related to H1 and GD1 traits. It is generally believed that GI promotes flowering under long days in a circadian clock-controlled flowering pathway [54]. It is also involved in numerous physiological processes, mediating rhythmic, phytochrome B signalling, carbohydrate metabolism and cold stress response [55]. Ding et al. have described the identification of the Populus (poplar) orthologs of GI and their critical role in short-day-induced growth cessation and concluded that GI controls seasonal growth cessation in Populus [56]. In our study, GI may have played an important role in the photoperiodic control of growth cessation in E. ulmoides.
The ABIL1 was predicted to be related to H1, GD1 and C8 traits in E. ulmoides. Based on a previous study, ABIL1 encodes a subunit of the WAVE (Wiskott-Aldrich syndrome verprolin homologous protein) complex and is required for the activation of the ARP2/3 complex, which is involved in the nucleation and branching of actin microfilaments [57]. Studies have shown that actin microfilaments are involved in several physiological processes in plants, including organelle movement [58], intracellular vesicle transport [59], plant stress resistance [60] and response to plant hormones [61]. In addition, actin is involved in the process of cell division, acting as an orbit for the spindle and assisting in its localization in the middle of mitosis. This indicates that ABIL1 is an important regulatory factor in the growth and development of E. ulmoides.
The filamentation temperatures-sensitive H (FtsH) gene encoded an ATP and Zn 2+ -dependent protein with ATPase activity, proteolytic activity and molecular chaperone activity [62]. It plays an important role in thermal shock, hyperpermeability, light stress, cold induction, disease and other stress conditions [63]. The m-AAA complexes consist of AtFtsH3 and/or AtFtsH10 in Arabidopsis thaliana (arabidopsis), participating in the maturation of the ribosomal subunit L32 by proteolysis. Hence, in the double FtsH3/10 mutants of A. thaliana, molecular regulatory mechanisms are impaired, including mitotic biogenesis, mitochondrial translation and the function of the OXPHOS (oxidative phosphorylation) complex [64]. It is speculated that FtsH is involved in regulating mitochondrial activity in E. ulmoides cells.
The eukaryotic translation initiation factor 4G (eIF4G) is a scaffold protein that organises the assembly of initiation factors needed to recruit the 40S ribosomal subunit to an mRNA [65]. The eIF4G expresses two isoforms (eIF4G, eIFiso4G) in plants. Chen et al. [66] have investigated the role of eIFiso4G in plant growth by using null mutants for the eIF4G isoforms in Arabidopsis. The eIFiso4G Forests 2020, 11, 311 14 of 18 loss-of-function mutants exhibited smaller cell, leaf, plant size and biomass accumulation, especially a reduction in chlorophyll levels, which correlated with its reduced photosynthetic activity. It is speculated that eIF(iso)4G affects the growth and development of E. ulmoides by regulating photosynthesis.
Although genes can potentially influence the growth and development of plants, further studies should be conducted to validate and prove the effects of genes in E. ulmoides on the growth-related traits.

Conclusions
In this study, the original genetic linkage maps were updated with transcriptome SSR data, and QTL analysis was performed on tree growth traits over 10 consecutive years in the E. ulmoides F1 population ("Xiaoye" × "Qinzhong No.1"). The integrated map was 1913.29 cM long, covering 94.10% of the estimated genome and with an average marker density of 2.20 cM. A total of 869 markers were mapped into 19 major independent linkage groups. Growth-related traits measured over 10 consecutive years showed a significant correlation, and 89 hypothetical QTLs were forecasted and divided into 27 distinct loci. Three traits for tree height, ground diameter and crown diameter detected 25 QTLs (13 loci), 32 QTLs (17 loci) and 15 QTLs (10 loci), respectively. Based on the BLASTX results, six candidate genes were obtained. It is important to explore the growth-related genetic mechanism and lay the foundation for the genetic improvement of E. ulmoides at the molecular level.

Conflicts of Interest:
The authors declare no conflict of interest.