Impact of Mixing on the Structural Diversity of Serbian Spruce and Macedonian Pine Endemic to Relict Forest Communities in the Balkan Peninsula

: The aim of this paper is to analyze the effect of different degrees of mixing on the diversity structure in stands left to spontaneous development. The research included two communities of species endemic to the Balkan Peninsula—the Serbian spruce ( Picea omorika Panˇci´c Purk.) and the Macedonian pine ( Pinus peuce Griseb). Data from eight sample plots were used in the research. The changes in diameter and height structure, spatial arrangement of trees, and diameter differentiation were analyzed. The analyzed parameters of structural diversity show relatively low to moderate values. Results showed an increase in mixing was reﬂected in the width and shape of distributions. A spatial analysis of stands with a higher degree of mixing showed a tendency towards a random to regular distribution of individuals, in contrast to stands with a lower degree of mixing which showed a tendency towards a clump distribution. The pronounced species’ dimensional and spatial diversity conﬁrms their importance to the condition of modern forest management. Signiﬁcant differences in the change of structure are shown by stands with a share of admixed species of above 20% by volume. The obtained results refer to stands left to spontaneous development, suggesting than an active research and management approach must be assumed to realize the goal of protecting rare forest ecosystems.


Introduction
The diversity of vegetation of the Balkan Peninsula, in a broader context, can be explained by its position where diametrically opposed climatic influences meet and affect the area [1]. The most important relict-endemic forest ecosystems of the Balkan Peninsula are the communities of Serbian spruce (Picea omorika) and Macedonian pine (Pinus peuce). Anthropogenic impacts in the past have reduced the range of Serbian spruce to small and inaccessible habitats along the middle course of the Drina River (Bosnia and Herzegovina, Serbia) and Macedonian pine to habitats of the highland region (Bulgaria, Serbia, Macedonia, Montenegro, Albania, and Greece) ( Figure 1) [2]. The character of habitats in which these communities develop, the passive approach to their protection, as well as the absence of economic interests, have resulted in poor research of these ecosystems [3,4].
Structural parameters have been assessed as very good indicators of biodiversity and are recommended as its direct measure [4,5]. Whittaker [6] defines different categories of biodiversity. The proposed division at the stand level implies a variety of tree positions, tree species, and their dimensions [7][8][9]. The study of these characteristics does not neglect the significant three-dimensional structure of stands [8,[10][11][12][13].
Today, unlike in the past when monocultures were of great importance in Europe, ever more attention is paid to heterogeneous mixed stands [4]. The heterogeneity of natural spontaneously developed stands is higher than in those managed by man [4,14].  [20]. The sample plots of Serbian spruce are located in the Bosnian municipality of Višegrad (Veliki Stolac) and the Macedonian pine in the Serbian municipality of Tutin (Beleg).  [20]. The sample plots of Serbian spruce are located in the Bosnian municipality of Višegrad (Veliki Stolac) and the Macedonian pine in the Serbian municipality of Tutin (Beleg).

Research Area
During 2019, four sample plots (SP) were established in the stands of Piceaomorika and another four in the stands of Pinus peuce (Figure 1). The sample plots were located The mixed stands of Piceo-OmorikaeAbieti fagetumČol.1965 occupy the edges and lower areas of the investigated locality (Veliki Stolac), where they spatially lean on mixed communities of spruce, fir, beech, and Austrian pine. The altitude ranges from 1300 m (SP 2) to 1370 m (SP 3). The exposure is north, and the slope is 30-42 • . The parent rock is limestone with shallowly developed black soils on it [21]. The average annual rainfall is 977 mm with an average annual temperature of 5 • C.
The Macedonian pine stands at the investigated locality (Beleg) form a vegetation belt above the forests of spruce and the forests of spruce and Macedonian pine Piceo-Pinetum peuces Lakš. 1965 [20]. They are located at an altitude ranging from 1825 m (SP 5) to 1910 m (SP 6). The exposures are different (SP 5 W-NW; SP 6 N-NW; SP 7 E; SP 8 N), and the slope of the terrain is 15-30 • . The parent rock is limestone, and the soil is shallow dystric humus-silicate soil. The average annual rainfall is 749 mm, and the average annual temperature is 5.2 • C.
The share of other tree species in the stands of Serbian spruce (SP 1-4) and Macedonian pine , expressed as % of volume per hectare (V × ha −1 ), is shown in Table 1.
The stands included in this research have the character of strictly protected natural communities that have not been exposed to direct anthropogenic influences that could have compromised their natural character. The preservation of these relict-endemic ecosystems has been contributed to by the difficult-to-access terrain, distance of the locality from populated areas, and lack of forest roads.

Data Collection
The positions of the sample plots were recorded with a GPS device. The diameters at breast height (dhb) and the total heights (h) of all trees with a diameter greater than 10 cm were measured in the sample plots. The trees were marked with numbers at breast height and their coordinates were recorded. The coordinates were recorded using a WILD T05 theodolite and a laser distance meter from the defined angles of the sample plots.

Data Processing
Complex structural analysis based on the spatial relationships of neighboring trees requires the use of the edge-correction method [22]. In order to obtain unbiased data, the plus-sampling method [23,24] was used. Trees from 44.7 m × 44.7 m (0.2 ha) subsamples from sample plots' centers were used to avoid the edge effect. All parameters of spatial analysis were calculated on the basis of a defined subsample. The mixing effect by the number of trees in the sample plots was defined by the measure of evenness E [25], the percentage ratio of (H) (the species diversity index [26]), and the maximum diversity (H max = lnS), where S is the number of different tree species in the sample plot. The mixing effect is also defined according to the percentage share of species by volume (Table 1). Tree volumes were calculated according to the local models V = f (dbh, h) for the analyzed tree species [27][28][29][30][31]. Volume was calculated for each individual tree based on input values of diameter and total height according to two-way volume tables. The volume thus obtained on the sample plots was calculated per hectare.
The measure of spatial segregation of species used is the mingling index M [32]. A structural group of four neighboring trees was used to determine this index [8,[33][34][35]. The index of structural groups is calculated according to the following formula: At the level of the sample plot, index M was calculated as the arithmetic mean of the individual values of the M i index of structural groups. The designations used in the function have the following meaning: i is the ordinal number of the reference tree (1 . . . N), N is the number of trees in the plot, j is the number of the neighboring tree (1 . . . 4), n is the observed number of neighboring trees (4), and m j is a parameter (1 if the reference and neighboring species are different, 0 if they are the same).

Horizontal Structure
The Clark-Evans aggregation index R [36] was used to analyze the spatial pattern of trees in the sample plots. The index was used according to the adjusted calculation method [37] of real (rA) and expected average distance between individuals (rE), as follows: Letters in the formula have the following meaning: A is the area (m 2 ), N is the number of trees, P is the perimeter (m), and r i is the distance from the i-th tree to the nearest neighbor.
The spatial pattern of trees was determined based on the Ripley's K(r) function [38]. The edge-correction effect [23,[39][40][41] implies an unbiased way of calculating the K(r) function: Interpretation of the results was performed using the L(r) function, in a form that in the simplest way shows the change in the L value with an increasing distance r. In the graphical representation of the spatial pattern of trees, we used the adapted form of the L (r) function [40,42].
The letters in the previous two functions have the following meanings: A is the area (m 2 ), N is the number of trees, u ij is the distance between the i-th tree and the j-th tree, I r (u ij ) is the indicator of the function that equals 1 if u ij ≤ r (or 0 vice versa), w ij is the edge-correction factor for the i-th tree and its neighboring j-th tree, and r is the analyzed reference distance. The height structure by sample plots is shown graphically based on the percentage of trees (N%) by height classes of 3.0 m width. In order to simplify the comparative analysis of the height structure of individual sample plots, statistical measures of central tendency (average height), dispersion (standard deviation, variation width), and distribution (skewness, kurtosis) were used.
Species profile index A [43] was used as a unique measure of the height structure of the sample plots, based on the adapted method of calculating diversity index H [26]. The comparatively significant relative value of the species profile index A rel was obtained from the ratio of absolute A and maximum A max value, as follows: The letters used have the following meanings: S is the number of tree species, Z is the number of vertical layers (3), p ij is the share of species in vertical layers (n ij /N), n ij is the number of trees of a certain species in the vertical layer, and N is the total number of trees.

Diameter Structure
The diameter structure of trees in the sample plots is graphically shown on the basis of the percentage of trees (N%) by diameter classes of 5 cm width. In order to facilitate the comparative analysis of diameter distributions in the sample plots, the statistical parameters listed for the vertical structure were used.
The diameter differentiation index Td [32] was used to analyze diameter diversity. The index quantifies the diversity of reference tree diameters and their nearest neighboring trees [5,32,44]. A structural group of four neighboring trees was used to determine this index [8,32,33,45]. Structural group indices were calculated according to the following equation: At sample plot level, the Td index is calculated as the arithmetic mean of the individual values of the Td i index of structural groups. The letters used in the formula have the following meaning: i is the number of the reference tree (1 . . . N), j is the number of the nearest tree (1 . . . n), n is the number of nearest trees (4), and N is the total number of trees in the plot.
In addition to the above parameters and indices, the Gini coefficient (GC) was used to estimate the structural unevenness of the sample plots [46,47]: The letters used have the following meaning: j is the rank of trees according to the diameter expressed in ascending order from 1 to n, n is the total number of trees in the sample plot, ba j is the basal area of the j-th tree by rank (m 2 ).
The hyperbolic tangent index S [48] was used to define the relationships between the species in the sample plot. The relationship between the diameters of reference trees of admixed species and dominant species that are spatially closest to the observed reference trees was analyzed. A structural group with four neighboring trees was used [48]. The structural group index S i is calculated as follows: The sample plot index S was calculated as the arithmetic mean of the structural group indices S i without edge-correction [22,46]. The letters in the previous formula have the following meaning: i is the number of the reference tree of the admixed species (1 . . . N), j is the number of the nearest dominant tree (1 . . . k), k is the number of nearest dominant trees to the reference tree (4), m i is the reference tree (observed characteristic-dbh i ), m j is the neighboring tree (observed characteristic-dbh j ), factor a ranges from 0.0 to 1.5 (the most appropriate is 1.0), and N is the total number of trees of admixed species in the sample plot.

Statistical Analysis
Testing of the distribution of the investigated characteristics of trees in relation to their normal distribution was performed using the Shapiro-Wilk test of normality. The empirical data of the distribution and normal distribution were compared. The Levene test was used to test the homogeneity of variance. A comparison of the empirical distributions of the analyzed characteristics was done based on the Mann-Whitney U test. The acceptance (rejection) of the hypotheses about the normality of data distribution, i.e., their mutual differences, was tested at the significance level of 0.05. The analysis of the relationship between the used structure parameters was performed using the Spearman's rank-order correlation. The testing of the spatial patterns was conducted with 199 Monte Carlo simulations of the null hypothesis of complete spatial randomness. The parts of the L (r) function above, below, and within confidence the intervals indicate clumped, regular, and random spatial patterns.
Different software packages [49,50] were used in the data processing, analysis, and testing. The spatial analysis of trees was performed using the "spatstat" [51] package in the R-3.5.1 program [52].

Dendrometric Characteristics of the Investigated Forest Communities
Forest communities of relict-endemic species, according to the basic dendrometric values (Table 2), have the following characteristics. In the stands of Serbian spruce, the mean diameter per basal area ranges from 24.7 (SP 2) to 29.9 cm (SP 3), and the mean height ranges from 25.8 (SP 2) to 27.9 m (SP 3). The number of trees ranges from 810 (SP 3) to 1285 pcs × ha −1 (SP 4). The basal area values range from 52.9 (SP 1) to 65.1 m 2 × ha −1 (SP 4), and volumes from 649 (SP 1) to 794 m 3 × ha −1 (SP 4). There is an obvious trend of increasing volume with an increasing degree of mixing of stands of this community (Va%). In the Macedonian pine stands, the mean diameter per basal area ranges from 31.8 (SP 5) to 44.2 cm (SP 8) and the height ranges from 17.6 (SP 5) to 20.7 m (SP 8). The number of trees ranges from 326 (SP 8) to 1015 pcs × ha −1 (SP 5), the basal area ranges from 50.1 (SP 8) to 80.0 m 2 × ha −1 (SP 5), and the volume ranges from 511 (SP 8) to 700 m 3 × ha −1 (SP 5). There is an evident decrease in the number of trees and volume, as well as an increase in the mean diameter and height with an increase in the mixing of Macedonian pine stands.  Differences at the community level are expressed in all parameters. The average values of heights (HQMD; HDo) and slenderness (HQMD/QMD) are higher for 35% (32%) and 100%, and average value of diameters (QMD; Do) are lower for 48% (52%) in the community of Serbian spruce. The average number of trees (1079 pcs × ha −1 ) and volume (728 m 3 × ha −1 ) are higher by 82% and 28%, respectively, in the Serbian spruce community compared to the Macedonian pine community. In contrast, the average value of basal area in the Macedonian pine community (63.2 m 2 × ha −1 ) is higher by 8% than in the Serbian spruce community (58.4 m 2 × ha −1 ). The different ratios of volume and basal area in these two forest communities are a consequence of large differences in slenderness, i.e., the height structure of the stands.

Spatial Pattern of Tree Distribution in the Investigated Stands
The average distance between neighboring trees (Table 3) in the community of Serbian spruce ranges from 1.44 m (SP 2) to 1.84 m (SP 3). The aggregation indices R have values close to the reference random distribution. The test results indicate that only the trees in the sample plot with the highest degree of mixing (SP 4) have a regular distribution (p < 0.05). The average distance between neighboring Macedonian pine trees ranges from 1.14 m (SP 1) to 3.41 m (SP 4). The aggregation indices R differ significantly in the sample plots. The testing indicates that only the trees from sample plot 7 are randomly distributed in space. In the sample plots with the lowest degree of mixing (SP 5, SP 6), a statistically significant clumping was found (p < 0.001, p < 0.01), while the trees in the sample plot with the highest degree of mixing (SP 8) had a regular distribution (p < 0.05). The spatial pattern of tree distribution in the studied communities was determined based on the Ripley's L(r) function ( Figure 2). In the community of Serbian spruce in the first two sample plots with a lower degree of mixing (SP 1, SP 2), the distribution of trees with a change in distance has different values. Based on the conducted Monte Carlo simulations of the null hypothesis of complete spatial randomness, at distances of up to 2.0 m (SP 1) the trees are randomly arranged, and at over 2.0 m the arrangement has a clump character. In sample plot 2, the random distribution of trees is at distances of up to 4.0 m, from 4.0 to 8.0 m the distribution is clumped, and at over 8.0 m it is at the upper limit of statistical randomness. In the other two sample plots (SP 3, SP 4) with a higher degree of mixing, the distribution of trees is within the limits of statistical significance of random distribution at all distances.
2.0 m (SP 1) the trees are randomly arranged, and at over 2.0 m the arrangement has a clump character. In sample plot 2, the random distribution of trees is at distances of up to 4.0 m, from 4.0 to 8.0 m the distribution is clumped, and at over 8.0 m it is at the upper limit of statistical randomness. In the other two sample plots (SP 3, SP 4) with a higher degree of mixing, the distribution of trees is within the limits of statistical significance of random distribution at all distances.
In the Macedonian pine community, in sample plots with a lower degree of mixing (SP 5, SP 6), the distribution of trees has different values with a change in distance. At distances of up to 3.0 m (SP 5) and 2.5 m (SP 6), the trees are clumped in space. As these distances increase, the trees have a random arrangement. In the sample plots with a higher degree of mixing (SP 7, SP 8), the arrangement of trees is within the limits of statistical significance of random distribution at all distances, with a significantly wider confidence level. In the stands of Serbian spruce the values of Ripley's L(r) function increase with an increasing distance, while with increasing mixing they have smaller absolute values at the observed distances. As the mixing increases, the spatial pattern of trees tends to change from a clumped to a random distribution. In the stands of Macedonian pine, the functions have slightly more complex relationships, but it is clear that with increasing mixing, the arrangement of trees has a trend from a clumped to a regular distribution.  In the Macedonian pine community, in sample plots with a lower degree of mixing (SP 5, SP 6), the distribution of trees has different values with a change in distance. At distances of up to 3.0 m (SP 5) and 2.5 m (SP 6), the trees are clumped in space. As these distances increase, the trees have a random arrangement. In the sample plots with a higher degree of mixing (SP 7, SP 8), the arrangement of trees is within the limits of statistical significance of random distribution at all distances, with a significantly wider confidence level.
In the stands of Serbian spruce the values of Ripley's L(r) function increase with an increasing distance, while with increasing mixing they have smaller absolute values at the observed distances. As the mixing increases, the spatial pattern of trees tends to change from a clumped to a random distribution. In the stands of Macedonian pine, the functions have slightly more complex relationships, but it is clear that with increasing mixing, the arrangement of trees has a trend from a clumped to a regular distribution.

Vertical Structure
The distributions of trees by height classes in the sample plots in the investigated forest communities are graphically presented in Figure 3. Based on the testing of empirical height distributions (Mann-Whitney U test), a statistically significant difference (p < 0.05) was found between the distributions in the sample plots.
The arithmetic mean height (hA) in the Serbian spruce stands ranges from 22.4 (SP 4) to 27.1 m (SP 3). The variation width of heights (hV) is even across the sample plots and ranges from 29.4 m (SP 1) to 31.7 m (SP 4), with a higher standard deviation in sample plots with a higher degree of mixing (SP 3, SP 4). The height structure in all sample plots shows a negative (right) skewness (α3 < 0), while the values of the kurtosis coefficient α4 < 3 in SP 1, SP 3, and SP 4 indicate kurtosis from above, and in SP 2 lateral kurtosis (α4 > 3). The concentration of the number of trees of 30% (SP 1) and 40% (SP 2) in one height class and lower values of deviation from the arithmetic mean indicate a slightly different structure compared to the sample plots with a higher degree of mixing (SP 3, SP 4). The height to 27.1 m (SP 3). The variation width of heights (hV) is even across the sample plots and ranges from 29.4 m (SP 1) to 31.7 m (SP 4), with a higher standard deviation in sample plots with a higher degree of mixing (SP 3, SP 4). The height structure in all sample plots shows a negative (right) skewness (α3 < 0), while the values of the kurtosis coefficient α4 < 3 in SP 1, SP 3, and SP 4 indicate kurtosis from above, and in SP 2 lateral kurtosis (α4 > 3). The concentration of the number of trees of 30% (SP 1) and 40% (SP 2) in one height class and lower values of deviation from the arithmetic mean indicate a slightly different structure compared to the sample plots with a higher degree of mixing (SP 3, SP 4). The height structures in these sample plots do not have such clearly expressed maximums and have a slightly larger variation width. In the Macedonian pine stands, hA values are lower than in the stands of Serbian spruce. They range from 15.5 (SP 5) to 20.5 m (SP 8) and increase with the degree of mixing in the stand. Thestandard deviation hV also shows higher values with increased mixing. In all Macedonian pine sample plots, height structures have a negative (right) skewness (α3 < 0), while in SP 5, SP 6, and SP 8 α4 < 3, and in SP 7 α4 > 3. Like in the Serbian spruce stands, stands with a lower degree of mixing (SP 5, SP 6) are clearly separated from stands with a higher degree of mixing (SP 7, SP 8) in terms of height structure.
The differentiation of the vertical structure is also expressed on the basis of the results of vertical profile fulfillment Arel (Table 4). In the stands of Serbian spruce, Arel values range from 43.5% (SP 1) to 58.9% (SP 4), and in the Macedonian pine stands from 45.8% (SP 5) to 68.6% (SP 8). The growing trend of Arel in stands of both tree species is positively correlated with the increase in the degree of mixing. The increase in mixing directly affected the complexity of vertical form of the investigated stands.  In the Macedonian pine stands, hA values are lower than in the stands of Serbian spruce. They range from 15.5 (SP 5) to 20.5 m (SP 8) and increase with the degree of mixing in the stand. Thestandard deviation hV also shows higher values with increased mixing. In all Macedonian pine sample plots, height structures have a negative (right) skewness (α3 < 0), while in SP 5, SP 6, and SP 8 α4 < 3, and in SP 7 α4 > 3. Like in the Serbian spruce stands, stands with a lower degree of mixing (SP 5, SP 6) are clearly separated from stands with a higher degree of mixing (SP 7, SP 8) in terms of height structure.
The differentiation of the vertical structure is also expressed on the basis of the results of vertical profile fulfillment A rel (Table 4). In the stands of Serbian spruce, A rel values range from 43.5% (SP 1) to 58.9% (SP 4), and in the Macedonian pine stands from 45.8% (SP 5) to 68.6% (SP 8). The growing trend of A rel in stands of both tree species is positively correlated with the increase in the degree of mixing. The increase in mixing directly affected the complexity of vertical form of the investigated stands.

Diameter Structure
Based on the testing of empirical diameter distributions (Mann-Whitney U test), a statistically significant difference (p < 0.05) was found between the sample plots. In the Serbian spruce stands, the arithmetic mean diameter (dbhA) ranges from 23.6 (SP 4) to 28.7 cm (SP 3). A significantly higher variation width (dbhV) and standard deviation were found in the sample plots with a higher degree of mixing (SP 3, SP 4). The diameter distribution in the sample plot with the lowest degree of mixing has a negative skewness (α3 = −0. According to the shape of distributions (Figure 4), the sample plots with a lower degree of mixing (SP 1, SP 2) show a narrower diameter distribution and a more pronounced concentration of trees within one diameter class. The shapes of distributions focused on thinner diameter classes indicate pronounced processes of selection and differentiation of higher diameter classes. Tree differentiation (SP 3, SP 4) becomes more intense with an increase in mixing in stands (low right part of distribution).  (Figure 4). The forms indicate the processes of selection (SP 5), initiated (SP 6), and pronounced (SP 7, SP 8) tree differentiation. The diameter structure in all sample plots has a clearly expressed single maximum. The appearance of several less pronounced peaks along the variation width indicates a more complex structural form of these stands compared to the classical structure of even-aged stands.

Parameters of Diameter Diversity
The spatial measure of mixing M clearly distinguishes between two categories of species segregation in stands ( Table 5). The first category with a higher degree of segregation is below 0.20 (SP 1, SP 2, SP 5, SP 6) and the lower one is above 0.  (Figure 4). The forms indicate the processes of selection (SP 5), initiated (SP 6), and pronounced (SP 7, SP 8) tree differentiation. The diameter structure in all sample plots has a clearly expressed single maximum. The appearance of several less pronounced peaks along the variation width indicates a more complex structural form of these stands compared to the classical structure of even-aged stands.

Parameters of Diameter Diversity
The spatial measure of mixing M clearly distinguishes between two categories of species segregation in stands ( Table 5). The first category with a higher degree of segregation is below 0.20 (SP 1, SP 2, SP 5, SP 6) and the lower one is above 0.20 (SP 3, SP 4, SP 7,  M-spatial mixing index, Td-diameter differentiation index, GC-Gini coefficient of homogeneity, S-hyperbolic tangent index of admixed species.

Spearman's Correlation Coefficient ρ
Based on the established relationship between stand characteristics and mixing, a correlation analysis was performed using the Spearman's correlation coefficient (Table 6). A strong positive correlation was found between the mixing index M and the aggregation index R ( = 0.87, p < 0.01), the mixing index M and the relative species profile index A rel ( = 0.86, p < 0.05), the aggregation index R and the relative species profile index A rel ( = 0.74, p < 0.05), and the diameter differentiation index Td and the Gini coefficient GC ( = 0.85, p < 0.01). A weak to moderately strong correlation, without statistical significance, was found among the other parameters of stand structure.

Differences in the Basic Characteristics between the Researched Stands and Other Stands of the Same Species
In previous studies of mixed stands of Serbian spruce in Serbia [53,54], lower values of all numerical (dendrometric) variables were found at a locality in the immediate vicinity of the investigated stands except for the number of trees per unit area. The structure of these stands is much more complex than in the investigated stands and ranges from uneven-aged stands to structures close to selection stands. The increase in total volume with increasing mixing was also confirmed by these studies.
Available research on Macedonian pine stands in Macedonia (Pelister) [27,55] and Serbia (Šar Mountain) [56,57], higher values of the mean diameter and height of trees were found, while the basal area and volumes are smaller than in the investigated stands resulting from the smaller number of trees. These differences are a consequence of different ecological conditions of stand development. The stands of Macedonian pine [27,[55][56][57] are located in optimal conditions of development [58], while the investigated stands are located on the border of the vertical (1825 to 1910 m above sea level) and horizontal (north) range of the species (Figure 1).

Effect of Mixing on the Spatial Pattern of Tree Distribution
The spatial distribution of trees in stands is influenced by numerous factors. It is most influenced by the management system, the origin and structure of stands, the composition of the species that build them, developmental stages, natural disturbances, habitat conditions and abiotic conditions to which they have adapted their development. In intensively managed stands, the distribution of trees tends towards a regular (random) distribution [59,60]. In stands that are not subjected to management treatments, active selection and differentiation processes influence the regulation of tree distribution with a tendency toward a regular distribution [40,41], i.e., random distribution [60,61], which was confirmed by this research. These pronounced processes in mixed stands of the studied species affected the random distribution of trees. The regularity of distribution is especially pronounced in trees with larger diameters [40,41], which is also the trend observed in this research (SP 7, SP 8).
The tendency towards clumping of trees at lower distances in the stands of Serbian spruce (SP 1, SP 2) can be explained by microhabitat conditions (form of the micro-relief conditioned by the parent rock, a mosaic shallow initial phase of the soil) and the relative absence of other species with different bio-ecological characteristics which would contribute to the process of intensive competition and differentiation. The formed homogeneous groups in the coves and at the crossings between ridges, which intersect the entire locality, condition the clumping of trees at relatively larger distances. The clumping of Macedonian pine trees at lower distances (SP 5, SP 6) can be explained by the specifics of the species development in high mountain climatic conditions. The relations of the spatial distribution of the trees can be divided into two groups; I. with a lower degree of mixing (SP 1, SP 2; SP 5, SP 6), where the trees tend to clump (group) and II. with a higher degree of mixing (SP 3, SP 4; SP 7, SP 8), where the trees tend to be randomly to regularly distributed. Wider distributions of dimension (Figures 3 and 4) indicate an increasingly pronounced influence of the separation of ecological niches of different species in mixed stands and the dimensional unevenness and intensive separation of taller and shorter trees into strata. This fact also affected the spatial arrangement of the trees.

Diversity Indices as an Appropriate Measure of Stand Differences (Td, GC)
The diameter differentiation index Td based on the researched references shows different values in different stand situations. The lowest values are typical of even-aged stands that are actively managed, 0.11 in taigas [62], 0.13-0.21 in pine cultures of different ages [41], 0.21 in an even-aged beech stand [8], 0.25 in young Douglas fir culture [8], 0.30 [59], and 0.34 [60] in even-aged stands of different forest types. Slightly higher average values of Td were found in uneven-aged stands of beech 0.27-0.42 [63], 0.36 [64], and 0.42 [8] in mixed stands of beech. In a series of permanent research plots in the Czech Republic [65][66][67][68] on a similar sample in protected mixed stands of spruce, fir and beech values were determined to be from 0.41 to 0.55 [67] and in managed mixed stands of beech of different origin they were from 0.36-0.49 [68]. In relict pure stands of Scots pine lower values (0.20-0.33) were found [65]. In research on mixed stands of Scots pine without management treatment [66], a trend of decreasing the degree of diameter differentiation (from 0.37-0.48 to 0.36-0.42) was determined in the analyzed period of 15 years on all permanent research plots. The inventory of different forest types in the Austrian Alps determined the average value of the Td index in stands of 0.37, and values significantly above this average (0.50) were found in two-story (multi-story) stands [69]. In the selected stands, the determined mean values of the Td index are 0.43 [64] and 0.42 [59]. The highest Td values are typical of virgin forests, and they range from 0.47 [70] to 0.76-0.78 [60].
By comparing the presented Td indices, the differentiation of diameters has characteristic higher or lower values in relation to the developmental and structural type of stands. These differences indicate, but do not determine, the stand type. In non-managed forests, the differences are the result of competition between species and within the species during development, and in forests that are regularly managed, they are a direct consequence of these relationships and the nature of applied management treatments.
The determined values of the index from 0.24 to 0.33 indicate a simpler type of stand, which only confirms the previous analysis of the height and diameter structure.
In the practical use of the Td index, it is important to pay attention to the distribution of index values of structural groups, as well as the use of appropriate measures of central tendency. The approximately normal values of the Td index distribution indicate the possibility of using mean values [62] as a representative parameter in the analysis of the relationship, which is the case here. For significant deviations from the normal distribution, the use of mode [62] or median [69] is adequate.
The Gini coefficient GC has the lowest values in young even-aged stands of 0.15-0.30 [71], 0.22 [60], and 0.25 in two-story stands [69]. Higher values are typical of mixed and uneven-aged stands of different compositions (e.g., 0.37-0.54 [19], 0.50-0.58 [71], 0.49-0.57 [72], and 0.35-0.52 [73]). Values above 0.60 [71] are characteristic of selective stands and virgin forests. The mean values in the selection stands of the Austrian Alps are 0.63 [59]. The highest values of the index were found in beech and spruce reserves and virgin forests of the southeastern Carpathians 0.69-0.71 [74], beech and hornbeam in the Czech Republic 0.67-0.75 [60], beech, fir, spruce in Bosnia and Herzegovina 0.67 [70], and beech in Serbia 0.45-0.52 [72]. The values presented in different stands justify the use of GC as an adequate measure of stand inequality, which also indicates but does not determine the stand type.
In relation to the distribution of the analyzed dimensions, the indices show certain regularities. Stands with a normal distribution have lower GC values compared to stands with the inverse J-shape [46,59] and regular distributions [46]. In this research, two stands with a statistically significant normal distribution (SP 1, SP 6) have the lowest GC values. This confirms the assumption that GC values increase with the deviation of the distribution from the normal shape [47,75]. The diameter distribution in uneven-aged stands is reflected in higher and more constant index values compared to even-aged stands [47]. In the research of different forms of pure and mixed stands of spruce, fir, and beech of the Eastern Carpathians [75] GC has been shown to be a useful tool in differentiating different structural types of stands. Regardless of the established justification and superiority in relation to other indices [46], the possibility of obtaining the same values in different structures (as with the previous index) limits its application [76]. That shortcoming requires the analysis of the flows of the Lorentz curve on which the index is based or the use of adequate additional parameters in determining the differences of the analyzed stands [76].
In the practical use of this index, as a useful tool in forest management, it is necessary to pay attention to the impact of the callipering threshold on its value [75], the type and size of the sample used [77], as well as the index value in the structures of transitional forms in the process of stand regeneration [47].

Impact of Admixed Species on Stand Structure
In order to determine the causes of changes in structure parameters with an increasing degree of mixing, the ratio of diameters between admixed and dominant species was analyzed in the investigated stands. The hyperbolic tangent index S in the stands of Serbian spruce has a growing trend in relation to mixing. The parameter shows that, with an increase in their share, admixed species have larger dimensions compared to the neighboring (competing) trees of dominant species. This fact can explain the almost regular changes in the analyzed parameters of the structure of Serbian spruce stands with an increase in the degree of mixing. The S index in Macedonian pine stands does not show a high dependence on the degree of mixing, but still there are small changes in structural parameters with its increase. The index determines the cause of the change in diversity well and determines equally well the competitive status of the relationships among the species. Serbian spruce in mixed stands can be characterized as extremely endangered by the competition of admixed tree species. The determined impact of mixing, considering the priority goal of the survival of this tree species (generally of species diversity), compromises the justification of establishing the strictest protection regimes of these stands, allowing them to develop spontaneously.

Correlation of the Investigated Parameters of Structural Diversity
There is not much literature data on the impact of mixing on the aggregation index. Contrary to the obtained results (Table 6), in stands with active management, the established ratio has a negative character of low intensity [62].
The pronounced link between the vertical fulfillment of the A rel profile and the mixing of stands is based on their mathematical correlation. The obtained values confirm the advantage of using the index in stands with different numbers of species [5], as well as the impact of deviation of structure from pure single-story stands on the growth of index values [78].
The impact of mixing on dimensional diversity indices is indirectly expressed through the impact on stand structure. Differences in the distributions of the investigated stands with increasing mixing had a positive effect on diameter differentiation indices (Td, GC). A study of different empirical and theoretical stand situations [46] found a strong impact of mixing on GC values. Significant differences in GC values were found between the pure and mixed even-aged stands of the most common species of Central Europe [79]. Mixed stands had on average significantly higher GC values (0.46) than pure stands (0.36), with an established correlation of medium strength.
The moderately strong positive correlation found in studies of the impact of mixing on the M and Td differentiation indices in virgin forests [70] and actively managed different forest types [62] was also confirmed by this study.
Prior research in stands of even-aged and almost even-aged structure found a strong correlation between Td and GC [80], which is also confirmed by this study. In contrast, the weak correlation of these indices in stands of more complex structure, virgin forests [70], diverse stands [80], and stands of different types [71], emphasizes the impact of structure on the change in their relationship. The reason for this should be sought in various mathematical assumptions of the above indices.
Since the hyperbolic tangent index S is a new index, no research has been performed on its correlation with other structural parameters. In addition to the tree characteristics diversity index, it has also been proposed as a useful tool for researching absolute and relative growth rates [48]. Considering the concept of using the index, the relationship between the S index and the mixing index M depends on the distribution of the observed characteristics of the admixed species in stands with different degrees of mixing. In this research, it has a weak positive correlation with the stand mixing index.

Conclusions
The structural diversity of the studied relict-endemic communities confirms their great importance in terms of biodiversity and reifies their need for protection and conservation in the Balkans.
Structural indices, whose values clearly vary in different stand situations, have proven to be useful and reliable indicators for defining stand characteristics, as well as for differentiating stands on the basis of them. Their mutual correlation was observed, as well as their correlation with the mixing of stands, which implies an increase in structural diversity with an increase in the degree of mixing in stands.
Compared to pure stands, mixed stands have a more pronounced skewness and kurtosis of the diameter structure, significantly wider variation widths of diameters and heights, a more pronounced vertical development, and higher values of structural indices, which all together result in a greater structural heterogeneity (diversity). These differences are significantly more pronounced in stands with a higher degree of mixing.
According to the obtained results, a 20% share of admixed species by volume and the spatial mixing of species M of 0.20, can be considered the lower limit of the impact of mixing on the change in structural characteristics of stands, with changes reflected in all parameters of horizontal and vertical structures.
The past development of the investigated stands without direct anthropogenic influence caused the formation of the initial phases of virgin stand forms. The established experiments contribute to a better understanding of structural relationships in the conditions of spontaneous stand development. The obtained low to medium values of structural diversity do not justify the existing management concept. In addition, they point to the need for an active management approach in order to preserve the biodiversity of rare and endangered forest communities. The dialectical principle of the development of forest communities, from the initial to the terminal phase, makes the establishment of permanent and strict protection regimes debatable from the aspect of biodiversity conservation. Finally, forest communities left to spontaneous development will certainly not go in the direction of conservation. Funding: This research was funded by the Ministry of Education, Science, and Technological Development of the Republic of Serbia within the project "Sustainable management of total forest potentials in the Republic of Serbia"-EVBR 37008.

Data Availability Statement:
The data presented in this study are available on request from the corresponding authors.