Genome-Wide Association Study With Growth-Related Traits and Secondary Metabolite Contents in Red- and White-Heart Chinese Fir

Chinese fir [Cunninghamia lanceolata (Lamb.) Hook] is an important evergreen coniferous tree species that is widely distributed in many southern provinces of China and has important economic value. The Chinese fir accounts for 1/4 and 1/3 of the total artificial forest area and stock volume, respectively. Red-heart Chinese fir is popular in the market because of its high density and red heartwood. The long-growth cycle hindered the breeding process of Chinese fir, while molecular marker-assisted breeding could accelerate it. However, Chinese fir, a perennial conifer species, has a large genome, which has not yet been published. In this study, the growth-related traits and secondary metabolite contents of red- and white-heart Chinese fir were measured and found to be different between them. There are extremely significant differences among growth-related traits (p < 0.001), but secondary metabolite contents have different correlations due to differences in chemical structure. Moreover, genotype effect analysis of the substantially correlated single nucleotide polymorphisms (SNPs) revealed that most of the loci related to each growth-related traits were different from each other, indicating a type specificity of the genes regulated different growth-related traits. Furthermore, among the loci related to secondary metabolite contents, nine loci associated with multiple metabolite phenotypes such as Marker21022_4, Marker21022_172, Marker24559_31, Marker27425_37, Marker20748_85, Marker18841_115, Marker18841_198, Marker65846_146, and Marker21486_163, suggesting the presence of pleiotropic genes. This study identified the potential SNP markers associated with secondary metabolites in Chinese fir, thus setting the basis for molecular marker-assisted selection.

Introduction

Chinese fir is a fast-growing, economical timber species widely distributed in southern China. It is well-known for its high quality and large quantity. In 2021, the artificial forest area and stock volume of Chinese fir both ranked first in China. Since 1976, breeders have been striving to breed Chinese fir varieties of high quality and quantity (Zheng et al., 2013). Among them, red-heart (RH) Chinese fir is very popular because of its red heartwood color, high density, and has higher content of cedarwood and sclareol when compared with the sapwood (Duan et al., 2015; Yang et al., 2016). In our earlier study, it was found that flavonoid metabolites were the main differential metabolites between red- and white-heart Chinese fir heartwood, especially luteolin-7-O-glucuronide-5-O-rhamnoside and luteolin-7-O-glucoside (Cao et al., 2022). Besides, in the previous study, we collected 700 Chinese fir clones, which preserved in the Longshan state forest farm of Guangdong Province. We assessed their genetic diversity, and moreover, constructed a core group of 300 representative clones by using 21 simple sequence repeat markers and wood-property traits with the aim of maximizing representation of the 700 population (Duan et al., 2017).

Metabolites refer to substances produced or consumed by metabolic process, generally including only small molecules. Metabolites are well known for a variety of functions, including defense, structure, fuel, signaling, and interactions with other organisms and themselves (as a cofactor to an enzyme; Fang et al., 2019). Metabolites are divided into primary and secondary metabolites. Unlike primary metabolites, directly involved in normal growth, development, and reproduction, the secondary metabolites often play important roles in plant coloration, resistance, and defense against microbial and insect (Celedon and Bohlmann, 2018). The core primary metabolic pathways are similar in most plant and non-plant species; however, much more secondary metabolites is specifically produced in plants (Alseekh and Fernie, 2018). The total number of plant metabolites is expected to exceed 200,000 right now, most of which are secondary metabolites or their derivatives (Schwab, 2003; Pichersky and Lewinsohn, 2011). Secondary metabolites in plants play a critical role in their interactions with the environment.

Among the secondary metabolites, polyphenolic metabolites are the most abundant and widely distributed in nature (Tohge et al., 2013). Flavonoids, as the most abundant kind of polyphenols, were estimated to contain more than 8,000 compounds, which play important roles in many biological processes, such as resistance to biotic/abiotic stress (Winkelshirley, 2002) and some chronic diseases, including cardiovascular diseases and certain cancers (Luo et al., 2008). Flavonoids are composed of a common C6-C3-C6 structure with two benzene rings (ring A and ring B) interconnected by a three-carbon heterocyclic pyran ring (ring C), which can be subdivided into six categories according to their chemical structures, namely flavones, flavonols, flavanones, flavanols, isoflavones, and anthocyanins (Figure 1), with a total of more than 5,000 species (Peng et al., 2017). Flavonoids and their derivatives are primarily synthesized by a variety of enzymes involved in phenylalanine and flavonoid biosynthesis pathways. Phenylalanine synthesized in the shikimate pathway is cleaved by phenylalanine ammonia lyase to generate trans-cinnamic acid, which can be used to synthesize lignin, lignans, and flavonoids (Fraser and Chapple, 2011; Maeda and Dudareva, 2012). Subsequently, a cytochrome 450 monooxygenase, catalyzed by cinnamic acid 4-hydroxylase hydroxylates the C4 position of cinnamic acid to produce 4-coumaric acid. Naringenin chalcone is synthesized by 4-coumarate CoA ligase and catalyzed by ATP consumption. It is used as a substrate to synthesize flavonoids and their derivatives by chalcone synthase, chalcone isomerase, flavonoid 3-hydroxylase, flavonoid synthase, flavonoid 3′-hydroxylase, flavonoid 3′5’-hydroxylase, flavonoid synthase I, and UDP-glycosyltransferase.

FIGURE 1

Figure 1. Chemical structures and classification of flavonoids. The skeleton structure consists of A, B, and C rings with groups positions numbered. Flavonoids are sub-divided according to the hydroxylation pattern and conjugation between the aromatic rings (indicated in purple), and the chemical structures of six simple flavonoids are showed.

Plants in nature contain only certain classes of metabolites, so most plant metabolites show species-specific accumulation (Fernie et al., 2004). The diversity of metabolites within plants makes them ideal models for dissecting the genetic basis of metabolite accumulation. Metabolic pathways have been dissected in several ways, including screening mutant libraries (Yonekura-Sakakibara et al., 2014), analyzing gene families (Sadre et al., 2016), and conducting comparative genomics studies (Huang et al., 2017). To gain a deeper understanding of the genetic and biochemical basis of diversity in the metabolome and to artificially synthesis natural metabolites, next-generation sequencing-based linkage and association maps have been applied to metabolomic studies (Luo, 2015; Alseekh and Fernie, 2018). Not only that the genome-wide association studies (GWASs) has also been applied to the study of plant metabolome diversity due to the development of genome sequencing technologies (Fang and Luo, 2019). Metabolite-based GWAS have the advantage of greatly improving large-scale gene-metabolite interactions, annotation and identification in plants, and provide insights into the basic genetic studies of plant metabolome (Chen et al., 2014). Metabolite-based GWAS has now become one of the most using tools for identifying genetic determinants of plant metabolic diversity (Fang and Luo, 2019).

Metabolites are considered to be the bridge between the genome and phenotype and, in some cases, can be the cause or marker of morphological features (Saito, 2013). As a result, increasing interest has been focused on studying the association of metabolites, an intermediate trait associated with the biochemical and physiological states of plants (Chen et al., 2014; Wen et al., 2014). Plant height and grain shape, as complex traits, are generally determined by many loci with less impact (Huang et al., 2010), while metabolite content, especially secondary metabolite content, is usually regulated by a few more influential loci (Chen et al., 2014). Plant primary and secondary metabolites have been extensively studied, and hundreds of quantitative trait loci have been detected in Arabidopsis, rice, and tomato (Keurentjes et al., 2006; Schauer et al., 2006; Chan et al., 2010; Matsuda et al., 2012). Association analysis of core oil biosynthesis with single nucleotide polymorphism (SNP) loci found that 74 of 1 million loci in 368 core clones were significantly associated with grain oil concentration and fatty acid composition (p < 1.8E-6; Li et al., 2013). Molecular marker assisted breeding selection is based on that the strong association of the molecular markers with the target trait genes. By identifying molecular markers, the genes responsible for the target trait can also be detected.

In this study, we used an existing core population of 300 Chinese fir clones, including 48 RH and 252 white-heart (WH) Chinese fir, to examine both phenotypes, such as growth-related traits and secondary metabolite content in heartwood, and genetic polymorphisms to identify correlations between RH and WH Chinese fir. The identification of genetic polymorphisms and phenotypes correlation loci provides a theoretical basis for future research on RH Chinese fir metabolites and lays a foundation for molecular marker assisted breeding of RH Chinese fir.

Materials and Methods

Plant Materials

Three hundred Chinese fir clones were collected from the ex-situ gene bank of the Longshan State Forest Farm, Guangdong Province, China (25°11′N, 113°28′E, 285–296 m above sea level). Six growth-related traits and six secondary metabolite contents in the heartwood were measured in all 300 clones, with at least three randomly selected ramets per clone. Growth-related traits included tree height (H), diameter at breast height (DBH), stem volume (V), percentage of heartwood (P), wood basic density (WBD), and hygroscopicity (HY). The secondary metabolites included dihydrokaempferol (DK), dihydroquercetin (DQ), luteolin, pinocembrin, apigenin, and naringenin. Three individuals were selected as biological replicates.

Growth-Related Traits and Secondary Metabolites Detection

H and DBH were measured during the field surveys and V, P, WBD, and HY were measured in the laboratory using methods described by Duan et al. (2016). V = 0.000058777042 × DBH^1.9699831 × H^0.89646157. P was measure by using the formula P = r²/R² × 100%, which r and R represent the length of heartwood and xylem, respectively. WBD and HY were measured using formula WBD = 1/(W1/W2 − 0.346) and HY = (W1 − W2)/W2, which W1 and W2 represent the water-saturated weight and the oven-dry weight, respectively. The contents of six secondary metabolites, including DK, DQ, luteolin, pinocembrin, apigenin, and naringenin were measured by using liquid chromatography with mass spectrometry (LC-MS/MS), and the detailed methods and parameters of LC-MS/MS were described by Cao et al. (2022). V_G (genetic variance), V_E (environment variance), and residual were measured by using SPSS 24.0, and H_B² (generalized heritability) was calculated by using formula below:

$H_B^2=\frac{V_G}{V_G+V_E+\mathrm{residual}}\times 100\%$

DNA Extraction and Library Construction

Three-hundred microgram of Chinese fir needles were used for DNA extraction. The DNA quality and concentration of the individuals in the group were detected, the DNA quality was detected by 1% agarose gel, and the electrophoresis results were observed and recorded by the gel imaging system, and stored in an ultra-low temperature refrigerator.

At present, no genome sequence information of Chinese fir has been released, and there are no genomes of other closely related species for reference, while the conifer species that have released genome sequences are white spruce, Norway spruce, loblolly pine and Chinese pine (Birol et al., 2013; Nystedt et al., 2013; Neale et al., 2014; Niu et al., 2022), therefore, based on the genome size and GC content of Chinese fir, the Norway spruce genome was finally selected as the reference genome for electronic digestion prediction. The Norway spruce (Picea abies) genome was used as a reference for electronic digestion. The assembled genome size was 12Gb, and the GC content was 37.87%. The download address is: http://congenie.org/.

SNPs Marker Development and Locus-Phenotype Association Study

Using the deepest sequence in each SLAF tag as the reference sequence, the reads obtained by sequencing were aligned to the reference sequence, and SNPs were developed using two software, GATK (McKenna et al., 2010) and SAMtools (Li et al., 2009). The intersection of SNP markers obtained by the two methods was used as the final reliable SNP marker dataset, according to the integrity ≥0.5, the minor genotype frequency (MAF) ≥ 0.05, the obtained SNPs were filtered, and finally the population SNPs with high quality and good consistency were obtained for further analysis.

Based on high-quality SNPs, association analysis was performed on growth-related traits and flavonoid contents traits of associated populations using the mixed linear model (MLM) of TASSEL 5.0 software (Bradbury et al., 2007). The formula is calculated as follows:

$Y=u+X\alpha +Q\beta +K\mu +e$

MLM uses Q + K, which is a comprehensive analysis method of population structure and kinship, and substitutes the kinship coefficient matrix into the model as a random effect. The influence of other factors on the association results. The genotype vector is used as X, the phenotype value vector is used as Y, the fixed effect value vector is used as u, α is the weight vector for each marker, β is the trait vector for each subgroup, μ is the association matrix, and e is the residual. MLM correction will be more stringent, so it also has a relatively high accuracy, reducing the occurrence of false positives, but some useful SNP sites may be lost due to too high standards.

Results

Variation in Growth-Related Traits and Secondary Metabolite Content in Heartwood

In this study 300 Chinese fir clones, including 48 RH and 252 WH Chinese fir, were analyzed. First, H, DBH, V, P, WBD, and HY were measured (Table 1). In RH, the H, DBH, and HY were 7.24 m, 12.64 cm, and 237.65%, respectively, which resulted lower than those obtained from the WH (8.08 m, 13.72 m, and 263.92%). In contrast, the V, P, and WBD in RH were 0.0620 m³, 22.26%, and 0.3393 g.cm⁻³, respectively, which resulted higher than those measured from the WH. Moreover, the amplitude of all six growth-related traits in RH were 4.50–11.67 (H), 5.77–21.73 (DBH), 0.0093–0.2273 (V), 13.43–40.61 (P), 0.2512–0.4783 (WBD), and 152.98–333.62 (HY), with the ratios of the maximum to the minimum being lower than those calculated for the WH, thus demonstrating that RH Chinese fir grows more evenly.

TABLE 1

Table 1. Differences in growth-related traits of red- and white-heart Chinese fir.

Furthermore, the content of the following metabolites was examined in both the RH and WH: DK, DQ, luteolin, pinocembrin, apigenin, and naringenin (Table 2). DK, DQ and naringenin in RH were 0.1831, 0.1425, and 3.7770 μg.g⁻¹, respectively, which resulted higher than contents of the three metabolites in WH (0.1021, 0.1385, and 2.5778 μg.g⁻¹). However, the contents of luteolin, pinocembrin, and apigenin in RH were 0.0485, 0.7882, and 0.0968 μg.g⁻¹, respectively, which were lower than those in WH (0.0523, 1.6461, and 0.2124 μg.g⁻¹). Among them, the contents of pinocembrin and apigenin in WH were more than 2-fold higher when compared with their content in RH. The amplitude of all six secondary metabolites showed a similar pattern to growth-related traits, indicating that the variation range of RH was lower than that of WH.

TABLE 2

Table 2. Differences in the content of secondary metabolites in the heartwood of red- and white-heart Chinese fir.

The variance components of each trait in the associated population are analyzed and showed in Supplementary Table S1. It was found that V_E of each trait was much lower than V_G. The generalized heritability values were greater than 40% except P, with only 15.15%.

Correlations Between Pairs of Traits

The overall correlation analysis of the selected 12 associated phenotypes is shown in Figure 2. A significant positive correlation among H, DBH, V, and HY was detected, whereas H, DBH, V, and HY were significantly negatively correlated with WBD (p < 0.001). Furthermore, no significant difference between WH and RH was observed. We detected high correlation coefficients (r > 0.80) among H, DBH, and V. Regardless of growth-related traits or secondary metabolite contents, there was no significant correlation between P and the other 11 traits, and the correlation coefficient (r) was only approximately 0.1. A significant negative correlation between WBD and the other growth-related traits, except P, was detected, indicating that the faster the growth of Chinese fir, the lower the density of wood, and the correlation coefficient with HY was the largest (r > 0.96). In the correlation analysis of H, DBH, and V, luteolin, pinocembrin, and apigenin showed the same trends. Although the contents of the three flavonoid metabolites were positively correlated with the three growth traits (p < 0.001), the correlation coefficient was low (r ≈ 0.25), indicating the lack of a linear correlation.

FIGURE 2

Figure 2. Phenotypic differences and correlations of six growth-related traits and the contents of six secondary metabolites in red-heart and white-heart Chinese fir. H, tree height; DBH, diameter at breast height; V, volume; P, percentage of heartwood; WBD, wood basic density; HY, water absorption; DK, dihydrokaempferol; DQ, dihydroquercetin; Pino: pinocembrin; and Na, Naringenin.

In addition, there were some correlation differences in the secondary metabolite contents (Figure 2). In general, the correlation of secondary metabolites was divided into two categories: an extremely significant positive correlation among DK, DQ, and naringenin, and an extremely significant positive correlation among luteolin, pinocembrin, and apigenin. Specifically, the contents of luteolin and the other five metabolites were significantly different, among which, there was a linear correlation between luteolin and DK (r = 0.76), DQ (r = 0.72), and naringenin (r = 0.70) only in RH. The higher correlation coefficients observed in RH as compared with those in WH suggest that a linear correlation between the compounds exists only in RH. Moreover, DQ was only significantly positively correlated with DK, luteolin, and naringenin, and the correlation coefficient in RH was significantly higher than that in WH, which also showed a linear correlation between these compounds in RH. In addition, there were extremely significant positive correlations between apigenin and pinocembrin only with luteolin and no significant correlations with other metabolites. There were different degrees of correlation between naringenin and the other five metabolites, but the correlation coefficient between naringenin and other metabolites in RH was significantly higher than that in WH, indicating that naringenin is more linearly related to other secondary metabolites in RH.

GWAS of Traits of Chinese Fir

In our earlier study, simplified genome sequencing of 300 Chinese fir core populations was performed using the SLAF-seq method, and 955,503 SNPs were detected. The obtained population SNPs were further filtered according to integrity ≥0.5 and MAF ≥ 0.05, and finally 166,646 high-quality SNP markers were developed (Duan et al., 2017). The SNP information statistics of each sample were shown in Supplementary Table S2. Six growth-related traits and six flavonoid secondary metabolite contents were measured, and then TASSEL5.0 software mixed linear models (MLM) were used to detect associations between SNP marker loci and traits.

Through the qq plot, it was found that the x-axis was a uniformly distributed value (theoretical value) and had undergone −log10 transformation, and the y-axis was the actual value of p (observed value) and had undergone −log10 transformation. It was caused by random drift, and then begins to deviate from the straight line due to the existence of random drift and the real site, indicating that the association analysis based on the MLM model was suitable for this study (Figures 3, 4). The Manhattan plot showed the inability of the simplified genomic analysis to assign SNP loci to a certain chromosome. Furthermore, in the Manhattan plot, the higher the point, the more significant the association between the marker loci and traits.

FIGURE 3

Figure 3. The qq-plot and manhattan of genome-wide association studies for growth-related traits based on mixed linear model (MLM). (A) H, (B) DBH, (C) V, (D) P, (E) WBD, and (F) HY.

FIGURE 4

Figure 4. The qq-plot and manhattan of genome-wide association studies for secondary metabolite contents based on MLM. (A) Luteolin, (B) pinocembrin, (C) apigenin, (D) naringenin, (E) DK, and (F) DQ.

Under the significance criterion of p < 1E-4, 62 SNP marker loci were detected that were significantly associated with each growth-related trait using a MLM. Under MLM, 9, 4, 32, 5, 9, and 3 SNP markers were significantly associated with H, DBH, V, P, WBD, and HY, respectively (Table 3). On average, there were 10 SNP markers that were significantly associated with each growth-related traits, with the most were associated with V. Statistical analysis of the association interpretation rate of growth-related traits revealed that the overall rate varied from 6.480% (Marker19980_201_SNP associated with DBH) to 12.141% (Marker49894_207_SNP associated with V), with an average interpretation rate of 8.023%. Among them, the SNP markers associated with H accounted for 7.207%–8.659% of phenotypic variation, the SNP markers associated with DBH accounted for 6.480%–9.244% of phenotypic variation, and the SNP markers associated with V accounted for 6.491%–12.141%, the SNP markers associated with P was 6.750%–10.239%, the SNP markers associated with WBD was 6.724%–9.169%, and the SNP markers associated with HY ranged from 6.899% to 7.714%.

TABLE 3

Table 3. SNPs with significant association to growth-related traits based on MLM genome-wide association study (GWAS).

Under the significant criterion of p < 1E-7, using MLM, a total of 163 SNP marker locus were detected that significantly associated with the content of various flavonoid metabolites. Under MLM, 15, 5, 40, 32, 37, and 34 SNP markers were significantly associated with DK, naringenin, DQ, luteolin, pinocembrin, and apigenin, respectively (Table 4), On average, there were 27 significantly associated SNP markers for each growth-related trait, with the most SNP markers associated with DQ. Statistical analysis of the associated interpretation rate of metabolite content revealed that the overall interpretation rate varied from 11.513% (Marker21373_159_SNP associated with taxifolin) to 44.750% (Marker21022_4_SNP associated with apigenin), with an average interpretation rate of 17.833%, higher than the results of association analysis related to growth-related traits. Among these markers associated with Chinese fir growth-related traits, 41 markers were associated with multiple traits, of which 27 were associated with the content of two metabolites and 14 were associated with the content of three metabolites, which was significantly higher than the results of the association analysis of growth-related traits, indicating that the genes regulating flavonoid metabolites were relatively concentrated. Correlation analysis between phenotypic traits also showed that there was a certain degree of correlation between the tested traits, suggesting the presence of pleiotropic genes, and the phenotypic correlation was verified at the molecular level. These SNPs jointly control different phenotypic traits and can be used as key molecular markers in future resource identification, tree selection, and other studies.

TABLE 4

Table 4. SNPs with significant association to flavonoid contents based on MLM GWAS.

Genotype Effects of Most Significantly Associated SNP With Each Traits

To further analyze the genotype effect of key loci, we selected the five SNP loci with the smallest p value, representing the five SNP loci with the most significant associations, and assigned them to perform genotype effect analysis (Figures 3, 4). By analyzing the genotype effect of the five associated SNP loci, we found that the homozygous and heterozygous genotypes showed different degrees of significant differences in association with growth-related traits and secondary metabolite contents.

The five SNP loci that were most significantly associated with each trait were selected and analyzed for genotype effects (Figures 5, 6). Overall, we identified differences among the SNP marker loci significantly associated with each growth-related trait, with only two loci found to be associated with multiple traits, including Marker46641_94 in H, DBH, and V, and Marker21403_124 in WBD and HY. For Marker46641_94, accessions with the homozygous TT genotype were associated with higher H, DBH, and V values, indicating a possible association with the growth of Chinese fir. For Marker21403_124, accessions with the heterozygous CT genotype were associated with a lower WBD and higher HY, indicating a possible association with wood quality. In association study with secondary metabolite contents, there were more single loci correlated with multiple traits, including Marker21022_4 and Marker21022_172 in luteolin, pinocembrin, and apigenin; Marker24559_31 and Marker27425_37 in luteolin and pinocembrin; Marker20748_85 in pinocembrin and apigenin; Marker18841_115, Marker18841_198, and Marker21486_163 in naringenin and DK; and Marker65846_146 in naringenin, DK, and DQ. For Marker21022_4 and Marker21022_172, accessions with homozygous CC and AA, respectively, were associated with higher luteolin, pinocembrin, and apigenin contents, and accessions with homozygous CC, TT, and AA in Marker24559_31, Marker27425_37, and Marker20748_85 showed the same trend. Marker18841_115, Marker18841_198, and Marker65846_146 accessions with heterozygous TC, AC, and CT, respectively, were associated with a higher content of naringenin and DK. In addition, accessions with heterozygous CT in Marker65846_146 were associated with higher DQ content. Moreover, for Marker21486_163, accessions with homozygous AA were associated with a higher content of DQ.

FIGURE 5

Figure 5. Box plots showing genotype effects of the top five single nucleotide polymorphisms (SNPs) associated with growth-related traits. (A) H, (B) DBH, (C) V, (D) P, (E) WBD, and (F) HY.

FIGURE 6

Figure 6. Box plots showing genotype effects of the top five SNPs associated with secondary metabolite contents. (A) Luteolin, (B) pinocembrin, (C) apigenin, (D) naringenin, (E) DK, and (F) DQ.

Discussion

Growth-related traits are important phenotypic traits of forest trees and are closely related to wood quality. Breeders have been working to improve wood quality by enhancing these traits. Secondary metabolites play a crucial role in heartwood, not only affecting the corrosion resistance and insect resistance of most woody plant heartwood, but also producing color substances through chemical reactions that bring coloration to the heartwood, further improving its marketability. Through the difference analysis of the 12 phenotypic traits in 300 Chinese fir clones, we reported that RH grew more slowly but the wood quality was higher when compared with that of WH Chinese fir. Moreover, analyzing the ratios of maximum and minimum measurements of the growth-related traits, we found that the phenotype of RH changes less, the characteristics are more stable, and the growth is more uniform. The determination of the correlation of the six traits with the content of secondary metabolites revealed that the contents of DK, DQ, and naringenin were higher in RH. DK, DQ, and naringenin belong to the middle and upper metabolites in the flavonoid biosynthesis pathway, and a higher content indicates more sufficient substrates for the synthesis of downstream products. In the middle and lower reaches of the pathway, the lower content of luteolin, pinocembrin, and apigenin in RH may be due to glycosylation and other processes, resulting in more stable glycosylation products, indicating that the glycosylation process in RH may also be more active than that in WH. Compared with WH, RH has more upstream substrates and downstream glycosylation products, which ultimately changes the color of heartwood. In a previous study, we found that the differentially expressed genes in phenylalanine pathway (PAL, C4H, and 4CL), upstream of flavonoid biosynthesis pathway (CHS, CHI and F3H) and glycosylation (UGT) were all up-regulated in the transition zone of RH compare to WH (Cao et al., 2022).

By analyzing the 12 phenotypic traits of the 300 core populations screened, it was shown that the traits of this population were little affected by the environment and can be inherited stably (Supplementary Table S1). H and DBH represent the growth rate of woody plants and V represents the yield of wood, which can clarify the growth advantage of the clones. WBD is another important indicator of wood growth-related traits that can affect the final wood quality. Wood density varies widely among different individuals and generally has strong heritability, which is weakly affected by environmental factors (Louzada and Fonseca, 2002). Percentage of heartwood (P), which represents the yield of heartwood, is an important indicator for evaluating the growth of Chinese fir. In the present study, no significant correlations were found with other growth-related traits. In previous studies, it was demonstrated that a large amount of metabolic activity occurred in the transition zone of xylem, with the death of ray parenchyma cells and the disappearance of starch grains, the transition zone would transform into heartwood (Nakaba et al., 2012; Cao et al., 2022). HY refers to the weight increase of the wood after it is completely soaked in water, which is affected by WBD. In this study, a significant negative correlation with WBD was detected, which was consistent with the results of other studies (Duan et al., 2016). We found that most growth-related traits had significant positive correlations between pairs, with the exception of WBD, which had a significant negative correlation compared with H, DBH, V, and HY, and was consistent with previous literature (Duan et al., 2016). In addition, there was no significant correlation between P and other phenotypes in this study, and similar trends have rarely been reported in other species.

GWAS are observational studies of a set of genome-wide genetic variants in the populations and their association with specific traits. Recent years, GWAS have been shown to be very powerful in analyzing the genetic basis of complex phenotypic variation, primarily applied to the study of physiological and agricultural traits in plants (Huang et al., 2010; Tian et al., 2011; Horton et al., 2012; Chen et al., 2016; Yano et al., 2016; Xiao et al., 2017). In a MLM, the population structure and kinship are comprehensively analyzed, making it the most common used GWAS method. Besides, it also has demonstrated high efficiency in correcting for inflation caused by many small genetic effects, and for controlling for bias in population stratification (Yu et al., 2006; Zhao et al., 2007). In our investigation, to reduce the occurrence of false-positive results in association analysis, the MLM model in TASSEL5.0, was used. The qq plot showed that MLM was suitable for the association analysis of Chinese fir growth-related traits and secondary metabolite contents; however, if the same threshold was used (p < 1E-4), the number of association sites for secondary metabolite contents was significantly higher than that for growth-related traits. Therefore, we set the significant association threshold for secondary metabolites to p < 1E-7 and obtained SNP sites with more significant associations. Using MLM, we obtained 62 SNP marker loci related to growth-related traits and 163 SNP marker loci related to secondary metabolite contents. Since there may be corresponding genes near the obtained association loci (Morris et al., 2013), the genes containing these SNPs loci may be involved in the growth of fir or the synthesis of flavonoid secondary metabolites. Palle et al. (2013) used the SNPs markers to analyze the association between the expression levels of genes potentially related to loblolly pine wood formation, thus contributing to a deeper understanding of phenotypic variation in Chinese fir. In addition, with the rapid development of sequencing technology, an increasing number of tree species, including Chinese fir, will undergo whole-genome sequencing, which will help us to further explain the associated loci obtained in this study.

The GWAS performed using MLM had additive or dominant effects or both. This indicated that the phenotypic of heterozygous Chinese fir may be better than those of homozygous samples, which was also proved by the SNPs-phenotype correlation analysis. For example, Marker42910_4 was related to H, Marker36988_54 was related to DBH, Marker49894_207 was related to V, Marker31544_181 was related to WBD, and Marker21169_24 was related to HY. It was found that for heterozygous individuals at these loci in the related population, the numerical mean of heterozygous was significantly higher than that of homozygous individuals, and these loci with dominant or additive effects may be the important genetic basis for heterosis in Chinese fir. We also observed that the same marker was associated with multiple traits, primarily secondary metabolites. For example, Marker21022_4 and Marker21022_172 co-existed in association with luteolin, pinocembrin, and apigenin, whereas Marker18841_115, Marker18841_198, Marker65846_146, and Marker21486_163 were also significantly associated with naringenin and DK, and most of the significant association loci in DQ were special. Therefore, we also divided the detected flavonoid secondary metabolites into three categories according to the level of significant association, which indicated that the flavonoid metabolites may have different regulatory genes owing to their different chemical structures. In recent years, there have been many reports of GWAS associated with metabolism. Xiao et al. (2021) identified eight quantitative trait loci significantly associated with the metabolic pathway of Salicylic acid (SA) biosynthesis in Populus. During GWAS of tea, it was found that 307 SNP marker were proved to be related to some flavor-related metabolites such as theanine, caffeine and catechins, and some of the markers were even pleiotropic (Fang et al., 2021). In the present study, some single markers are associated with multiple traits and illustrate the close correlation between secondary metabolite traits, indicating pleiotropic effects, and reflecting the complexity of physiological and biochemical metabolic processes. This phenomenon has been reported in several studies. For example, Tian et al. found that the SNP10 locus located in the PtGA20Ox gene was significantly associated with one wood growth trait and three texture traits, simultaneously (Tian et al., 2012). Yang et al. (2015) also found that the SNP11 locus located in the Pto-Wuschela gene was significantly associated with tree height and volume per plant. The sites associated with multiple traits may be more effective in molecular marker-assisted selection and may also be helpful for mining important genomic regions.

Conclusion

This study found that there were differences in growth between white- and red-heart Chinese fir. Red-heart Chinese fir grows more slowly but the wood quality is higher. The 12 traits of this population were little affected by the environment and can be inherited stably. Furthermore, several single SNP loci which significantly associated with secondary metabolite contents has been found, indicating pleiotropic effects. Chinese fir is a perennial conifer with a huge genome. As of this manuscript writing, its genome has not been published. At present, it is impossible to locate SNP loci in chromosomes and genes by simplified GWAS. In the future research, the loci can be accurately located with reference to the genome to promote the process of molecular marker-assisted breeding of Chinese fir.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: We have uploaded the data on GSA data base at https://ngdc.cncb.ac.cn/, and the GSA accession number is CRA007096.

Author Contributions

YL and HZ conceived and designed the experiments. SC wrote the article. SC, HD, and BW performed the experiments. SC, YS, and RH analyzed the data. JL and WD participated in and helped to complete the experiments. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the Key-Area Research and Development Program of Guangdong Province (no. 2020B020215001), the Science and Technology Research Project of Beijing Forestry University (2018WS01), and the National Natural Science Foundation of China (no. 31972956).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.922007/full#supplementary-material

Supplementary Table S1 | The variance components and the generalized heritability of various characters in associated group.

Supplementary Table S2 | The SNP information.

References

Alseekh, S., and Fernie, A. R. (2018). Metabolomics 20 years on: what have we learned and what hurdles remain? Plant J. 94, 933–942. doi: 10.1111/tpj.13950

PubMed Abstract | CrossRef Full Text | Google Scholar

Birol, I., Raymond, A., Jackman, S. D., Pleasance, S., Coope, R., Taylor, G. A., et al. (2013). Assembling the 20 Gb white spruce (Picea glauca) genome from whole-genome shotgun sequencing data. Bioinformatics 29, 1492–1497. doi: 10.1093/bioinformatics/btt178

PubMed Abstract | CrossRef Full Text | Google Scholar

Bradbury, P. J., Zhang, Z., Kroon, D. E., Casstevens, T. M., Ramdoss, Y., and Buckler, E. S. (2007). TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics 23, 2633–2635. doi: 10.1093/bioinformatics/btm308

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, S., Deng, H., Zhao, Y., Zhang, Z., Tian, Y., Sun, Y., et al. (2022). Metabolite profiling and transcriptome analysis unveil the mechanisms of red-heart Chinese fir (Cunninghamia lanceolata (lamb.) hook) heartwood coloration. Front. Plant Sci. 13:594. doi: 10.3389/fpls.2022.854716

CrossRef Full Text | Google Scholar

Celedon, J. M., and Bohlmann, J. (2018). An extended model of heartwood secondary metabolism informed by functional genomics. Tree Physiol. 38, 311–319. doi: 10.1093/treephys/tpx070

PubMed Abstract | CrossRef Full Text | Google Scholar

Chan, E. K., Rowe, H. C., Hansen, B. G., and Kliebenstein, D. J. (2010). The complex genetic architecture of the metabolome. PLoS Genet. 6:e1001198. doi: 10.1371/journal.pgen.1001198

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., Gao, Y., Xie, W., Gong, L., Lu, K., Wang, W., et al. (2014). Genome-wide association analyses provide genetic and biochemical insights into natural variation in rice metabolism. Nat. Genet. 46, 714–721. doi: 10.1038/ng.3007

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, W., Wang, W., Peng, M., Gong, L., Gao, Y., Wan, J., et al. (2016). Comparative and parallel genome-wide association studies for metabolic and agronomic traits in cereals. Nat. Commun. 7, 1–10. doi: 10.1038/ncomms12767

CrossRef Full Text | Google Scholar

Duan, H., Cao, S., Zheng, H., Hu, D., Lin, J., Cui, B., et al. (2017). Genetic characterization of Chinese fir from six provinces in southern China and construction of a core collection. Sci. Rep. 7, 1–10. doi: 10.1038/s41598-017-13219-0

CrossRef Full Text | Google Scholar

Duan, H., Cao, S., Zheng, H., Hu, D., Lin, J., Lin, H., et al. (2016). Variation in the growth traits and wood properties of Chinese fir from six provinces of southern China. Forests 7:192. doi: 10.3390/f7090192

CrossRef Full Text | Google Scholar

Duan, H., Hu, R., Wu, B., Chen, D., Huang, K., Dai, J., et al. (2015). Genetic characterization of red-colored heartwood genotypes of Chinese fir using simple sequence repeat (SSR) markers. Genet. Mol. Res. 14, 18552–18561. doi: 10.4238/2015.December.28.2

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, C., Fernie, A. R., and Luo, J. (2019). Exploring the diversity of plant metabolism. Trends Plant Sci. 24, 83–98. doi: 10.1016/j.tplants.2018.09.006

CrossRef Full Text | Google Scholar

Fang, C., and Luo, J. (2019). Metabolic GWAS-based dissection of genetic bases underlying the diversity of plant metabolism. Plant J. 97, 91–100. doi: 10.1111/tpj.14097

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, K., Xia, Z., Li, H., Jiang, X., Qin, D., Wang, Q., et al. (2021). Genome-wide association analysis identified molecular markers associated with important tea flavor-related metabolites. Hortic. Res. 8, 42. doi: 10.1038/s41438-021-00477-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernie, A. R., Trethewey, R. N., Krotzky, A. J., and Willmitzer, L. (2004). Metabolite profiling: from diagnostics to systems biology. Nat. Rev. Mol. Cell Biol. 5, 763–769. doi: 10.1038/nrm1451

PubMed Abstract | CrossRef Full Text | Google Scholar

Fraser, C. M., and Chapple, C. (2011). The phenylpropanoid pathway in Arabidopsis. Arabidopsis Book 9:e0152.

Google Scholar

Horton, M. W., Hancock, A. M., Huang, Y. S., Toomajian, C., Atwell, S., Auton, A., et al. (2012). Genome-wide patterns of genetic variation in worldwide Arabidopsis thaliana accessions from the RegMap panel. Nat. Genet. 44, 212–216. doi: 10.1038/ng.1042

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, A. C., Kautsar, S. A., Hong, Y. J., Medema, M. H., Bond, A. D., Tantillo, D. J., et al. (2017). Unearthing a sesterterpene biosynthetic repertoire in the Brassicaceae through genome mining reveals convergent evolution. Proc. Natl. Acad. Sci. 114, E6005–E6014. doi: 10.1073/pnas.1705567114

CrossRef Full Text | Google Scholar

Huang, X., Sang, T., Zhao, Q., Feng, Q., Zhao, Y., Li, C., et al. (2010). Genome-wide association studies of 14 agronomic traits in rice landraces. Nat. Genet. 42, 961–967. doi: 10.1038/ng.695

PubMed Abstract | CrossRef Full Text | Google Scholar

Keurentjes, J. J., Fu, J., De Vos, C., Lommen, A., Hall, R. D., Bino, R. J., et al. (2006). The genetics of plant metabolism. Nat. Genet. 38, 842–849. doi: 10.1038/ng1815

CrossRef Full Text | Google Scholar

Li, H., Peng, Z., Yang, X., Wang, W., Fu, J., Wang, J., et al. (2013). Genome-wide association study dissects the genetic architecture of oil biosynthesis in maize kernels. Nat. Genet. 45, 43–50. doi: 10.1038/ng.2484

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, R., Yu, C., Li, Y., Lam, T.-W., Yiu, S.-M., Kristiansen, K., et al. (2009). SOAP2: an improved ultrafast tool for short read alignment. Bioinformatics 25, 1966–1967. doi: 10.1093/bioinformatics/btp336

PubMed Abstract | CrossRef Full Text | Google Scholar

Louzada, J. L. P., and Fonseca, F. M. (2002). The heritability of wood density components in Pinus pinaster Ait. and the implications for tree breeding. Ann. For. Sci. 59, 867–873. doi: 10.1051/forest:2002085

CrossRef Full Text | Google Scholar

Luo, J. (2015). Metabolite-based genome-wide association studies in plants. Curr. Opin. Plant Biol. 24, 31–38. doi: 10.1016/j.pbi.2015.01.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Luo, J., Butelli, E., Hill, L., Parr, A., Niggeweg, R., Bailey, P., et al. (2008). AtMYB12 regulates caffeoyl quinic acid and flavonol synthesis in tomato: expression in fruit results in very high levels of both types of polyphenol. Plant J. 56, 316–326. doi: 10.1111/j.1365-313X.2008.03597.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Maeda, H., and Dudareva, N. (2012). The shikimate pathway and aromatic amino acid biosynthesis in plants. Annu. Rev. Plant Biol. 63, 73–105. doi: 10.1146/annurev-arplant-042811-105439

CrossRef Full Text | Google Scholar

Matsuda, F., Okazaki, Y., Oikawa, A., Kusano, M., Nakabayashi, R., Kikuchi, J., et al. (2012). Dissection of genotype–phenotype associations in rice grains using metabolome quantitative trait loci analysis. Plant J. 70, 624–636. doi: 10.1111/j.1365-313X.2012.04903.x

PubMed Abstract | CrossRef Full Text | Google Scholar

McKenna, A., Hanna, M., Banks, E., Sivachenko, A., Cibulskis, K., Kernytsky, A., et al. (2010). The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303. doi: 10.1101/gr.107524.110

PubMed Abstract | CrossRef Full Text | Google Scholar

Morris, G. P., Ramu, P., Deshpande, S. P., Hash, C. T., Shah, T., Upadhyaya, H. D., et al. (2013). Population genomic and genome-wide association studies of agroclimatic traits in sorghum. Proc. Natl. Acad. Sci. 110, 453–458. doi: 10.1073/pnas.1215985110

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakaba, S., Begum, S., Yamagishi, Y., Jin, H. O., Kubo, T., and Funada, R. (2012). Differences in the timing of cell death, differentiation and function among three different types of ray parenchyma cells in the hardwood Populus sieboldii × P. grandidentata. Trees 26, 743–750. doi: 10.1007/s00468-011-0640-0

CrossRef Full Text | Google Scholar

Neale, D. B., Wegrzyn, J. L., Stevens, K. A., Zimin, A. V., Puiu, D., Crepeau, M. W., et al. (2014). Decoding the massive genome of loblolly pine using haploid DNA and novel assembly strategies. Genome Biol. 15:R59. doi: 10.1186/gb-2014-15-3-r59

CrossRef Full Text | Google Scholar

Niu, S., Li, J., Bo, W., Yang, W., Zuccolo, A., Giacomello, S., et al. (2022). The Chinese pine genome and methylome unveil key features of conifer evolution. Cell 185, 204.e214–217.e214. doi: 10.1016/j.cell.2021.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Nystedt, B., Street, N. R., Wetterbom, A., Zuccolo, A., Lin, Y.-C., Scofield, D. G., et al. (2013). The Norway spruce genome sequence and conifer genome evolution. Nature 497, 579–584. doi: 10.1038/nature12211

PubMed Abstract | CrossRef Full Text | Google Scholar

Palle, S. R., Seeve, C. M., Eckert, A. J., Wegrzyn, J. L., Neale, D. B., and Loopstra, C. A. (2013). Association of loblolly pine xylem development gene expression with single-nucleotide polymorphisms. Tree Physiol. 33, 763–774. doi: 10.1093/treephys/tpt054

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, M., Shahzad, R., Gul, A., Subthain, H., Shen, S., Lei, L., et al. (2017). Differentially evolved glucosyltransferases determine natural variation of rice flavone accumulation and UV-tolerance. Nat. Commun. 8, 1–12. doi: 10.1038/s41467-017-02168-x

CrossRef Full Text | Google Scholar

Pichersky, E., and Lewinsohn, E. (2011). Convergent evolution in plant specialized metabolism. Annu. Rev. Plant Biol. 62, 549–566. doi: 10.1146/annurev-arplant-042110-103814

CrossRef Full Text | Google Scholar

Sadre, R., Magallanes-Lundback, M., Pradhan, S., Salim, V., Mesberg, A., Jones, A. D., et al. (2016). Metabolite diversity in alkaloid biosynthesis: a multilane (diastereomer) highway for camptothecin synthesis in Camptotheca acuminata. Plant Cell 28, 1926–1944. doi: 10.1105/tpc.16.00193

PubMed Abstract | CrossRef Full Text | Google Scholar

Saito, K. (2013). Phytochemical genomics—a new trend. Curr. Opin. Plant Biol. 16, 373–380. doi: 10.1016/j.pbi.2013.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Schauer, N., Semel, Y., Roessner, U., Gur, A., Balbo, I., Carrari, F., et al. (2006). Comprehensive metabolic profiling and phenotyping of interspecific introgression lines for tomato improvement. Nat. Biotechnol. 24, 447–454. doi: 10.1038/nbt1192

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwab, W. (2003). Metabolome diversity: too few genes, too many metabolites? Phytochemistry 62, 837–849. doi: 10.1016/S0031-9422(02)00723-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, F., Bradbury, P. J., Brown, P. J., Hung, H., Sun, Q., Flint-Garcia, S., et al. (2011). Genome-wide association study of leaf architecture in the maize nested association mapping population. Nat. Genet. 43, 159–162. doi: 10.1038/ng.746

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, J., Du, Q., Chang, M., and Zhang, D. (2012). Allelic variation in PtGA20Ox associates with growth and wood properties in Populus spp. PLoS One 7:e53116. doi: 10.1371/journal.pone.0053116

PubMed Abstract | CrossRef Full Text | Google Scholar

Tohge, T., Watanabe, M., Hoefgen, R., and Fernie, A. R. (2013). The evolution of phenylpropanoid metabolism in the green lineage. Crit. Rev. Biochem. Mol. Biol. 48, 123–152. doi: 10.3109/10409238.2012.758083

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, W., Li, D., Li, X., Gao, Y., Li, W., Li, H., et al. (2014). Metabolome-based genome-wide association study of maize kernel leads to novel biochemical insights. Nat. Commun. 5, 1–10. doi: 10.1038/ncomms4438

CrossRef Full Text | Google Scholar

Winkelshirley, B. (2002). Biosynthesis of flavonoids and effects of stress. Curr. Opin. Plant Biol. 5, 218–223. doi: 10.1016/S1369-5266(02)00256-X

CrossRef Full Text | Google Scholar

Xiao, L., Du, Q., Fang, Y., Quan, M., Lu, W., Wang, D., et al. (2021). Genetic architecture of the metabolic pathway of salicylic acid biosynthesis in Populus. Tree Physiol. 41, 2198–2215. doi: 10.1093/treephys/tpab068

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, Y., Liu, H., Wu, L., Warburton, M., and Yan, J. (2017). Genome-wide association studies in maize: praise and stargaze. Mol. Plant 10, 359–374. doi: 10.1016/j.molp.2016.12.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Wei, Z., Du, Q., Chen, J., Wang, Q., Quan, M., et al. (2015). The genetic regulatory network centered on Pto-Wuschela and its targets involved in wood formation revealed by association studies. Sci. Rep. 5, 1–14. doi: 10.1038/srep16507

CrossRef Full Text | Google Scholar

Yang, H., Wen, S., Xiao, F., and Fan, G. (2016). Study on the chemical compositions of ethanol extraction from Chenshan red-heart Chinese fir. South China For. Sci. 44, 35–37. doi: 10.16259/j.cnki.36-1342/s.2016.01.010

CrossRef Full Text | Google Scholar

Yano, K., Yamamoto, E., Aya, K., Takeuchi, H., Lo, P.-C., Hu, L., et al. (2016). Genome-wide association study using whole-genome sequencing rapidly identifies new genes influencing agronomic traits in rice. Nat. Genet. 48, 927–934. doi: 10.1038/ng.3596

PubMed Abstract | CrossRef Full Text | Google Scholar

Yonekura-Sakakibara, K., Nakabayashi, R., Sugawara, S., Tohge, T., Ito, T., Koyanagi, M., et al. (2014). A flavonoid 3-O-glucoside: 2″O-glucosyltransferase responsible for terminal modification of pollen-specific flavonols in Arabidopsis thaliana. Plant J. 79, 769–782. doi: 10.1111/tpj.12580

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, J., Pressoir, G., Briggs, W. H., Vroh Bi, I., Yamasaki, M., Doebley, J. F., et al. (2006). A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat. Genet. 38, 203–208. doi: 10.1038/ng1702

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, K., Aranzana, M. J., Kim, S., Lister, C., Shindo, C., Tang, C., et al. (2007). An Arabidopsis example of association mapping in structured samples. PLoS Genet. 3:e4. doi: 10.1371/journal.pgen.0030004

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, H., Hu, D., Wang, R., Wei, R., Yan, S., Zhu, H., et al. (2013). Plus tree resource survey and genebank construction for Cunninghamia lanceolata. J. Southwest. For. U. 33, 22–26. doi: 10.3969/j.issn.2095-1914.2013.01.004

CrossRef Full Text | Google Scholar