This study was carried out at the experimental farm of CSIR-Institute of Himalayan Bioresource Technology, Palampur (32°06′05″N; 76°34′10″E,1393 m altitude) India, during the cropping seasons of 2016 and 2017. The mean annual rainfall of the experimental region is approximately 250 cm, and the mean annual temperature is 18°C. The weather parameters, viz., maximum and minimum temperature, relative humidity, and rainfall, during the crop seasons are shown in Figure 1 . As per the USDA soil taxonomy classification system, the soil of the experimental region (Palampur) belongs to Alfisols (Sharma and Kumar, 2003). The physicochemical properties of the initial soil samples (0-15 cm depth) were studied before experimenting. The soil type of the experimental plots is silty clay, acidic in reaction (pH 5.52) having high organic carbon (1.49%) and medium available N (285.38 kg ha⁻¹), available P (14.56 kg ha⁻¹), and available K (256.59 kg ha⁻¹).

The agronomic traits recorded for all plots were plant height (cm), number of primary branches per plant (no. plant⁻¹), and fresh biomass yield (t ha⁻¹). The plant number in a unit area (no. m⁻²) has also been recorded. The crop was harvested at the end of September when all plots reached the full flowering stage. The full flowering stage is the indicator for harvesting (IHBT and Himgold, 2001). The whole plants (upper ground parts) are used for the extraction of essential oil. At harvest, the 1-m² area from the center of each plot was used for recording the agronomic traits and extraction of essential oil. Plant height was recorded by measuring the distance between the base of the stem and the top of the inflorescence from the five randomly selected plants from 1.0 m² quadrats. Likewise, the number of branches (no. plant⁻¹) was also recorded from the same plants. The minor variability of the population within the plots was observed in both the SM. In the case of broadcasting SM, fresh biomass yield was estimated by harvesting a 1.0-m² quadrat in each plot from three places. The area of 1 m from the plot boundaries was not used for sampling. In the case of the line sowing, the yield was estimated by harvesting three lines of 166.66 cm in length. Each line of 166.66 cm represents a 166.66 × 60 cm quadrat (1.0 m² quadrat). Finally, the fresh biomass yield (t ha⁻¹) was estimated from these 1.0 m² quadrats. For the extraction of essential oil and estimation of oil concentration in fresh biomass, 1-kg subsamples were taken from each plot.

The essential oil was extracted from the fresh aerial parts of the plants harvested from each plot. The fresh samples were placed in a 5-L volumetric flask with water in a 1:2 (w/v) ratio and then hydro-distilled in a Clevenger-type apparatus for 4 h. For essential oil extraction, 1 kg of fresh aerial biomass was used. After the completion of distillation, the quantity of essential oil was recorded for the further calculation of oil concentration in the fresh aerial biomass. The extracted oil was hydrated with anhydrous sodium sulfate and then stored in sealed glass vials at 4°C in a dark place for further analysis.

The essential oil extracted from all the treatments was further analyzed on a Shimadzu GC-2010 gas chromatograph (Shimadzu, Tokyo, Japan). The instrument was equipped with a flame ionization detector (FID). A DB-5 fused silica capillary column (30 m × 0.25 mm, film thickness 0.25 µm) was used with nitrogen as a carrier gas with a velocity of 1.05 ml min⁻¹. The oven temperature was programmed to hold at 70°C for 4 min and to increase to 220°C at the rate of 4°C min⁻¹ and finally held at 220°C for 5 min. The injector and detector temperatures were set at 240°C and 250°C, respectively. The quantification of the compounds was done based on the peak area normalization, and for each component, the response factor was considered equal to 1. Gas chromatography-mass spectrometry (GC-MS) analyses of all extracted oil samples were carried out by a Shimadzu QP2010 GC-MS system (Shimadzu, Tokyo, Japan), which was equipped with an autoinjector AOC-5000 and a DB-5 capillary column (30 m × 0.25 mm i.d. fused silica and film thickness 0.25 µm). Here, the carrier gas was helium, the flow rate was 1.1 ml min⁻¹, and the other operating conditions were the same as described in Section 2.5. For compound identification, the retention indices (RIs) for all volatile compounds on the homologous series of n-alkanes (C8–C24) were calculated. The volatile components of the essential oil of T. minuta were identified based on a comparison of their retention indices and mass spectra with the National Institute of Standards and Technology mass spectral (NIST-MS) library (Stein, 2005).

The data on agronomic traits, essential oil, and compositions of essential oil were subjected to analysis of variance (ANOVA) by Statistica 7 software (Stat Soft Inc., Tulsa, OK, USA). All the analyses were conducted based on the Statistica 7 software manual. For all the traits, replications were considered random effects, whereas SM and SR were considered fixed effects during ANOVA by Statistica 7 software. To understand the biomass yield response to SR under different SMs, the regression models were also developed between the independent variable as different SRs and the dependent variable as fresh biomass yield. The correlation matrix among the agronomic traits (number of plants, plant height, number of branches, and biomass yield) and essential oil yield of T. minuta was also developed. Principal component analysis (PCA), a multivariate analysis technique, was conducted on the chemical profiling of essential oil to assess the influences of treatment combinations and the nature of variations among the treatment combinations by the same statistical software. This multivariate analysis explains the relationship among the variables and the nature of the influences of the treatments on the variables. The degree of the relationship between any two variables is measured by the correlation coefficient value that is calculated by the cosine of the angle between their vectors in the PCA biplot (Yan and Rajcan, 2002; Dehghani et al., 2008). Simultaneously, heat maps were prepared based on the chemical profiling of essential oil by the XLSTAT software (version 2018.2). The respective F values and significant level (p-value) obtained from Statistica 7 are also presented in Table 1 to detect the significance of the main and interaction effects. In all parameters ( Tables 2 , 3 , Figure 2 ), the pairwise mean differences were assessed by the least significant difference (LSD) when the F test in ANOVA corresponding to treatments and interactions was significant at p ≤ 0.05.

The analyzed data revealed that the number of plants per unit area (no. m⁻²) was significantly (p ≤ 0.05) influenced by the SM and SR during both years ( Table 2 ). Irrespective of SR, broadcasting SM registered a significantly (p ≤ 0.05) higher number of plants per unit area (20.20 and 20.46 no. m⁻²) compared with the line SM (14.00 and 16.64 no. m⁻²) during both years. Among the SRs, the SR at 6 kg ha⁻¹ produced a significantly (p ≤ 0.05) higher number of plants (22.00 and 23.18 no. m⁻²) compared with the other SRs in 2016 but remained statistically at par with the SR at 5 kg ha⁻¹ in both years. Averaged across the SM, the SR at 6 kg ha⁻¹ produced approximately 2.1–2.2 times higher number of plants per unit area compared with the SR at 2 kg ha⁻¹, which produced the least numbers of plants (9.83 and 10.82 no. m⁻²) in both years. Broadcasting sowing plants registered significantly (p ≤ 0.05) higher plant height (172.36 cm) compared with the plants sown in line (161.77 cm) during 2017. Nevertheless, in contrast to plant height, line sowing plants produced a significantly (p ≤ 0.05) higher number of branches (10.52 no. plant⁻¹) compared with the plants sown by broadcasting during 2017. Averaged across the SM, the SR did not alter the plant height and number of branches during the first cropping season. In the second year, the effects of SR on these parameters were significant (p ≤ 0.05), and the maximum height (173.70 cm) and branches (12.33 no. plant⁻¹) were attained with higher SR (at 6 kg ha⁻¹) and lower SR (at 2 kg ha⁻¹), respectively.

Fresh biomass yield (Mg ha⁻¹) of T. minuta was significantly (p ≤ 0.05) affected by the SM and SR during both years ( Table 2 ), and the trend of the results was consistent over the duration of this study. Broadcasting SM produced significantly (p ≤ 0.05) higher biomass yield compared with line SM by approximately 15.8% and 7.6% during the first and second cropping seasons, respectively. Fresh biomass per plant decreased as SR increased up to 6 kg across the years. However, overall biomass yield (Mg ha⁻¹) was increased linearly with SR up to 4 kg and thereafter declined ( Figures 3A, B ). As SR increased from 2 to 6 kg, the individual plant yields declined by approximately 33% in both years ( Figures 3A, B ). Averaged across the SM, the overall biomass yield response to SR varied significantly (p ≤ 0.05) as SR increased from 2 to 6 kg ha⁻¹, and the utmost values (26.31 and 28.13 Mg ha⁻¹) were recorded with the SR at 4 kg ha⁻¹ in both years ( Figures 3A, B ; Table 2 ). As the level of SR increased from 2 to 4 kg, the total biomass yields increased by approximately 36%–56% ( Table 2 ). On the other hand, the higher SR (6 kg ha⁻¹) reduced biomass yield compared with the SR at 4 kg ha⁻¹. The interactions between SM and SR on all the agronomic traits were insignificant (p ≥ 0.05) in both years.

The correlation matrix was developed among the agronomic traits and essential oil yield to understand how yield attributes (plant density, plant height, and individual plant yield) contribute to the total biomass yield of T. minuta in this study ( Table 4 ). The number of branches is an important yield component. However, it was not included in the correlation matrix since a weak correlation was made with biomass yield. In this experiment, significant positive correlations (r = 0.90, p ≤ 0.01 in the first year; r = 0.64, p ≤ 0.05 in the second year) were found between the number of the plant (no. m⁻²) and total biomass yield (Mg ha⁻¹) in both years ( Table 4 ). There were also significant but negative correlations (r = −0.97, p ≤ 0.01 in the first year; r = −0.91, p ≤ 0.05 in the second year) between individual plant yield (g plant⁻¹) and the number of plants (no. m⁻²). Individual plant yield (g plant⁻¹) was significantly but negatively correlated with total biomass yield (r = −0.81, −80; p ≤ 0.01) and essential oil yield (r = −0.82, −0.79; p ≤ 0.01). However, significant and positive correlations (r = −0.99, p ≤ 0.01 in the first year; r = −0.93, p ≤ 0.01 in the second year) were documented between total biomass yield and essential oil yield. Thus, mixed correlations were demonstrated in this experiment. Overall, essential oil yield increased with increasing biomass yield, and total biomass yield increased with plant density.

The regression models have also established between the independent variable as different SRs and the dependent variable as fresh biomass yield ( Figures 3C, D ). The quadratic model (second-degree polynomial) best explained the biomass yield–SR relationship based on the statistical model selection criterion of R ² (R ^{2 =} 0.997 and 0.924 in the first year; R ^{2 =} 0.949 and 0.951 in the second year) values ( Figures 3C, D ). The biomass yield–SR relationship in the two SMs followed a similar trend in both years ( Figures 3C, D ) despite significant differences in actual yield levels between the two SMs. For the first year, the equations are

For the second year, the equations are

Where y and X indicate biomass yield (Mg ha⁻¹) and SR (kg ha⁻¹), respectively. The biomass yield response to SR varied significantly as the SR increased from 2 to 4 kg ha⁻¹ ( Table 2 , Figures 3C, D ). During the first cropping season, the average yield increased significantly, with a slope of 4.80 Mg ha⁻¹ for line sowing and 8.21 Mg ha⁻¹ for broadcasting, with an SR increase from 2 to 3 kg ha⁻¹ ( Figure 3C ). A moderate biomass yield increase, 3.86 Mg ha⁻¹ for line sowing and 2.05 Mg ha⁻¹ for broadcasting, was documented when SR was increased from 3 to 4 kg ha⁻¹, but with a less relative increase, 1.35 Mg ha⁻¹ for line sowing, when SR was increased from 4 to 5 kg ha⁻¹. In the case of the broadcasting sowing, yield declined by 1.70 Mg ha⁻¹. During the second year, a similar biomass yield increasing trend was recorded for the broadcasting SM with a slope of 6.91 and 1.43 Mg ha⁻¹ when the SR was increased from 2 to 3 kg ha⁻¹ and 2 to 3 kg ha⁻¹, respectively ( Figure 3D ). For the higher SR, when SR was increased from 4 to 6 kg ha⁻¹, yield declined by 2.21 Mg ha⁻¹ for line sowing and 3.04 Mg ha⁻¹ for broadcasting.

The SM and SR did not alter the concentration (%) of essential oil in the aboveground fresh biomass of T. minuta in the present study, and the trend was consistent over the duration of this study ( Figures 2A, B ). However, higher SR registered a marginally higher concentration of essential oil (0.24% in 2016 and 0.33% in 2017) compared with low SR, and the overall concentration of essential oil in biomass was slightly higher in 2017( Figure 2B ). The total yield of essential oil was significantly (p ≤ 0.05) influenced by SM in the first cropping season, and broadcasting SM produced approximately 23.4% and 4.4% higher essential oil compared with line SM during 2016 and 2017, respectively ( Figure 2C ). The effect of SR on essential oil yield was significant (p ≤ 0.05) in both years, and the utmost values (63.52 and 93.04 L ha⁻¹) were recorded with SR of 4 kg ha⁻¹ in both years ( Figure 2D ). Intermediate SR (at 4 kg ha⁻¹) produced approximately 56.9%–61.3% and 7.2%–12.3% higher oil yield compared with low SR (at 2 kg ha⁻¹) and high SR (at 6 kg ha⁻¹), respectively.

Qualitative evaluation of the essential oils of T. minuta under different SMs and SRs was done through the identification and quantification of the compounds ( Table 3 ). In Table 3 , we presented only 13 compounds, which were identified in all the treatments and replications. The compounds, which were not identified in all the treatments and replications, were not presented in this table. These 13 compounds contributed approximately 88.66%–96.01% of the total volume of essential oil samples obtained in this study. The data revealed that only a few compounds (E-caryophyllene, α-humulene, Z-β-ocimene, E-β-ocimene, and tagetone in the second year only) were affected by the SM in this study. Broadcasting SM produced a significantly (p ≤ 0.05) higher amount of E-caryophyllene and α-humulene in both years ( Table 3 ). In T. minuta essential oil, one of the major and important compounds is Z-β-ocimene, which was significantly (p ≤ 0.05) influenced by the SM in the first cropping season, and maximum Z-β-ocimene content (35.69%) was recorded with broadcasting SM. However, the trend was not consistent between the years ( Table 3 ). In the case of E-β-ocimene, the trend over the years was similar, and a significantly (p ≤ 0.05) higher value (35.69%) was recorded with broadcasting SM during 2017. The SR did not influence the composition of essential oil in this study; however, the maximum (36.60% and 35.38%) values of Z-β-ocimene were documented with the application of higher SR (at 6 kg ha⁻¹) ( Table 3 ).

Variations in the chemical profile of the essential oil under the interaction effects between SM and SR have also been illustrated by heat maps ( Figures 4A, B ). A total of 13 compounds were used in the heat maps, which depict the changes in the chemical profiling of the essential oil of T. minuta due to the interaction between the SM and SR. Clustering was done for the chemical profile of essential oil as well as for treatment combinations by developing dendrograms ( Figures 4A, B ). The heat map exhibited that the accumulation patterns of most of the compounds were uneven over the years ( Figures 4A, B ). During the 2016 season, maximum (38.54%) and minimum (27.92%) Z-β-ocimene concentrations were recorded with M_BS₅ (broadcasting SM with SR at 6 kg ha⁻¹) and M_LS₄ (line SM with SR at 5 kg ha⁻¹), respectively. However, in 2017, M_LS₅ (line SM with SR at 6 kg ha⁻¹) and M_BS₃ (broadcasting SM with SR at 4 kg ha⁻¹) registered maximum (35.50%) and minimum (25.22%) quantities of Z-β-ocimene, respectively ( Figures 4A, B ). Dihydrotagetone is another major compound that also shows inconsistent accumulation patterns under different treatment combinations over the years. In 2016, the maximum value (29.14%) was recorded with M_LS₄ (line SM with SR at 5 kg ha⁻¹), whereas in 2017, the maximum value (33.81) was recorded with M_BS₂ (broadcasting SM with SR at 3 kg ha⁻¹). Moreover, the heat map also indicated that treatment combinations were divided into two major clusters in both years ( Figures 4A, B ).

PCA analyses were performed taking into account the identified compounds in essential oil under different SMs and SRs for both years ( Figures 5A–D ). In this study, the PCA biplot exhibited that the first component (PC1) and the second component (PC2) jointly explained 72.99% and 72.02% of the total variations for 2016 and 2017, respectively ( Figure 5 ). The PCA biplot technique confirmed the relationship among the variables (compounds) and the nature of influences of the treatment combinations on the variables. In the PCA biplot, the vectors from the origin to each marker of the compounds enable visualization of the relationships between and among the compounds ( Figures 5A, C ). The strength of the relationship between any two variables is assessed by the correlation coefficient, which is estimated by the cosine of the angle between their vectors in the PCA biplot (Yan and Rajcan, 2002; Dehghani et al., 2008). Two variables are independent, positively correlated, and negatively correlated when the angle between their vectors is 90°, less than 90°, and greater than 90°, respectively. Similarly, variables (compounds in this study) with longer and shorter vectors have strong and weak correlations with the PCs, respectively. In 2016 ( Figure 5A ), strong positive associations between dihydrotagetone and tagetone; among Z-β-ocimene, E-β-ocimene, and alloocimene; between Z-ocimenone and E-ocimenone; and between E-caryophyllene and α-humulene were observed as exhibited by the narrow acute angles between their vectors (e.g., r = cos 0 = +1). Also, Figure 5A demonstrates the negative association between dihydrotagetone and E-caryophyllene, between dihydrotagetone and α-humulene, and between E-ocimenone and E-β;-ocimene, since large obtuse angles between their vectors were observed. Similarly, in 2016, strong positive associations between E-β-ocimene and limonene, between sabinene and dihydrotagetone, between Z-ocimenone and tagetone, and between E-caryophyllene and α-humulene were observed as exhibited by the narrow acute angles between their vectors ( Figure 5C ). Furthermore, Figure 5C demonstrates the negative association between E-caryophyllene and E-β-ocimene, between E-ß-ocimene and tagetone, and between E-ocimenone and E-β-ocimene, since large obtuse angles between their vectors were observed. PCA analysis based on 13 volatile compounds revealed a clear separation into four and three treatment combination groups for 2016 and 2017, respectively ( Figures 5B, D ). In 2016, line SM with SR from 2 to 5 kg ha⁻¹ developed a single group and was placed in the positive end of PC1 ( Figure 5B ). Moreover, broadcasting SM with SR of 3 and 6 kg ha⁻¹ developed an independent group and was separated by the PC2 ( Figure 5B ). In 2017, the PCA biplot ( Figure 5B ) showed that the broadcasting SM with SR 2, 3, and 4 kg ha⁻¹ came under the same group in terms of essential oil composition; however, all three treatment combinations were placed in the negative end of PC1.