Pedigree reconstruction and spatial analysis for genetic testing and selection in a Larix kaempferi (Lamb.) Carrière plantation

Background Larix kaempferi is one of the major timber species in Northeast Asia. Demand for the reforestation of the species is rising in South Korea due to an increase in large timber production and utilization. However, progeny trials for the species have not been explored, making it challenging to foster advanced generations of tree improvement. In the present study, genetic testing and selection for diameter growth were conducted using pedigree reconstruction and phenotypic spatial distribution analysis in a plantation of L. kaempferi. The aim of the present study was to select the superior larch individuals using the pedigree reconstruction and phenotypic spatial distribution to substitute progeny trials. The plantation of seed orchard crops was established in 1990 and one-hundred and eighty-eight trees were selected as the study material. Genetic variation was investigated first to validate its adequacy as breeding material. Genetic testing was carried out using a model considering pedigree information and spatial autoregression of the phenotypes. Results The expected heterozygosity of the mother trees and offspring were 0.672 and 0.681 presenting the corresponding level of genetic variation between two groups. The pedigree reconstruction using maternity analysis assigned one to six progenies to ninety-two candidate mothers. The accuracy of genetic testing was exceedingly increased with the animal model considering AR1 ⊗ AR1 structure compared to the animal model only. The estimated genetic variance of the former was 9.086 whereas that of the latter was 4.9E-5 for DBH. The predicted breeding values of the offspring for DBH were ranged from -5.937 cm to 5.655 cm and the estimated heritability of diameter growth was 0.344. Conclusions The genetic testing approach based on pedigree reconstruction and phenotypic spatial distribution analysis was considered a useful analytical scheme that could replace or supplement progeny trials. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03530-y.


Background
The genus Larix, a deciduous tree species in the family Pinaceae, provides timber with high economic value as well as great ecological value in forests. Larch tree improvement programs have been carried out for the major timber species and hybrids in the genus [1][2][3][4]. Larix kaempferi (Lamb.) Carrière has been one of the major timber species in Korea since its introduction in 1904. Large scale plantations of the species were cultivated until the 1970s [5]. However, the plantation areas had since been reduced due the low timber production and difficulties encountered in larch wood processing for a time [6]. Conversely, demand for larch timber and seeds has been increasing since 2000s consistent with increase in timber production and developments in wood processing technologies [6]. Today, forestry policies that could expand the larch plantation and revitalize their use in Korea are being pursued in response to increasing demand. Over the past five years, the plantation area under L. kaempferi has increased considerably, contrary to the area under other major species. It reached 273,000 ha of forest area, which is the second largest after the area under Pinus densiflora (2015) [7].
Several studies have been conducted only on partial stages of the L. kaempferi tree improvement program in Korea. Provenance tests and adaptability tests had been performed to introduce the superior provenance of larch until 1987 [6]. 145 plus trees were selected during 1959-1980 and 306. 2 ha of seed orchard has been established to produce improved seeds of the species since 1968. [6]. In addition, studies to promote flowering and selection of multi-flowering clones have been carried out to address a key constraint in the unstable larch seed production [8,9]. However, the tree improvement via progeny testing of the species was suspended due low levels of historic utilization, as well as a putatively limited gene pool caused by the introduction of only a small number of individuals to Korea.
Tree improvement is one of the essential activities in forest management, which is used to enhance the economic value of forests [10]. It aims to increase genetic gain in a gradual manner with recurrent testing and selection of advanced generations [11,12]. The selection of advanced generation by genetic testing of larch has become necessary to ameliorate economic value of the species which has been actively used in plantation forestry in Korea. To shorten the time-consuming progeny trial approach of genetic testing, the alternative methods to the approach are required to meet current demand.
Pedigree reconstruction in plantations using seed orchard crops is one of the feasible approaches that could substitute the progeny trial method. Pedigree reconstruction analysis was first applied in tree improvement activities to detect the paternity arising from polymix breeding [13], followed by Breeding without Breeding [14], quasifield trials [15], ad hoc breeding [16], and in situ genetic evaluation [17]. Pedigree reconstruction analysis is carried out to support existing tree improvement programs and to improve the accuracy of genetic parameter and breeding value estimates by substituting artificial crossing or wind pollination [18][19][20][21][22].
The accurate distinction of environmental effects is essential for proper genetic testing. Genetic testing considering spatial effects on the target traits has been attempted in previous studies in progeny trials or experimental forests [23]. Analytic models accounting for spatial structure could reduce the standard error associated with typical interfamily differences by 5-20% when compared to random complete block design experiments in genetic testing for tree height [24]. In addition, the fitness of the analysis was improved using the spatial model when compared to the incomplete block model in testing the growth properties of the Pinus pinaster [25]. The model has also been demonstrated to increase the accuracy of the log-likelihood and genetic value predictions based on restricted maximum likelihood [26]. The variation due to the spatial factors could be modelled effectively complementing the spatial autocorrelation. The supplementation of the spatial autocorrelation factor was more effective in the modeling of the spatial variation, significantly reducing the residual variation in a progeny trial of Cunninghamia lanceolata [23]. The genetic factor was efficiently explained considering the residuals caused by spatial autoregression of the target traits even in progeny trials based on a block design [26][27][28][29][30].
The objective of the present study was to perform genetic testing aided by pedigree reconstruction analysis and phenotypic spatial distribution analysis. The distinctive challenges of the present study to previous studies were the absence of a priori information on the family and the experimental design of the study site. The study was conducted in phases with two detailed subjects to analyze the genetic and environmental effects at the individual level, including 1) genetic variation analysis to determine the adequacy of the individuals as breeding material followed by pedigree reconstruction analysis, and 2) genetic testing based on pedigree information and phenotypic spatial autoregression. The overall aims of the present study were to present schematic of genetic testing in plantations that could be applied 1) to acquire basic genetic variation data and pedigree information on breeding materials of seed orchard crops and 2) to ameliorate the accuracy of genetic testing in plantations while considering phenotypic spatial distribution.
The following hypotheses were formulated: 1) the plantation of the offspring from the plus trees in the seed orchard would retain a significant amount of the genetic variation even though information about the introduced population is missing and 2) the genetic and environmental effects would be detectable from growth differences within and between the local groups in the plantation.

Genetic variation in plantation of seed orchard crops of Larix kaempferi
The number of alleles and levels of heterozygosity in the three groups of 145 plus trees (PT) and offspring (HC and HS) were estimated using seven microsatellite markers. The average number of alleles (N A ) and effective alleles (N E ) were 9.7 and 3.6, respectively ( Table 1).
The estimated genetic parameters of the mother tree group were analyzed by the National Forest Seed and Variety Center [8]. Genetic variation in the L. kaempferi was slightly lower than that in the Japanese natural population (0.717-0.762) based on the average expected heterozygosity (H E ) [31]. However, the difference was considered to be due to the low number of polymorphic markers used, as the estimates were in the 0.518-0.860 range. The inference was supported by the even lower H E in the introduced 140 plus trees and the offspring in China with 17 microsatellite markers, at 0.525 and 0.557, respectively [32]. The plantation of the seed orchard crops, with the considerable genetic variation compared to that in the group of mother trees, was considered the appropriate base population.
N A and N E were similar between PT and HC, with averages of 10.0 and 3.7, respectively ( In a comparison between the two offspring groups (HC and HS), the differences in genetic variation were not significant, although some of the estimates were slightly lower in HS. The latter had lower N A , N E , H O , H E values in except in the analysis with bcLK033 and LK4178. N A in bcLK033 was similar between the two groups, but N E was somewhat higher in HS. Both estimates of LK4178 were higher in HS than in HC, exhibiting distinctiveness from the other markers.

Pedigree reconstruction of individuals in the L. kaempferi plantation
The average polymorphic information content (PIC) was 0.638 and the non-exclusion probability of identity was as low as 4.1E-07, based on the seven microsatellite markers. Each individual was identifiable as none of them showed the coincident genotypes in the seven loci. The non-exclusion probability of the first parent was 0.055. Ninety trees corresponding to 27% of the overall trees were assigned at relaxed confidence level ( Table 2). The assignment of the other 244 trees was accomplished but did not fall within the confidence level. The proportions of assignments within the relaxed confidence level were 23% and 25% in HC and HS, respectively ( Table 2). The low confidence level in the present study was considered to be derived from the characteristics of the used markers rather than the studied materials regarding the similar assignment patterns in the two discrete plantations. The maternity analysis assigned the offspring to 92 and 79 candidate mothers in HC and HS respectively (Table  S1). The total number of the candidate mothers was 112. The number of mutually exclusive mother trees was 33 in HC and 20 in HS, displaying the diverse constitution of the mother trees. The results potentially indicated that the genetic constitution would be distinguished by the reproduction year considering the samples of the seed orchard crops were produced in different years. The number of progenies by each mother tree were 1-6 trees in HC. The mother trees with six progenies were CN05, JPR-4, and CN09, and the 42 mother trees had only the one offspring. The number of progenies by each mother tree were 1-5 trees in HS. The assignment revealed that the mother trees GB12 and GB27 each had five progenies, and 36 mother trees had one offspring in HS. The offspring by mother trees were relatively evenly dispersed in the plantation (Fig S1-S5). Similarly, majority of the mother trees had one offspring according to the results of pedigree reconstruction analyses for Abies nordmanniana and Larix occidentalis [15,17].

Genetic testing in the plantation of L. kaempferi
The prediction of breeding values was performed for each mother tree and the individual tree in the plantation although the confidence level of the maternal assignment was rather low. The estimated variance of the genetic effect was infinitesimal, at 4.9E-05, and the estimate for the residual was 25.07. Consequently, differences among the individuals were imperceptible (Table S2). The maximum predicted breeding value was 0.00003 cm for hc1485, with a DBH 29.7 cm, and the minimum estimates was -0.00002 cm for hc1442, with a DBH of 6.6 cm. The correlation coefficient of the observed and predicted value was 0.94, which could be explained with predictions considering not only observations but also maternal information.
To separate the genetic effects with increased accuracy, the environmental effects were investigated to be accounted in the model of genetic testing. The spatial autocorrelation of phenotype was significant using Moran's I index (p < 7.4E-12). The results of seven types of spatial models, including exponential, gaussian, spherical, linear, rational, AR1, AR1 ⊗ AR1, and the basic animal model were compared to assess the model that would most adequately explain the phenotypic spatial structure of the plantation. The distribution of the phenotype was rasterized for use in the spatial analysis with AR1 ⊗ AR1 (Fig. 1). The rasterization was cautiously processed to prevent data loss due to adjacency between individuals, considering their irregular distribution. The integration of the animal and AR1 ⊗ AR1 model (animal + AR1 ⊗ AR1) yielded the highest fitness followed by the basic animal model, based on the log likelihood and AICc ( Table 3). The AICc of the animal + AR1 ⊗ AR1 model was largely distinct from those of the other models, highlighting its superior fitness. The fitness of the model with spatial AR1 was very low, suggesting the need for the use of a two-dimensional firstorder autoregression model. The fitness and residual of the animal model and animal + AR1 ⊗ AR1 model were confirmed in the present study (Fig. 2, Fig. 3). The residual was reduced in the animal + AR1 ⊗ AR1 model when compared with in the animal model.
The variance component of genetic factor was 9.086, and the heritability was estimated as 0.344 ( Table 4). The accuracy of the genetic testing improved following the elimination of the environmental effects, with a heritability estimate of 0.224 for L. kaempferi DBH [1]. The predicted breeding value (PBV) of each individual also reflected the enhanced power of explanation of the animal + AR1 ⊗ AR1 model (Table S2). The maximum and minimum PBV were 5.655 cm and -5.937 cm, respectively, showing the apparent difference between the animal model and the animal + AR1 ⊗ AR1 model. The Spearman correlation coefficient (r) between the observation and the predicted breeding value using the  animal + AR1 ⊗ AR1 model was 0.87 (Fig. 4). The weakened correlation was considered to be due to the consideration of the environmental effects when evaluating the phenotypes. The r between the prediction by the animal model and the animal + AR1 ⊗ AR1 model was 0.89, revealing the dissimilar estimates obtained by the two models. The breeding values predicted by the two models were compared to confirm the frequent changes in rank (Table S2) reflecting the practicality of spatial analysis approaches in genetic testing. PBVs of mother trees were estimated from the fitted model (Fig. 5), with the highest value obtained in KW02, followed by JB10, NO05, NO03, and GB25 (Table S3). The PBV of GB23 was the lowest followed by JB08, CN15, and KW56. In the case of KW02, which had the highest PBV, five trees were assigned to its family in HC (Table  S1, Table S3). All the five individuals, including hc1327, hc1427, hc1485, hc1538, and hc1549, had positive PBVs of 1.507 cm, 4.846 cm, 3.621 cm, 2.171 cm, and 0.117 cm, respectively. JB10, which had the second highest PBV, was the mother tree of hc1335, which showed the highest PBV among the offspring, at 5.655 cm. The progenies of NO05 were hc1436 and hc1552, with PBVs of 3.971 cm and 2.604 cm, which were the 6th and 24th highest PBVs ). GB23, the mother tree with the lowest PBV, had two progenies, hc1509 and hc1516, with PBVs of -3.812 cm and -5.937 cm, respectively. The latter had the lowest breeding value among the 188 trees. The PBVs of mother trees were considered to reflect the PBVs of their progenies, however the accuracy of analysis was limited due to the insufficient number of mother trees and their progenies.

Genetic variation in plantation of seed orchard crops of Larix kaempferi
Forest trees, which are long-lived plants that undergo outcrossing frequently, generally have high genetic diversity, excluding in the cases of endangered or rare species. The L. kaempferi seed orchard offspring was confirmed to retain considerable genetic variation in the present study. High genetic diversity was detected in the limited natural distribution range of L. kaempferi in Mt. Fuji even on fragmented habitats [31]. Frequent gene flow between populations was inferred from the low level of genetic differentiation despite the geographical distance in both Japanese larch and European larch [31,34].
Studying the change in the genetic variation using seed orchard crops is necessary to detect and preclude decline in genetic diversity with increased inbreeding [35][36][37]. Similar levels of genetic diversity in seed orchards and their crops are frequently reported, dispelling concern over narrowed genetic base in tree improvement programs. The genetic diversity of the progeny trial corresponded to that of a naturally regenerated stand in Picea sitchensis and Picea abies [35,36]. In addition, the genetic diversity of a plantation with seed orchard crops was not affected by a tree improvement program in Cryptomeria japonica [37].
However, genetic variation in larch should be assessed cautiously. Genetic variation in larch is influenced by the mating system, reproductive characteristics, and changes in descendants. Deficit of heterozygosity has been detected in Larix sibirica and Larix gmelinii, and attributed to inbreeding [38]. Furthermore, in a study on the larch mating system, the distance of pollen dispersal was reported to be relatively short [32,39,40], and the genetic structure of European larch over its entire distribution was more obvious than those of other conifers [41].
In summary, larch generally retain considerable genetic variation, similar to most forest trees. However, the species could be experience decreased genetic variation due to inhibited gene flow or severe fluctuation in flowering and seed production activities [6,8,42].
The first prerequisite in tree improvement is understanding the genetic variation in the natural population and in the selected breeding population [43]. Genetic variation is the key factor influencing natural and artificial selection activities [44] related with the coefficient of variance of the genetic factor. Genetic variation determines the genetic gain and the potential of improvement by influencing adaptability and resistance to biotic or abiotic stress [45]. Studying tree genetic variation is also essential for sustainable forest management and conservation [46], which could enhance the understanding of variation in adaptable traits including growth, wood characteristics, and resistance to disease or drought [47]. Understanding genetic variation is particularly essential prior to genetic testing in L. kaempferi plantation considering the reproduction and mating characteristic of the species. The additional introduction to enrich the gene pool of larch would assure the genetic diversity of breeding materials in long-term tree improvement program.

Pedigree reconstruction of individuals in the L. kaempferi plantation
Using an increased number of polymorphic markers would be necessary to enhance the confidence level of the maternity assignment in L. kaempferi. A study investigated the accuracy of pedigree reconstruction using a microsatellite marker as the principal marker in the field of research [48]. The success rate of pedigree reconstruction was improved with an increase in the number of markers in the pedigree reconstruction analysis for European sturgeon (Acipenser sturio) [49]. The low confidence level of parental assignment was thought to be related with allele frequency in a Abies procera study, where the microsatellite marker revealed insufficient variation [48]. The separate experiment and genotyping of the plus trees and the offspring was the other plausible reason of a decrease in the confidence in this study. The simultaneous analysis of the overall sample would reduce the scoring error of genotyping [50].
The application of molecular biology technologies in tree improvement can be classified into two approaches [12,48], including association between target traits and markers, such as the genetic dissection or genomic selection, and enhancing tree improvement efficiency using DNA markers, which includes pedigree reconstruction analysis [51,52]. The development of the DNA markers with high resolutions increases the usefulness of the pedigree reconstruction [53]. Furthermore, it is expected that the confidence of assignment would be enhanced by the rapid development and commercialization of automated and high-throughput genotyping approaches, including SNP marker analysis [48].
The small number of assigned progenies by parent has been discussed as one of the main challenges of genetic testing with pedigree reconstruction [15], and it impedes researchers' ability to obtain sufficient numbers of samples required to establish statistical significance in genetic testing. Therefore, forward selection is considered desirable compared to the backward selection in the case. The latter has been recommended considering the relative ease of propagation following tree improvement acceleration [16]. However, forward selection was suggested in this study due to the possibility of deterioration of accuracy when performing genetic testing on mother trees.

Genetic testing in the plantation of L. kaempferi
Genetic gain following selection is expected in the genus Larix considering previous genetic testing [2,17,[54][55][56][57]. The heritability of L. decidua height was 0.23 in the 25-35 age group [17], and the heritability values for height and DBH estimates of L. kaempferi were 0.549 and 0.224, respectively in the genetic testing of artificial crosses [1]. The infinitesimal genetic variance in the analysis with the animal model highlighted the need for isolating the environmental factor. The plantation was located in a highly heterogeneous site, which hampered the analysis of environmental effects. Prior pedigree reconstruction studies have mainly been carried out in experimental forests under intensive management, with a focus on the development of analytical methodologies [14-17, 21, 22, 58]. In addition, few substantive genetic testing studies have been conducted in the genus Abies [15,16,48], and the genetic resources of L. decidua have been evaluated based on locations [17]. The effective analysis of environmental effects is a major obstacle in genetic testing. Notably, the application of pedigree reconstruction in practical genetic testing has been accomplished in the study of genus Abies by supplementation with polymorphic SNP marker and spatial analyses.
The individual trees and their mother trees were evaluated with PBVs based on the animal + AR1 ⊗ AR1 model in this study. Considering the spatial autoregression as the products of innumerous biotic and abiotic factors proved useful in the plantation genetic testing, in combination with the pedigree reconstruction analysis. The inferior fitness of the spatial model except AR1 ⊗ AR1 highlighted the importance of selecting the adequate spatial structure model beyond the consideration of the existence of spatial autocorrelation. The enhancement of genetic testing via spatial analysis with AR1 ⊗ AR1 has been reported [26][27][28][29][30]. Genetic testing with the model considering AR1 ⊗ AR1 structure in pedigree reconstruction analysis has been used in family selection in a Abies bornmülleriana progeny trial [16].
It is impossible to control the number of samples by family in genetic testing in plantations, contrary to the case in traditional progeny testing. Only one progeny was assigned to the multiple family in the present study as in the pedigree reconstruction analysis in L. occidentalis and A. nordmanniana [15,17]. The spatial analysis was considered to be essential in the selection without the experimental design including the replication of the family. Spatial analysis is also considered to be useful in the progeny trials [23,27,29], in particular, under high environmental heterogeneity, for example in mountainous sites in Korea [59].

Conclusions
Genetic testing was conducted in a L. kaempferi plantation to establish a systematic scheme that could supplement progeny trials in tree improvement activities. Pedigree reconstruction and phenotypic spatial distribution were applied in the analysis of genetic and environmental effects on tree growth. The pedigree reconstruction was performed for individual trees following the assessment of genetic variation to identify the adequacy of breeding population. The genetic testing was implemented considering the spatial autoregression to dissect the environmental factors influencing the growth characteristics of a plantation without experimental design.
Several factors should be taken into account in the course of interpreting the results of the present study. The planting and management records of the study site would save the additional labor for the verification of the data. In addition, assessing the genetic diversity is necessary for preventing a decrease in genetic variation and inbreeding depression in advanced generations. Furthermore, appropriate selection strategies should be adopted based on the limited sample sizes obtained in the progeny. Forward selection was suitable in the present study considering the small progeny size per mother tree. The effective application of the pedigree reconstruction in genetic testing requires the combined analysis of phenotypic spatial distribution as well as the influence of numerous environmental factors.

Study sites and sampling
The L. kaempferi plantation of seed orchard crops, owned by the state, was selected as the study site based on plantation history. The reforestation was undertaken on the preexisting larch plantation in Hongcheon-gun (HC), Gangwon-do (37°35′24.0'N, 127°53′24.0'E) with 3,000 larch seedlings / ha over an area of 1.3 ha in 1990. The remaining number of larch trees at present was 271 including the preexisted old trees. Therefore, the study materials were selected through screening to take the field circumstances into consideration preceding the analysis and in accordance with the national forest management record. The study complied with the guidelines of the National Institute of Forest Science (NIFoS) and the legislation of plant ethics in South Korea.
The growth data and the empirical yield  [60]. The outlying 83 trees with DBH > 30 cm were excluded from further analyses.
The overlapping existence of different age group of trees was sustained by deviation from the normal distribution pattern of the growth data (Fig. 6). The rest of the 188 were considered to be the subsequent generation from the seed orchard based on the scarce chance of larch seedling occurrence in the larch forest. The screening based on DBH was possible in the present study because the differences in growth between the excluded and remnant trees were rather distinct. One-hundred and forty-six trees from a five-year-old seed orchard planted in Hwaseong-si (HS), Gyeonggi-do, were additionally included in the analysis for the comparison with the HC plantation.
Three types of data were collected for the phenotype, genotype, and spatial distribution analyses of the individuals. Tree height and DBH, cambium or leaf samples, and GPS coordinates (Trimble Geo 7X, Trimble Inc., Sunnyvale, CA, USA) were collected simultaneously.

Genetic variation and pedigree reconstruction of L. kaempferi plantation
The isolated DNA of 334 L. kaempferi trees were genotyped using seven microsatellite markers, which were used in a previous analysis on the mother trees in the seed orchard [61,62] (Table 5). The genotype data of the mother trees used in the present study were assessed by the National Forest Seed and Variety Center (NFSV). The polymerase chain reaction (PCR) was carried out using FAM-labelled forward primers of each marker followed by genotype analysis using an ABI 3730xl Genetic analyzer and GeneMapper 4.0 (Applied Biosystems, Foster City, CA, USA).
Genetic variation in the larch plantation of seed orchard crops was estimated based on the number of alleles (N A ), the number of the effective alleles (N E ), observed heterozygosity (H O ), expected heterozygosity (H E ), and fixation index (GenAlEx 6.5) [63]. The data on the mother trees were provided by the NFSV [6].
The maternity analysis was conducted for the pedigree reconstruction of the L. kaempferi plantation. The individuals in the plantation were assigned to the mother trees in the seed orchard using CERVUS 3.0 based on the maximum likelihood approach [64,65]. The confidence levels were set at 95% (restricted) and 80% (relaxed), with   default values in the program. The error rate of genotype identification was specified as 0.01 and the confidence level was decided based on 10,000 simulations.
Genetic testing using pedigree reconstruction and spatial distribution of phenotype The growth correlations among adjacent individuals were observed in the study site. Spatial autocorrelation was tested based on DBH using Moran' I prior to genetic testing. As the existence of the spatial autocorrelation of DBH was identified, the spatial model was explored using seven spatial structures, including AR1, exponential, Gaussian, spherical, linear, rational, and AR1 ⊗ AR1 [66,67] to verify the utility of AR1 ⊗ AR1. The best fitted model was selected by comparing the log-likelihood values and the second-order Akaike Information Criterion (AICc) of the eight linear mixed models, including the animal model [68,69]. The rasterization of the plantation was conducted to form 450 rows and 450 columns for AR1 ⊗ AR1 analysis. The accuracy of the genetic testing was improved using the AR1 ⊗ AR1 model, which is the Kronecker product of the autoregression matrix for rows and columns [70,71]. y = Xb + Zu + ξ + η. where y is the vector of data; b is a vector of fixed effects; u is a vector of random effects (genetic or family effect); ξ is spatially dependent residuals with a variancecovariance matrix of σ 2 ξ [AR1(ρ col ) ⊗ AR1(ρ row )] and η is spatially independent residuals (nugget effects). AR1(ρ col ) and AR1(ρ row ) are the first-order autoregressive correlation matrices in column and row respectively.
The breeding values of mother trees and the individuals in the plantation were predicted by the best fitted model using R program [72]. Existence of spatial autocorrelation was tested with package ape [73] and the rasterization was conducted using sp [74] and raster [75]. The analyses of models with spatial structures were conducted with package nlme [66], breedR [67] compared using MuMIn [76]