Philippa Kaur

Synergy between S. cerevisiae and R. etli in biofilm formationWhen S. cerevisiae Mat α Σ1278h and R.etli CE3 were grown in minimal medium with low glucose concentrations (0.1%), these species adhered to abiotic surfaces to form biofilms (Fig. 1). Interestingly, R. etli and S. cerevisiae formed a mixed biofilm whose biomass was ~ 3 times greater than that of either single-species biofilm (Fig. 1a). In addition, at 24 h, the number of colony-forming units (CFU)/cm2 of R. etli CE3 in the mixed biofilm was higher than that in the pure biofilm (Supplementary Fig. 1). Confocal laser scanning microscopy of biofilms stained with the Live/Dead Kit (propidium iodide and SYTO9) showed that in the mixed biofilm, the yeast cells formed patches, and the bacterial cells covered most of the surface (Fig. 1b). In contrast, monospecies biofilms of R. etli and S. cerevisiae had lower structural complexity and contained a greater (80%) number of dead cells, and their individual densities were lower than their populations in the mixed biofilm (Fig. 1b). These results suggest that in mixed biofilms, S. cerevisiae promotes bacterial growth.Figure 1The interaction between S. cerevisiae and Rhizobium etli CE3 results in the formation of a structurally complex and more productive biofilm in terms of biomass. (a) Biofilm formation of R. etli CE3 and S. cerevisiae Σ1278h Mat α and biofilm growth over time in minimal dextrose medium. The data are representative of 3 independent experiments +/− the S.D. values. (b) Top view and cross section of confocal micrographs of the S. cerevisiae-R. etli mixed biofilm and the single-species biofilms. Magnification 40 × . The images are representative of 3 independent experiments. Biofilms were developed on glass microscope slides and stained with a LIVE/DEAD viability kit. Red fluorescence indicates dead cells, and live cells are colored green. Images were acquired 24 h after inoculation.Full size image
S. cerevisiae secretes dicarboxylic acids that promote R. etli growth and biofilm formationWe found that the R. etli colonies that grew close to S. cerevisiae on solid glucose minimal medium were larger than those growing far from yeast colony (Fig. 2).Figure 2Yeast cells produce dicarboxylic acids that promote the growth of R. etli. (a) R. etli growth in coculture with S. cerevisiae BY4741 mutants (aco1Δ, fum1Δ, sdh1Δ and mdh1Δ) that accumulate dicarboxylic acids and a BY4741 strain with blockade of the aerobic respiratory chain (rho-). (b) Test on solid medium showing that S. cerevisiae BY4741 (*) secretes compounds that promote bacterial growth (  >). In contrast, BY4741 rho- cells (ρ), which do not produce dicarboxylic acids, do not promote the growth of R. etli CE3. R. etli CE3 cells were spread over MMD agar, and yeast cells were spotted in the center. (c) Top view of light micrographs of dual-species biofilms; S. cerevisiae (arrowhead) and R. etli (arrow). Biofilms were developed on glass microscope slides and stained with crystal violet. Magnification 20 × . The images are representative of 3 independent experiments. (d) Growth of R. etli strains in coculture with S. cerevisiae BY4741. The growth of the rhizobium strains was estimated at 24 h. R. etli CE3 strains: wild-type (wt), dctA- containing an empty expression plasmid (dctA-) and dctA- containing a plasmid expressing dctA (dctA-/dctA). The data are representative of 3 independent experiments +/− the S.D. values.Full size imageWe used a visual growth promotion assay on solid medium to screen for S. cerevisiae knockout strains (YKO library) that influenced bacterial growth. 159 yeast mutants were unable to promote R. etli CE3 growth (Supplementary Table 3). In general, these mutants were defective in mitochondrial function. Interestingly, we found that 5 strains with mutations in genes coding for enzymes involved in the TCA cycle showed an enhanced ability to promote bacterial growth compared to that of the wild-type strain (Fig. 2a).To determine how the S. cerevisiae mutants may affect the fungal-bacterial interaction, we analyzed factors that may be altered in mutants with mitochondrial function defects and a compromised TCA cycle.We compared the production of TCA intermediates between the wild-type and mutant yeast strains. Mutants defective in mitochondrial function (mef1Δ, gep5Δ, sdh2Δ, ppa2Δ, imp1Δ, cox7Δ, cyc1Δ and cyc2Δ) produced low amounts of tricarboxylic acids (Supplementary Fig. 2a). In contrast, the aconitase mutant (aco1Δ) produced 60% more citrate and succinate; the fumarase mutant (fum1Δ) resulted in fumarate accumulation; the succinate dehydrogenase mutants (sdh1Δ and sdh4Δ) produced 80% more succinate; and the mitochondrial malate dehydrogenase mutant (mdh1Δ) produced 60% more malate and succinate (Supplementary Fig. 2b). These results suggested that the large quantities of tricarboxylic acids secreted by the mutant yeast played a role in promoting bacterial growth in the cocultures.We analyzed the biomass of mixed biofilms formed by yeast cells defective in mitochondrial function (Σ1278B petit mutant). The ability of the wild-type and the petit mutant strains to form a monospecies biofilm was similar (Supplementary Fig. 3). In contrast, the mixed biofilm formed by yeast cells defective in mitochondrial function was significantly lower in biomass than that formed by the wild-type yeast strain (Fig. 2c). Also, Σ1278B petit mutant produced low amounts of tricarboxylic acids (Supplementary Fig. 2a).We next measured the biomass of the mixed biofilm formed by S. cerevisiae and a Rhizobium mutant unable to take up C4-dicarboxylic acids (dctA-). This evaluation revealed that C4-dicarboxylate uptake by R. etli is necessary to form mixed biofilms with high biomass (Fig. 2d).A symbiotic plasmid is involved in the phenotypic plasticity of R. etli.
The genome of Rhizobium etli CE3 is composed of a chromosome and 6 plasmids (pA, pB, pC, pD, PE and pF)11. To determine whether elements encoded by these replicons can participate in the establishment of commensalism, we evaluated the formation of biofilms by yeast and R. etli strains lacking these replicons12. We found that lack of pA, pB, pC or pF did not affect the ability of bacteria to coexist with yeast (Fig. 3a). Interestingly, a strain cured of plasmids pA-/pD- could not coexist with S. cerevisiae to form a mixed biofilm and obtain the benefits provided by the fungus (Fig. 3a).Figure 3Plasmids pA and pD encode proteins performing functions that are necessary for the coexistence of bacterial cells with yeast. Growth of R. etli strains in biofilms with S. cerevisiae S1278B. (a) Growth in mixed biofilms of R. etli strains lacking the plasmids; pA, pB, pC, pF and in one case of two plasmids, pA-/pD-. The growth of the rhizobia strains was assessed at 24 h. (b) Scheme of the genes contained in a cosmid that partially complements the growth of the pA-/pD- strain in mixed biofilms. Here, 3, 2 and only one gene was amplified to generate the plasmids AD1, AD2 and AD3, respectively, as indicated in the figure. (c) Growth of R. etli strains in mixed biofilms. Strains AD1 and AD2 are R. etli pA-/pD- cells that carried plasmids AD1 and AD2, respectively. The growth of rhizobium strains in mixed biofilms was estimated at 24 h. The data are representative of 3 independent experiments +/− the S.D. values.Full size imageTo determine the genetic elements from the symbiotic plasmid involved in the interaction with yeast, we complemented the R. etli pA-/pD- strain with a cosmid library containing fragments of partial digestion (EcoRI) of the R. etli CE3 genome13. We found that a cosmid containing 9 ORFs from plasmid pD (GenBank: U80928.5) partially restored the ability of R. etli pA-/pD- to form a mixed biofilm (Fig. 3b). This cosmid contains 7 insertion sequences (IS) and a predicted operon encoding a probable peptide pheromone/bacteriocin exporter (RHE_PD00332) and a probable bacteriocin/lantibiotic ABC transporter (RHE_PD00333) (Fig. 3b).The complete operon or only the ABC transporter gene, including its endogenous promoter and terminator regions, was cloned into plasmid pBBR1MCS-3, and the resultant plasmids were named AD1, AD2 and AD32, respectively (Supplementary table 1 and 2). We found that complementation with the complete operon (plasmid AD2) partially restored the ability of R. etli pA-/pD- to form a mixed biofilm with yeast (Fig. 3c). In contrast, complementing with the RHE_PD00332 gene (plasmid AD3) does not restore the phenotype. It is necessary to complement only with the RHE_PD00333 gene to determine if its product is involved in the phenotypic plasticity of R. etli. These results suggest that the ABC transporter gene (RHE_PD00333) is involved in the fungal-bacteria interaction.
S. cerevisiae produces a small molecule that affects R. etli growthTo determine how S. cerevisiae affects the growth of R. etli pA-/pD- (Fig. 4a), we evaluated the inhibitory activity of methanol extracts of S. cerevisiae culture supernatants.Figure 4S. cerevisiae s1278B produces a small molecule that only affects the growth of R. etli strains that do not harbor the symbiotic plasmid and plasmid A. (a) S. cerevisiae and R. etli strains were inoculated in close proximity onto MMD soft agar. R. etli pA-/pD- grew, forming a swarm far from the yeast colony. (b) Inhibition of R. etli pA-/pD- growth by 5 µg/mL of a purified compound from the yeast supernatant, which we named Sc2A. (c) Proposed molecular structure of Sc2A.Full size imageInterestingly, we found that the methanol extract inhibited R. etli pA-/pD- growth but had no activity against wild-type R. etli (Fig. 4b). We investigated the chemical constituents of the S. cerevisiae culture supernatants. After succesive organic solvent extractions, the methanolic extract was fractionated by HPLC and 8 fractions were obtained. Each fraction was tested for its determine its effect on the growth of R. etli pA-/pD-. Only a fraction with the ability to inhibit the growth of R. etli pA-/pD- was identified. This resulted in ~ 90% pure sophoroside, judging by its appearance as a dominant peak in the mass spectra obtained by Fast Atom Bombardment Mass Spectroscopy (FAB). As a result, a new sophoroside with bacteriostatic activity, named Sc2A, was isolated (Fig. 4c). The structure of Sc2A was elucidated by a combination of extensive spectroscopic analyses, including 2D NMR and HR-MS.Sc2A was isolated as a crystalline powder with a positive optical rotation ([α]D25 + 13.7°, c0.58, H2O). The molecular formula of Sc2A was determined to be C30H50O24 from its positive-mode FAB data (m/z 794.26 [M + H]+), which was consistent with the 13C NMR data. RMN1H (CD3OD, 400 MHz) data for Sc2A: δ 5.1 d (J = 3.6 Hz), 4.4 d (J = 8 Hz), 4.23 dd (J = 9, 4.8 Hz), 3.79 t (J = 10.8, 14.4 Hz), 3.73 m, 3.67 m, 3.639 m, 3.63 dd (J = 8, 9.2 Hz), 3.53 dd (J = 5.6, 5.2 Hz), 3.36 dd (J = 3.6, 4 Hz), 3.31 dd (J = 8, 8 Hz), 3.10 dd (J = 8, 7.6 Hz), 2.77 dd (J = 4.4, 6.8 Hz), 2.61 m, 2.46 m, 2.33 m, 2.12 m. RMN13C-DEPT (CD3OD, 400 MHz) data for Sc2A: δ 181.2 (C), 175.9 (C), 98.1(CH), 93.8 (CH), 78.05 (CH), 78.02 (CH), 76.30 (CH), 74.92 (CH), 73.80 (CH), 73.11 (CH), 71.78 (CH), 71.72 (CH), 64.37 (CH2), 62.87 (CH2),62.72 (CH2), 57.24 (CH), 30.70 (CH2), 26.19 (CH2), 28.21 (CH2).The IR spectrum of Sc2A displayed characteristic absorptions of 3416.34 cm-1 (O–H), 1642.10 (C = O), 1405.44 (C–OH), 1242.93 (C–O–C), 1040.36 (C-H), and 598.48 (O-C-O).Sc2A possesses a sophorose linked by 2,5 hexanedione to another molecule of sophorose (Fig. 4c).Sc2A induces the expression of genes involved in symbiosisExpression from the nifH and fixA promoters was studied in R. etli monocultures and cocultures with yeast by monitoring GUS activity in living cells. Cells were grown on solid PY-D medium for 1 day, and monitoring of GUS expression showed that the nifH promoter was strongly induced when R. etli was grown with yeast in liquid medium and on solid medium (Fig. 5).Figure 5The expression of Rhizobium etli genes involved in symbiosis is induced in cocultures with yeast or by exposure to the small molecule Sc2A. (a) Activity of different R. etli promoters in monoculture (Re) or in coculture with yeast (+ Sc). Cells were cultured for 24 h in 1 ml of PY-D in 1.5-mL tubes. The tubes were kept closed to generate an environment with a low oxygen concentration. (b) Activity of the nifH promoter in R. etli cells grown alone (Re) or in coculture with yeast (+ Sc) on PY-D agar. (c) Effect of Sc2A on the expression of the nodA gene in R. etli cells grown in liquid culture. Cells stimulated with the flavonoid naringenin were included as a positive induction control. The data are representative of 3 independent experiments +/− the S.D. values.Full size imageAt the beginning of the symbiosis, the legume roots exude flavonoids, which induces in R. etli the expression of a group of genes (nod) involved in the synthesis of lipochitooligosaccharides, also called nodulation factors (NFs). Recognition of NFs by the host plant triggers both rhizobial infection and initiation of nodule organogenesis14. NodA protein is involved in N-acylation of the chitooligosaccharide backbone of NFs. Given the participation of nodA in the interaction of R. etli with a eukaryote, we decided to evaluate the expression of this gene in response to exposure to 5 µg/mL of Sc2A (this concentration is similar to that found in cocultures). We found that Sc2A induces the expression of nodA (Fig. 5c). However, the levels of induction of nodA were moderated compared to the values obtained upon naringenin induction (Fig. 5c). More

Morphology and light-emitting behavior of the undescribed Japanese Terebellidae wormIn 2016, some of the present authors were exploring shallow coastal waters (depth less than 1 m) to observe the ecological behaviors of marine animals in the Noto Peninsula, when they discovered unknown violet-light-emitting worms. At the sampling point, the worms were living in small holes (a few centimeters in diameter) or in cracks in rocks covered by sand at the shallow sea bottom (Supplementary Fig. S1). We successfully video-recorded their emission of violet light from the whole tentacle stretching into sea water when stimulated by air bubbling at night (Fig. 1A–C; Supplementary Videos S1 and S2). The violet-light emission consisted of rapid flashes with variable duration in the order of milliseconds (Supplementary Video S3), as observed for the worm P. perplexus in response to stimulation17. From our morphological observation, we identified the violet-light-emitting worm as a member of Polycirrus on the basis of the following characteristics18: (1) a sheetlike prostomium covering the upper lip; (2) avicular unicini on some neuropodia; (3) no branchiae. The specimens also have the following characteristics: (1) neurochaetae beginning on last notochaetigerous segment, chaetiger 14; (2) uncini with a long neck and concave base; (3) notopodial pre- and post-chaetal lobes both similar shape. These characters are also found in Polycirrus disjunctus Hutchings and Glasby18; however some of the characters in parapodial lobes and chaetae have differentiation. Thus, we concluded that this species should be treated as an undescribed species. Further comparative observation is needed to describe the species. At this time we treated the Polycirrus species observed in this study as Polycirrus sp. ISK. Application of an electric pulse also caused clear light emission from the tentacles of the living worm (Fig. 1D; Supplementary Video S3), and the luminescence spectrum showed that its λmax was 444 nm or slightly longer, depending on the individual (Fig. 2A). We also found that light emission was efficiently induced by the addition of KCl solution and observed the time course of light emission with rapid fluctuations with variable duration in the order of milliseconds for up to 30 s (Fig. 2B). The flash pattern was similar to that observed in a study of P. perplexus17. In the genus Polycirrus, P. medius and P. nervosus in Japan have been described18,19. However, the morphological features of the species in the present study differed from these species on the basis of our observations described above.Figure 1Photographs of Polycirrus sp. ISK. (A) Polycirrus sp. ISK in its natural habitat with bright-field illumination. (B) Bioluminescence of Polycirrus sp. ISK in its natural habitat without bright-field illumination. The worms were stimulated by air bubbling from SCUBA gear. (C) A single worm with stretched tentacles. Tentacles are indicated by white arrows. (D) The worm with light emission at the tentacles. This worm was stimulated by an electric shock. Scale bars = 100 mm for A and B, 10 mm for C and D. Each photograph was extracted from the videos recorded with the following settings: sensitivity, ISO 51200 or 11 lx; white balance, 4300 K or 5800 K; shutter speed, 1/30 s or 1/60 s; iris, F1.8-3.5; frame rate, 29.97 fps or 60 i; frame size, 1920 × 1080 pixels. Original high quality videos are available at https://youtu.be/KEsU0kWAEfg and https://youtu.be/24dxvPlBDB0Full size imageFigure 2Luminescence spectra and KCl-induced light emission of Polycirrus sp. ISK. (A) Spectrum analysis of Polycirrus sp. ISK using a living worm stimulated by an electric shock. The luminescence spectra were obtained from two different individuals. The λmax represented in closed circles and open circles were 444 nm and 446 nm, respectively. (B) Typical light-emission signal of a living worm soaked in 667 mM KCl. The black line indicates luminescence intensity after adding KCl solution, and the gray line indicates luminescence intensity before adding KCl solution.Full size imageJapanese Polycirrus spp. have not been described as luminous worms according to our review of the literature and web pages. In addition, the number of reports for new Polycirrus spp. from all over the world has been increasing, but a limited number of species are known to emit light13,17,18. Our finding of KCl-induced light emission from Polycirrus sp. ISK suggested that we can easily test the light-emitting ability of Polycirrus spp. by luminescence measurement just after adding KCl solution. A spectrum pattern has been reported for only one species, P. perplexus collected in California17, and it would be necessary for further understanding of these species to examine the light-emitting abilities and to compare light-emitting behaviors and spectrum patterns. The color of bioluminescence is often related to habitat, and light in the blue range is typical for pelagic species20. Thus, one of the points to be focused on is the ecological function of the violet-light emission of this worm inhabiting in a shallow coastal water environment. In P. perplexus, deterring predation is a possible function of luminescence based on that species’ habitat and its violet-light emission17,21. As shown in Supplementary Videos S1, S2, which are the first video records of in situ light emission of a Polycirrus species, the air bubble-stimulated luminescence of Polycirrus sp. ISK in its natural habitat also seemed to deter predation, but this explanation is still speculative.Differentially expressing genes between the tentacles and the rest of bodyA few years after discovering this worm, we found it difficult to collect enough of them to conduct common biochemical and chemical analyses because we did not find a place densely inhabited by hundreds of the worms whose wet weight was a few tens of milligrams (e.g. 16.5, 29.8, or 31.8 mg). Next, we conducted RNA-Seq analysis. In luminous animals with strong light emission, such as firefly or syllid polychaetes (Syllidae), luciferase expression is high especially at the luminous organ or in the whole body22,23. On the other hand, the light emission of Polycirrus sp. ISK was not so strong compared to that of fireflies, and the light-emitting area was limited to the tentacles. In addition, the genetic information related to the tentacles responsible for various ecological functions is still limited. Thus, in the present study we decided to purify RNA from the tentacles and the rest of body separately (Fig. 1C) and performed RNA-Seq analysis followed by a computational analysis using the MASER pipeline24. By de novo assembly, 110,775 contigs were predicted; 26.1% of them showed more than twice the expression level in the tentacles than in the rest of body, whereas 20.8% showed more than twice the expression in the rest of body than in the tentacles. When we performed a blastX search to the NCBI nr database for the contigs longer than 300 bp, 35.6% showed significant homology with registered genes with e-values of less than 1e−10. The average length for these contigs was 1384 bp, and half of them were in the range of 463–1863 bp (Supplementary Fig. S2). In the assembled sequence, we found the cytochrome oxidase subunit I (COI) gene and tried to construct a phylogenetic tree. However, the obtained phylogenetic tree was unreliable due to the low bootstrap values as shown in Supplementary Fig. S3.To focus on the tissue-specific genes, we first picked up genes with high expression levels based on high fpkm values (over 1000) and then ranked these genes based on the tissue-specificity judged by the comparison of fpkm values in tentacles and the rest of body. In tentacle-specific genes, we found that some genes coding for lectin(-like) domains were ranked in the top eight as shown in Supplementary Table S1. Of the top eight genes in the rest of body-specific genes (Supplementary Table S2), seven exhibited no similarity to any genes, and the remaining gene exhibited significant similarity to a hypothetical protein of Capitella teleta, which is a Polychaetes species with whole-genome information available25. Recently, TPM is preferably used to normalize expression level, and the value is used for statistical differential expression analysis26, and we also calculated TPM for tissue-specific genes (Supplementary Table S3).As we were unable to conduct statistical differential expression analysis due to no biological/technical replication resulted from difficulties in the sample collection, we simply compared TPM value between the tentacle and the rest of body samples. The ratio of TPM (tentacle/rest of body) was calculated, and then top 100 genes (Fig. 3A), which were highly expressing in the tentacle, were selected. Similarly, top 100 genes highly expressing in the rest of body were selected using the ratio of TPM (rest of body/tentacle) (Fig. 3B). These gene lists were annotated by gene ontology (GO) terms and analyzed using WEGO program27. WEGO results showed different expression patterns for the tentacle and the rest of body. In the tentacle, GO terms including cell adhesion, biological adhesion, small molecular binding, positive regulation of biological process, regulation of response to stimulus, carbohydrate binding, and immune response were significantly higher (Fig. 3C, D). In the rest of body, GO terms including hydrolase activity, catalytic activity, localization, and establishment of localization were significantly higher. In the top 100 genes highly expressing in the tentacle, we found 21 genes annotated as a gene coding for fucolectin by blast search (Supplementary Table S4). Fucolectin is a fucose-binding lectin involved in the innate immunity of diverse invertebrate species28. However, its function in invertebrates remains unclear, and no information is available for Terebellidae, including sequence information. Fucolectin was first identified in eel with mRNA distribution mainly in liver and gill28. In sea cucumber, expression of the fucolectin gene is confirmed in respiratory trees, muscle, and tentacle29. We were not able to see whether this gene was expressed in the respiratory organ of Polycirrus sp. ISK because a characteristic of the genus Polycirrus is the absence of branchiae18. Nevertheless, the tentacle-specific expression of fucolectin was consistent with the observation in sea cucumber, and the high expression of such proteins involved in innate immunity seemed reasonable because tentacles stretching out of their bodies can be damaged by attack of predators and thus are threatened by infectious bacteria and other pathogens11, as is the respiratory organ. In addition, localization of antimicrobial compounds in Terebellidae worms is suggested to be of antiseptic importance in damage by predation14. This study would provide indispensable information about the ecological meaning of Polycirrus sp. ISK’s life in future genetic studies.Figure 3WEGO analysis of highly expressing genes in the tentacle and the rest of body. (A) Box plot graph for the distribution of TPM value for top 100 genes highly expressing in the tentacle. Corresponding genes in each part are colored in the same gradation color according to the TPM value (red to blue form higher to lower value). (B) Box plot graph for the distribution of TPM value for top 100 genes highly expressing in the rest of body. Each gene is colored as in (A). (C) WEGO analysis of top 100 genes highly expressing in the tentacle (orange bar) and the rest of the body (blue bar). (D) P-values from Chi-square tests obtained by WEGO analysis. CC cellular component, MF molecular function, BP biological process.Full size imageTranscripts coding for luciferase-like genes in the wormTo find genes similar to the known luciferase, which is an enzyme oxidizing a specific compound called luciferin to emit light, from related species in polychaetes, we performed a blastX analysis against the Odontosyllis luciferase sequence using our RNA-Seq data. We found a gene coding for a protein that exhibited similarity to Odontosyllis luciferase, but the e-value was more than 1e−10 (Supplementary Fig. S4). In addition, the top hit for this gene analyzed by blastX was annotated to code an uncharacterized protein of Saccoglossus kowalevskii (Hemichordata), and its specific function was not predicted. Other hits were for genes from Chordata, Mollusca, and other phyla but there was no hit from Annelida. This result would suggest that the light-emission system of Polycirrus sp. ISK differs from that of the genus Odontosyllis, although further experiments using high purity Odontosyllis luciferase and the substrate will be necessary to confirm this. In further blastX analyses of representative luciferases, photoproteins, and a putative luciferase [luciferases from the ostracod Cypridina noctiluca (Accession number: BAD08210.1), the copepod Gaussia princeps (AAG54095.1), the deep-sea shrimp Oplophorus gracilirostris (BAB13775.1 and BAB13776.1), the firefly Photinus pyralis (AAA29795.1), the sea pansy Renilla reniformis (AAA29804.1); photoproteins from the hydrozoan jellyfish Aequorea victoria (AAA27720.1), the hydroid Clytia gregaria (CAA49754.1), the hydroid Obelia geniculate (AAL86372.1); a putative luciferase from the tunicate Pyrosoma atlanticum30 sequences using our RNA-Seq data], we found some tissue-nonspecific genes whose sequences exhibited similarity to firefly luciferase (FLuc) or Renilla luciferase-like protein (RLuc-like) sequences with an e-values of less than 1e−10 and percent identity of more than 50%. FLuc is a member of the acyl-adenylate-forming superfamily of enzymes responsible for firefly luciferin-dependent bioluminescence, which is found in terrestrial luminous beetles emitting light ranging from green to red31. Previously, a putative acyl-CoA synthetase protein was found in the luminous organ of firefly squid emitting blue light32, but there is no clear biochemical evidence that such protein is responsible for firefly squid’s bioluminescence. On the other hand, RLuc is responsible for coelenterazine-dependent bioluminescence, which is found in marine luminous organisms belonging to various taxa. An RLuc-like protein is found to be localized in luminous organs of the brittle star Amphiura filiformis, as revealed by taking advantage of the cross reactivity of anti-RLuc antibody to A. filiformis RLuc-like protein33. A recent study reported that recombinant RLuc-like protein found in P. atlanticum exhibited luciferase activity to coelenterazine30. However, an RLuc-like protein from sea urchin Strongylocentroutus purpuratus is confirmed to exhibit dehalogenase activity to various substrates but no luciferase activity to coelenterazine34. Therefore, it is suspected that Polycirrus sp. ISK possesses a luminescence system using an RLuc-like enzyme.Coelenterazine content in the wormTo investigate whether Polycirrus sp. ISK possesses not only a Renilla luciferase homologous gene but also coelenterazine, we analyzed an ethanolic extract of Polycirrus sp. ISK by UPLC with a UV–visible detector (Fig. 4). The obtained UPLC chromatogram did not show a peak corresponding to that of authentic coelenterazine. When further checking the chromatogram, we found the peak at a retention time similar to those of authentic coelenteramide and coelenteramine, which can be formed from coelenterazine. However, the absorption spectrum obtained by UPLC analysis and the mass spectrum obtained by MS/MS analysis were not identical to those of authentic coelenteramide or coelenteramine (Fig. 4 and Supplementary Figs. S5 and S6). In addition, when the worm extract was mixed with a recombinant RLuc, we did not detect luminescence using a luminometer. These results suggested that the luminescence system in the worm was independent of coelenterazine, although a RLuc homologous gene was found. Similarly, the existence of an RLuc homologous gene was reported in P. atlanticum, which has been suggested to use a coelenterazine-independent luminescence system relying on bacterial bioluminescent symbionts30,35. We also mixed the worm extract with a recombinant cypridinid luciferase, but we did not detect luminescence using a luminometer. This result was consistent with Harvey’s observation for P. caliendrum16. To examine whether the luminescence system is based on luciferin–luciferase reaction, which is found in various luminous animals including some syllid Odontosyllis spp.23,36,37,38,39, we prepared two different extracts of the worm using 100 mM HEPES–NaOH buffer (pH 7.4) and methanol, and subsequently subjected a mixture of the two to luminescent measurement. As a result, no light emission was detected from the mixture of the buffer and methanolic extracts of the worm. This result was also consistent with Harvey’s observation for P. caliendrum16. However, there is still a possibility that the light emission is based on luciferin–luciferase reaction, because luciferin–luciferase reaction found in fireflies or luminous mushrooms requires a cofactor such as ATP or NADPH, and we did not test all possible conditions due to the limitation of the number of collected specimens. In addition, extraction of luciferin and luciferase in the active form is sometimes difficult, as shown in previous studies37. Further studies using hundreds or more of the specimens must be performed to elucidate the mechanism underlying the violet-light emission.Figure 4Comparison of the ethanolic extract of Polycirrus sp. ISK with CTZ, CTMD, and CTM. (A) UPLC analysis of (a) the extract, (b) authentic CTZ, (c) authentic CTMD, and (d) authentic CTM using a multiwavelength detector. The black solid line indicates detection at 333 nm, and the blue solid line indicates detection at 435 nm. The compound between the red vertical dashed lines was collected for MS/MS analysis. (B) Absorption spectra of the compound from the extract, CTZ, CTMD, and CTM obtained at retention times of (a) 9.65, (b) 10.89, (c) 9.47, and (d) 9.27 shown in (A). CTZ coelenterazine, CTMD coelenteramide, CTM coelenteramine. These chemical structures are shown in Supplementary Fig. S5.Full size image More

1.Hansen, L. J. & Hoffman, J. R. Climate Savvy: Adapting Conservation and Resource Management to a Changing World, Austral Ecology (Island Press, 2011).Book 

Google Scholar 
2.Huang, Q., Fleming, C., Robb, B., Lothspeich, A. & Songer, M. How different are species distribution model predictions? Application of a new measure of dissimilarity and level of significance to giant panda Ailuropoda melanoleuca. Ecol. Inform. 46, 114–124 (2018).Article 

Google Scholar 
3.Giorgi, F. & Mearns, L. O. Probability of regional climate change based on the Reliability Ensemble Averaging (REA) method. Geophys. Res. Lett. 30, 1–4 (2003).Article 

Google Scholar 
4.Palmer, T. N., Doblas-Reyes, F. J., Hagedorn, R. & Weisheimer, A. Probabilistic prediction of climate using multi-model ensembles: From basics to applications. Philos. Trans. R. Soc. London B Biol. Sci. 360, 1991–1998 (2005).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
5.Dormann, C. F. et al. Prediction uncertainty of environmental change effects on temperate European biodiversity. Ecol. Lett. 11, 235–244 (2008).PubMed 
Article 

Google Scholar 
6.Spence, M. A. et al. A general framework for combining ecosystem models. Fish Fish. 19, 1031–1042 (2018).Article 

Google Scholar 
7.Janssen, P. H. M., Petersen, A. C., van der Sluijs, J. P., Risbey, J. S. & Ravetz, J. R. A guidance for assessing and communicating uncertainties. Water Sci. Technol. 52, 125–131 (2005).CAS 
PubMed 
Article 

Google Scholar 
8.Frewer, L. The public and effective risk communication. Toxicol. Lett. 149, 391–397 (2004).CAS 
PubMed 
Article 

Google Scholar 
9.Leiserowitz, A. A., Malbach, E. W., Roser-Renough, C., Smith, N. & Dawson, E. Climategate, public opinion, and the loss of trust. Am. Behav. Sci. 57, 818–837 (2012).Article 

Google Scholar 
10.Hyder, K. et al. Making modelling count: Increasing the contribution of shelf-seas community and ecosystem models to policy development and management. Mar. Policy 61, 291–302 (2015).Article 

Google Scholar 
11.Cartwright, S. J. et al. Communicating complex ecological models to non-scientist end users. Ecol. Model. 338, 51–59 (2016).Article 

Google Scholar 
12.Doyle, E. E. H., Johnston, D. M., Smith, R. & Paton, D. Communicating model uncertainty for natural hazards: A qualitative systematic thematic review. Int. J. Disaster Risk Reduct. 33, 449–476 (2019).Article 

Google Scholar 
13.Kamal, A. et al. Recent advances and challenges in uncertainty visualization: A survey. J. Vis. https://doi.org/10.1007/s12650-021-00755-1 (2021).Article 

Google Scholar 
14.Spiegelhalter, D. J. & Riesch, H. Don’t know, can’t know: Embracing deeper uncertainties when analysing risks. Philos. Trans. R. Soc. A 369, 4730–4750 (2011).ADS 
MathSciNet 
Article 
MATH 

Google Scholar 
15.Brodlie, K., Osorio, R. A. & Lopes, A. A review of uncertainty in data visualization. In Expanding the Frontiers of Visual Analytics and Visualization (eds Dill, J. et al.) 81–109 (Springer, 2012).Chapter 

Google Scholar 
16.MacEachren, A. M. et al. Visualizing geospatial information uncertainty: what we know and what we need to know. Cartogr. Geogr. Inf. Sci. 32, 139–160 (2005).Article 

Google Scholar 
17.Ibrekk, H. & Morgan, M. G. Graphical communication of uncertain quantities to nontechnical people. Risk Anal. 7, 519–529 (1987).Article 

Google Scholar 
18.Hawkins, E. The cascade of uncertainty in climate projections. https://www.climate-lab-book.ac.uk/2014/cascade-of-uncertainty/. 6 Feb 2014.19.Wilby, R. L. & Dessai, S. Robust adaptation to climate change. Weather 65, 180–185 (2010).ADS 
Article 

Google Scholar 
20.Daron, J., Lorenz, S., Taylor, A. & Dessai, S. Communicating future climate projections of precipitation change. Clim. Change. 166, 23 (2021).ADS 
Article 

Google Scholar 
21.Kinkeldey, C., MacEachren, A. M. & Schiewe, J. How to assess visual communication of uncertainty? A systematic review of geospatial uncertainty visualisation user studies. Cartogr. J. 51, 372–386 (2014).Article 

Google Scholar 
22.van Vuuren, D. P. et al. The representative concentration pathways: An overview. Clim. Change 109, 5–31 (2011).ADS 
Article 

Google Scholar 
23.Breslow, N. E. & Clayton, D. G. Approximate inference in generalized linear mixed models. J. Am. Stat. Assoc. 88, 9–25 (1993).MATH 

Google Scholar 
24.Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
25.Christensen, R.H.B. ordinal—Regression models for ordinal data. R package version 2018.4-19. https://CRAN.R-project.org/package=ordinal Deposited 19 Apr 2018.26.Peterson, B. & Harrell, F. E. Partial proportional odds models for ordinal response variables. Appl. Stat. 39, 205–217 (1990).Article 
MATH 

Google Scholar 
27.Momeni, A., Pincus, M. & Libien, J. Introduction to Statistical Methods in Pathology (Springer, 2018).Book 

Google Scholar 
28.Turner, H. & Firth, D. Bradley-Terry models in R: The BradleyTerry2 package. J. Stat. Softw. 48, 1–21 (2012).Article 

Google Scholar 
29.Akaike, H. Information theory and an extension of maximum likelihood principle. In Proceedings of the 2nd International Symposium on Information Theory. Tsahkadsor, Armenia, USSR, 2–8 September (1973).30.Rubinstein, R. The cross-entropy method for combinatorial and continuous optimization. Methodol. Comput. Appl. Probab. 1, 127–190 (1999).MathSciNet 
Article 

Google Scholar 
31.Kendall, M. G. A new measure of rank correlation. Biometrika 30, 81–93 (1938).Article 
MATH 

Google Scholar 
32.Venables, W. N. & Ripley, B. D. Modern Applied Statistics with S (Springer, 2002).Book 
MATH 

Google Scholar 
33.Firth, D., qvcalc: Quasi variances for factor effects in statistical statistical models. R package version 0.9-1. https://CRAN.R-project.org/package=qvcalc. Deposited 19 Sept 2017.34.Pihur, V., Datta, S. & Datta, S. RankAggreg: Weighted rank aggregation (2018). R package version 0.6-4. https://CRAN.R-project.org/package=RankAggreg. Deposited 18 Mar 2018.35.Belia, S., Fidler, F., Williams, J. & Cumming, G. Researchers misunderstand confidence intervals and standard error bars. Psychol. Methods 10, 389 (2005).PubMed 
Article 

Google Scholar 
36.Cumming, G., Fidler, F. & Vaux, D. L. Error bars in experimental biology. J. Cell Biol. 177, 7–11 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Hullman, J., Resnick, P. & Adar, E. Hypothetical outcome plots outperform error bars and violin plots for inferences about reliability of variable ordering. PLoS ONE 10, e0142444 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
38.Correll, M. & Gleicher, M. Error bars considered harmful: Exploring alternate encodings for mean and error. IEEE Trans. Vis. Comput. Graph. 20, 2142–2151 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
39.Lipkus, I. M. & Hollands, J. G. The visual communication of risk. JNCI Monogr. 25, 149–163 (1999).Article 

Google Scholar 
40.Spiegelhalter, D., Pearson, M. & Short, I. Visualizing uncertainty about the future. Science 333, 1393–1400 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
41.Few, S. Solutions to the problem of over-plotting in graphs. Visual Business Intelligence Newsletter (2008).42.Daron, J. D., Lorenz, S., Wolski, P., Blamey, R. C. & Jack, C. Interpreting climate data visualisations to inform adaptation decisions. Clim. Risk Manag. 10, 17–26 (2015).Article 

Google Scholar 
43.Heer, J. & Agrawala, M. Multi-scale banking to 450. IEEE Trans. Vis. Comput. Graph. 12, 701–708 (2006).PubMed 
Article 

Google Scholar 
44.Lorenz, S., Dessai, S., Paavola, J. & Forster, P. M. The communication of physical science uncertainty in European National Adaptation Strategies. Clim. Change 132, 143–155 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
45.Quispel, A., Maes, A. & Schilperoord, J. Graph and chart aesthetics for experts and laymen in design: The role of familiarity and perceived ease of use. Inf. Vis. 15, 238–252 (2016).Article 

Google Scholar 
46.Saary, M. J. Radar plots: A useful way for presenting multivariate health care data. J. Clin. Epidemiol. 61, 311–317 (2008).PubMed 
Article 

Google Scholar 
47.Vaughan, N. E. & Gough, C. Expert assessment concludes negative emissions scenarios may not deliver. Environ. Res. Lett. 11, 95003 (2016).Article 
CAS 

Google Scholar 
48.Peltier, J. Excel Charting Dos and Don’ts (Peltier Technical Services, 2013).
Google Scholar 
49.Lohse, G. L. The role of working memory on graphical information processing. Behav. Inf. Technol. 16, 297–308 (1997).Article 

Google Scholar 
50.Grainger, S., Mao, F. & Buytaert, W. Environmental data visualisation for non-scientific contexts: Literature review and design framework. Environ. Model. Softw. 85, 299–318 (2016).Article 

Google Scholar 
51.Few, S. Heatmaps: to bin or not to bin? Visual Business Intelligence Newsletter (2017).52.Moreland, K. Diverging color maps for scientific visualization. In Advances in Visual Computing (eds Bebis, G. et al.) 92–103 (Springer, 2009).Chapter 

Google Scholar 
53.Alrehiely, M., Eslambolchilar, P. & Borgo, R. Evaluating different visualization designs for personal health data. In Proceedings of the 32nd International BCS Human Computer Interaction Conference, Belfast, UK, 4–6 July (2018).54.Saket, B., Endert, A. & Demiralp, C. Task-based effectiveness of basic visualizations. IEEE Trans. Visual Comput. Graph. 25, 2505–2512 (2019).Article 

Google Scholar 
55.Dahshan, M., Polys, N. F., Jayne, R. S. & Pollyea, R. M. Making sense of scientific simulation ensembles with semantic interaction. Comput. Graph. Forum. 39, 325–343 (2020).Article 

Google Scholar  More

Water useMany potato physiological features (photosynthesis intensity, leaf water potential) morphological and agronomic features as the Soil Plant Analysis Development (SPAD) and dry matter content can be used as indicators of potato water stress. In this result water consumption and the average daily amount of water used for irrigation differed over the growing season, but differences also occurred between varieties and the humidity level (Table 1). When irrigating the Julinka variety at all stages of the growing season, regardless of the established pF values, water consumption per pot was higher. The average dose of water supplied per pot was 9.7%, 30.7% and 26.6% greater than for the Denar variety, at humidity levels 1, 2 and 3, respectively. The highest water consumption was observed during the potato growth period from BBCH 40/400 to 69/609 and ranged from 0.39 l/pot /day (level 1) to 0.99 l/ pot/day (level 3).Table 1 Water consumption per pot within potato growing stages (in liters) and average consumption of water per pot (in brackets).Full size tableThe highest water consumption in both potato varieties occurred in July (11–18 July). Analyzing the remaining two months of the irrigation period, it can be seen that in June the plants used less water than in July. Seasonal irrigation doses in mid-early potato of studies of Rolbiecki et al. (2015)9 ranged from 40 to 170 mm, and the highest daily values of field water consumption (over 3 mm) occurred in July, similar to the results in this research.Depending on the irrigation system, water consumption efficiency in potato varies from 5.4 to 12 kg m−316,24. Drip irrigation is one of the most effective methods and ranged from 6.3 to 8.6 kg m−3 (Sharma 2007)25. Different values for average WUE index’ in potato cultivation were obtained by Ati et al. (2012)26, and indicated value ranged from 5.9 to 12.2 kg m−3. In present research, average WUE index’ for the Denar variety was from 0.00 l day−1 in the 1st period to 0.79 l day−1 in the 5th harvest period, while for the Julinka it was from 0.49 to 0.92 l day−1, respectively.In the research by Zin El Abedin et al. (2019)27 the amount of water used for irrigating potato amounted to 1505 mm and 1062 mm for FI (full irrigation) and PRD (partial root zone drying) variants, respectively. The use of 50% of water consumption in the PRD reduced water productivity (WP), as compared to water stress in the form of excess FI and deficit irrigation (DI). A large amount of water in conditions of water deficit causes losses due to evaporation and leads to degradation of the soil environment. In turn, in this research the highest water consumption in both varieties was found at level 3, 39.60 l for the Denar variety and 50.15 l for the Julinka variety.Pszczółkowski et al. (2009)28 showed that early potato varieties water requirements in the period from May 1 to August 31 amounted to 336.4 mm, with greatest requirements in July (108—119.6 mm). In our research, the amount of water used depended on the assumed humidity level and amounted from 19.60 × 103 to 39.60 × 103 cm3 for the Denar variety and between 21.50×103 to 50.15 × 103 cm3 for the Julinka (Table 1).Total potato and tuber massThe total weight of plants aboveground—(stems with leaves) and underground (tubers, stolons and roots) was greater in water humidity level 1 than in humidity levels 2 and 3 (Table 2). Administration of increased amounts of water in the later stages of potato growth resulted in inhibition of biomass growth, mainly for the Julinka variety. At the 5th harvest time, at humidity level 3, the total weight and the weight of tubers were 59.2% and 54.7% lower than those obtained at level 1, respectively. At the same time, the difference for Denar was 11.9% and 18.8%, respectively. Begum et al. (2015, 2018)16,22 and Reyes-Cabrera et al. (2016)5 showed that the production of total and commercial tuber yield was strongly dependent on the total biomass production and its structure.Table 2 Potato total biomass and tuber increase depending on water humidity level (g per plant).Full size tableA three-factor analysis of variance showed that the total weight as well as the weight of potato tubers differed significantly by the humidity level and the variety. A significant effect was found for humidity level on the total weight and tuber weight for the Denar variety and tuber weight for the Julinka variety (Table 3).Table 3 Variance analysis for total biomass and tuber of potato depending on factors (significance verified by the Fisher test).Full size tableAnalysis of variance showed a significant impact of the variety on potato plant weight, while it did not show significant interaction of weight and weight of tubers between measurement dates. No significant effect was obtained for interaction between the factors studied (Table 3).Wang et al. (2009)29, concluded that the use of irrigation significantly contributed to an increase total and commercial tubers of medium-early Folva variety yield and its quality. Ossowski et al. (2013)30, shown that irrigation had a significant effect on medium-early potato varieties: Barycz, Mors, Triada tuber yield. When using drip irrigation, yield increased by 26%. In turn, Mazurczyk et al. (2007)31 showed that drip irrigation increased the tuber yield from 29.4–37.5 to 45.1–54.4 t·ha−1.Over the period from the 1st to the 5th harvest date, the total plant biomass increased from 3.5-fold (Julinka—level 3) to 7.2-fold (Julinka—level 1). On the first harvest, Denar did not produce tubers at levels 1 and 3, and for level 2 its weight was the lowest (6 g from a pot). The increase in tuber weight to the last harvest date was the highest for level 2: 23.9- and 22.9-fold, in Denar and Julinka varieties, respectively. At level 3, the growth dynamics of tubers was the lowest: 11.7 times for the Julinka and 9.1 times for the Denar variety (measured from the second harvest date). The highest total biomass increases and tuber weight was found between the 3rd, 4th and 5th dates when humidity was at levels 1 and 2, and between the 3rd and 4th dates at level 3.Kumari et al. (2011, 2018)1,2 concluded that drip irrigation significantly contributed to an increase in potato tuber yield 18% greater than with other irrigation methods. Xu et al. (2010)32 achieved higher yields using the same irrigation system (40–48 t ha−1), and potato tuber weight was reduced under the slight water stress. Potato reacts to stress when soil water tension exceeds 20 kPa24. In a study by Amer et al. (2016)33 potato tuber yield also decreased with the application of excessive irrigation, resulting in greater stress, increased vegetative growth and potential leaching of nutrients from the root zone.Changes also occur in the quality of potato tubers, such as the shape, skin smoothness and chemical composition34.In the research carried out by Zin El-Abedin et al. (2019)27 differences were found in potato tuber yield depending on the irrigation variant. At FI, the highest tuber yields of 31.77–35.91 Mg ha−1 were obtained. Water deficiency reduced tuber yield, in DI variants, by 53.24–65.15% as compared to the FI. Similar results were obtained by Kumari et al. (2011)1. In the present research, the tuber weight of the Denar variety in the fifth term in level 1, increased by 26% compared to the irrigation at level 2 and was a 24% increase for the Julinka variety under similar conditions. At humidity level 3 there was a decrease in total biomass by 12% and 59% (for Denar and Julinka, respectively) in comparison obtained at level 1. In the research Liu et al. (2006)35 the aboveground biomass reached the highest values in excess water conditions.Potato varieties react differently to the humidity of the soil. Mahmood et al. (2016)36 response of potato varieties diversity to soil water deficit, also Hassanapanah (2010)17 showed the reaction of potato varieties to stress conditions. In our study, a higher total and tubers weight was found for the Julinka variety than for the Denar variety.Regardless of the humidity level and variety, the trends in the biomass yield structure were similar (Fig. 3). A downward trend from the 1st to 5th harvest period was shown for roots and stolons. This varied from 5 to 18% at the beginning of the study to 2–5% by the 5th period. It should be noted that under level 3, especially for the Denar variety, the percentage of roots and stolons was at a constant, low level. The percentage of stems with leaves decreased from 68–90% at the first harvest time to 40–55% at 5th. The dynamics of the decline in the share of stems and leaves was highest at humidity level 3. The tuber percentage was from 0 to 20% for the 1st period to 40–60% for the 5th period.Figure 3Potato biomass structure changes depending on humidity level and tuber harvest term (percentage).Full size imageThe Denar variety, regardless of the humidity level, was characterized by a greater share of stems and leaves. For the Julinka, the tuber percentage at the last harvest was at the same or higher than in the case of stems and leaves. At humidity levels 2 and 3, tubers accounted up to 60% of the harvested biomass.The growth of stem and stolon biomass was noticeable at all stages of potato development (Table 2); greater dynamics were found in the growth of tuber mass (Fig. 3). Under level 3, the growth of the biomass of stems with leaves and stolons was slower than in level 2 of water was used.Water use efficiencyAverage daily doses of water used for the Denar and Julinka varieties in potato harvesting periods are shown in Fig. 4. The volume of water was determined each time for the corresponding level of humidity (1, 2 and 3). Based on the data obtained, a proportional increase in water consumption was found for both potato varieties. The most intensive increase in water consumption was noted at humidity level 3. The W index corresponding to the average daily dose of water calculated for the Denar variety varied from 0.40 l day−1 in the 1st period (O1) to 0.79 l day−1 in the 5th harvest period (O5), whereas for the Julinka it was from 0.49 l day−1 (O1) to 0.92 l day−1 (O5). The W values for the level 3 changed for the Denar variety from 0.23 l day−1 in (O1) to 0.38 l day−1 (O5), while for the Julinka from 0.28 l day-1 (O1) to 0.28 l day-1, respectively (O5). The difference in the intensity of water consumption increase for humidity levels was expressed by varying the values of simple directional coefficients approximating empirical data. The highest values of these coefficients were obtained for the humidity level 1. The directional coefficient for the Denar was 0.0077 day−1, and for the Julinka variety 0.009 day−1. For humidity level 3, these values are 4 and 6 times lower: 0.002 day−1 (Denar) and 0.0014 day−1 (Julinka), respectively.Figure 4Average daily water consumption for potato varieties, at three soil humidity levels (1, 2, 3) and in each of five growing stages (O1), (O2), (O3), (O4), (O5).Full size imageThe average daily water consumption throughout the growing season calculated from potato planting is shown in Fig. 5. The average daily water use was the highest for both varieties at humidity level 3. Index W1 for the Denar was 0.53 l day−1, while for Julinka was higher—0.70 l day−1. The water consumption for the humidity level 1 was about 2 times lower: for the Denar—0.27 l day−1 and for Julinka—0.29 l day−1.Figure 5Average daily water consumption for potato varieties, at three soil humidity levels (1, 2, 3), cumulative calculation from potato planting.Full size imageAhmadi et al. (2017)37 used various irrigation schedule strategies for water demand measurements at evapotranspiration. Water demand has been fully or partially satisfied in static and dynamic modes. The research presents dynamics of vapor pressure deficit (VPD) throughout the growing season. The value of VPD in the first days after planting the potato was about 0.5 kPa while in 70 days maximum value was noted (2.5 kPa), and at the end of the growing season (after 150 days) about 1.5 kPa. Due to the shorter potato growing season in present research, no decrease in water demand was noticed up to about 70 days and, as in the results of the research presented by Ahmadi et al. (2017)37, a steady increase in water demand was noted. Similar results were obtained by King et al. (2020)38 and the largest water deficit was found in the middle of vegetation, after 70–80 days after planting35,39.Values for average daily increase in potato tuber weight (index W2) in individual vegetation periods are presented below (Fig. 6). No approximation of functional models to empirical data is possible; hence, the conclusions are based on a description. In the 1st period, i.e. until day 24 (O1), tuber weight gains were smaller than in the other periods. Depending on the humidity level, these amounted to 2.0 to 3.5 g day−1 for the Denar variety, and 2.7 to 3.9 g day−1 for the Julinka. The differences for Denar were 1.5 g day−1 and for Julinka 1.2 g day−1. In the 2nd irrigation period (O2), average daily increase in potato tuber weight was the highest, from 5.9 g day−1 for level 2 to 7.9 g day−1 for level 3. Average daily tuber weight gain was 13% higher for level 1 than for level 2.Figure 6Average daily potato varieties tuber increase, at three soil humidity levels (1, 2, 3), in each of five potato growing stages (O1), (O2), (O3), (O4), (O5).Full size imageThe average daily weight gain of tubers of potato varieties (W3), calculated incrementally from the beginning of the experiment (Fig. 7). For the entire growing season, this indicator for the Denar variety was the highest for the humidity level 1st (5.7 g day−1), at the level 3rd (5.1 g day−1) and the lowest at the level 2nd (4.3 g day−1). The average daily weight gain of potato tubers of the Julinka was definitely highest for the first humidity level (8.1 g day−1).Figure 7Average daily potato varieties tuber increase, at three soil humidity levels (1, 2, 3), cumulative calculation from potato planting.Full size imageThe ratio of the average daily water consumption to the average weight gain of potato tuber (W4) for individual periods is given in Fig. 8. For humidity level 1 for Denar and Julinka varieties, the values decreased with the growing period of vegetation. In the period (O1), 0.079 l of water was used for the Denar variety and 0.075 l for the Julinka for an increase in potato tuber weight of 1 g. In the next stages of the growing season, this index ranged from 0.35 to 0.45 l g−1 for the Denar variety, for the Julinka it was definitely smaller and range from 0.25 to 0.34 l g−1. At humidity level 1, Julinka used less water than Denar to produce the same weight of tubers. At humidity level 2, the volume of water used at the beginning of growth was also the largest for the Denar variety (0.159 l g−1). This amount was two times higher than the volume at level 1. In subsequent periods, the indicator changed and ranged from 0.059 to 0.105 l g−1. For the Julinka variety, water consumption varied in individual periods from 0.085 to 0.113 l g−1 and showed no trend. At humidity level 3, Denar used the greatest amount of water, as compared to levels 1 and 2, and showing no trend. The Julinka variety used even more water at the same humidity level. This amount ranged from 0.164 to 0.298 l g−1 and, unlike in previous cases, it showed an upward trend with plant development.Figure 8Ratio of average daily water consumption to average daily tuber mass increase dependent on three soil humidity levels (1), (2), (3), in each of five potato growing stages (O1), (O2), (O3), (O4), (O5).Full size imageJovanovic et al. (2010)40 divided the potato growing season into five stages related to growth phases. There were no increases in the weight of leaves and stems, while the tuber weight, regardless of the irrigation method (PRD and FI), increased steadily. The weight of tubers in the last harvest, as compared to the first, increased five-fold. A similar relationship was obtained in the work of Shahnazari et al. (2007)41. This research also took account of different levels of humidity using the strategies of PRD and FI, also considering soil retention characteristics (pF curve). The research showed a clear steady increase in potato tuber weight in each harvest.The ratio of the average daily water consumption to the average weight gain of potato tuber varieties calculated cumulatively from the planting (Fig. 9). The W5 value (0.114 l g−1) for the Denar variety at the end of the growing season was the highest for the 3rd humidity level and was about two times higher than at level 1. Water consumption efficiency for the Denar variety was the highest at humidity level 1. The sequence of W5 values is similar for the Julinka, with the difference that for the 3rd level it was 0.205 l g−1; i.e. six times higher than the indicator for level 1. Water consumption efficiency for the Julinka variety was definitely highest at humidity level 1.Figure 9Ratio of average daily water consumption to average daily tuber mass increase dependent on three soil humidity levels (1), (2), (3), cumulative calculation from potato planting.Full size imageBadr et al. (2010)42 analyzed the tuber yield, using two irrigation systems: surface and subsurface drip line. The total volume of water applied during the growing season was the differentiating factor. Results showed that as the volume of water applied during the growing season increased, the yield increased. When the subsurface line was used, applying 75 mm of water during the growing season, the total yield was approx. 27.5 t ha−1, and 32.5 t ha−1 for 325 mm. The effect of water amount on increase in yield was greater for the surface drip line. After applying 75 mm, the yield was 17.5 t ha−1, and 40 t ha−1 (for 325 mm). Similar results were obtained in the work of Linker et al. (2016)43. Regardless of the frequency, amount and total size of irrigation treatments, a proportional increase in the size of crops was observed with increasing doses of water.Shahnazari et al. (2007)41 planned several harvest dates (H0–H4) throughout the entire growing season, analyzing the irrigation efficiency indicator (average WUE index’). Regardless of the irrigation technique, and taking into account, above all, the amount of water administered, the value of the average WUE index’ indicator was the highest in the period H2–H3, similar results were found in our own research. More

1.Parmesan, C. & Yohe, G. A globally coherent fingerprint of climate change impacts across natural systems. Nature 421, 37–42 (2003).CAS 
Article 

Google Scholar 
2.Root, T. L. et al. Fingerprints of global warming on wild animals and plants. Nature 421, 57–60 (2003).CAS 

Google Scholar 
3.Scheffers, B. R. et al. The broad footprint of climate change from genes to biomes to people. Science 354, aaf7671 (2016).
Google Scholar 
4.Cohen, J. M., Lajeunesse, M. J. & Rohr, J. R. A global synthesis of animal phenological responses to climate change. Nat. Clim. Change 8, 224–228 (2018).
Google Scholar 
5.Piao, S. et al. Plant phenology and global climate change: current progresses and challenges. Glob. Change Biol. 25, 1922–1940 (2019).
Google Scholar 
6.Merilä, J. & Hendry, A. P. Climate change, adaptation, and phenotypic plasticity: the problem and the evidence. Evol. Appl. 7, 1–14 (2014).
Google Scholar 
7.Chevin, L. M. & Hoffmann, A. A. Evolution of phenotypic plasticity in extreme environments. Philos. Trans. R. Soc. Lond. B Biol. Sci. 372, 20160138 (2017).
Google Scholar 
8.Fox, R. J., Donelson, J. M., Schunter, C., Ravasi, T. & Gaitán-Espitia, J. D. Beyond buying time: the role of plasticity in phenotypic adaptation to rapid environmental change. Philos. Trans. R. Soc. Lond. B Biol. Sci. 374, 20180174 (2019).
Google Scholar 
9.Thackeray, S. J. et al. Trophic level asynchrony in rates of phenological change for marine, freshwater and terrestrial environments. Glob. Change Biol. 16, 3304–3313 (2010).
Google Scholar 
10.Kharouba, H. M. et al. Global shifts in the phenological synchrony of species interactions over recent decades. Proc. Natl Acad. Sci. USA 115, 5211–5216 (2018).CAS 

Google Scholar 
11.Radchuk, V. et al. Adaptive responses of animals to climate change are most likely insufficient. Nat. Commun. 10, 3109 (2019).
Google Scholar 
12.Visser, M. E. & Gienapp, P. Evolutionary and demographic consequences of phenological mismatches. Nat. Ecol. Evol. 3, 879–885 (2019).
Google Scholar 
13.Kharouba, H. M. & Wolkovich, E. M. Disconnects between ecological theory and data in phenological mismatch research. Nat. Clim. Change 10, 406–415 (2020).
Google Scholar 
14.Samplonius, J. M. et al. Strengthening the evidence base for temperature-mediated phenological asynchrony and its impacts. Nat. Ecol. Evol. 5, 155–164 (2021).
Google Scholar 
15.Charmantier, A. et al. Adaptive phenotypic plasticity in response to climate change in a wild bird population. Science 320, 800–803 (2008).CAS 

Google Scholar 
16.Tomotani, B. M. et al. Climate change leads to differential shifts in the timing of annual cycle stages in a migratory bird. Glob. Change Biol. 24, 823–835 (2018).
Google Scholar 
17.Moyes, K. et al. Advancing breeding phenology in response to environmental change in a wild red deer population. Glob. Change Biol. 17, 2455–2469 (2011).
Google Scholar 
18.Lane, J. E., Kruuk, L. E. B., Charmantier, A., Murie, J. O. & Dobson, F. S. Delayed phenology and reduced fitness associated with climate change in a wild hibernator. Nature 489, 554–557 (2012).CAS 

Google Scholar 
19.Todd, B. D., Scott, D. E., Pechmann, J. H. K. & Whitfield Gibbons, J. Climate change correlates with rapid delays and advancements in reproductive timing in an amphibian community. Proc. R. Soc. Lond. B Biol. Sci. 278, 2191–2197 (2011).
Google Scholar 
20.Taylor, S. G. Climate warming causes phenological shift in pink salmon, Oncorhynchus gorbuscha, behavior at Auke Creek, Alaska. Glob. Change Biol. 14, 229–235 (2008).
Google Scholar 
21.Mills, L. S. et al. Camouflage mismatch in seasonal coat color due to decreased snow duration. Proc. Natl Acad. Sci. USA 110, 7360–7365 (2013).CAS 

Google Scholar 
22.Lameris, T. K. et al. Arctic geese tune migration to a warming climate but still suffer from a phenological mismatch. Curr. Biol. 28, 2467–2473.e4 (2018).CAS 

Google Scholar 
23.Singer, M. C. & Parmesan, C. Phenological asynchrony between herbivorous insects and their hosts: signal of climate change or pre-existing adaptive strategy? Philos. Trans. R. Soc. Lond. B Biol. Sci. 365, 3161–3176 (2010).
Google Scholar 
24.Charmantier, A. & Gienapp, P. Climate change and timing of avian breeding and migration: evolutionary versus plastic changes. Evol. Appl. 7, 15–28 (2014).
Google Scholar 
25.Keogan, K. et al. Global phenological insensitivity to shifting ocean temperatures among seabirds. Nat. Clim. Change 8, 313–317 (2018).
Google Scholar 
26.Both, C. & Visser, M. E. Adjustment to climate change is constrained by arrival date in a long-distance migrant bird. Nature 411, 296–298 (2001).CAS 

Google Scholar 
27.Both, C., van Asch, M., Bijlsma, R. G., van den Burg, A. B. & Visser, M. E. Climate change and unequal phenological changes across four trophic levels: constraints or adaptations? J. Anim. Ecol. 78, 73–83 (2009).
Google Scholar 
28.Cresswell, W. & McCleery, R. How great tits maintain synchronization of their hatch date with food supply in response to long-term variability in temperature. J. Anim. Ecol. 72, 356–366 (2003).
Google Scholar 
29.Visser, M. E., Van Noordwijk, A. J., Tinbergen, J. M. & Lessells, C. M. Warmer springs lead to mistimed reproduction in great tits (Parus major). Proc. R. Soc. Lond. B Biol. Sci. 265, 1867–1870 (1998).
Google Scholar 
30.Sanz, J. J., Potti, J., Moreno, J., Merino, S. & Frías, O. Climate change and fitness components of a migratory bird breeding in the Mediterranean region. Glob. Change Biol. 9, 461–472 (2003).
Google Scholar 
31.Marrot, P., Charmantier, A., Blondel, J. & Garant, D. Current spring warming as a driver of selection on reproductive timing in a wild passerine. J. Anim. Ecol. 87, 754–764 (2018).
Google Scholar 
32.Burgess, M. D. et al. Tritrophic phenological match–mismatch in space and time. Nat. Ecol. Evol. 2, 970–975 (2018).
Google Scholar 
33.Visser, M. E., Holleman, L. J. M. & Gienapp, P. Shifts in caterpillar biomass phenology due to climate change and its impact on the breeding biology of an insectivorous bird. Oecologia 147, 164–172 (2006).
Google Scholar 
34.Pearce-Higgins, J. W., Yalden, D. W. & Whittingham, M. J. Warmer springs advance the breeding phenology of golden plovers Pluvialis apricaria and their prey (Tipulidae). Oecologia 143, 470–476 (2005).CAS 

Google Scholar 
35.Nussey, D. H., Clutton-Brock, T. H., Elston, D. A., Albon, S. D. & Kruuk, L. E. B. Phenotypic plasticity in a maternal trait in red deer. J. Anim. Ecol. 74, 387–396 (2005).
Google Scholar 
36.Husby, A. et al. Contrasting patterns of phenotypic plasticity in reproductive traits in two great tit (Parus major) populations. Evolution 64, 2221–2237 (2010).
Google Scholar 
37.Matthysen, E., Adriaensen, F. & Dhondt, A. A. Multiple responses to increasing spring temperatures in the breeding cycle of blue and great tits (Cyanistes caeruleus, Parus major). Glob. Change Biol. 17, 1–16 (2011).
Google Scholar 
38.Fisher, J. I., Mustard, J. F. & Vadeboncoeur, M. A. Green leaf phenology at Landsat resolution: scaling from the field to the satellite. Remote Sens. Environ. 100, 265–279 (2006).
Google Scholar 
39.Duparc, A. et al. Co-variation between plant above-ground biomass and phenology in sub-alpine grasslands. Appl. Veg. Sci. 16, 305–316 (2013).
Google Scholar 
40.Hinks, A. E. et al. Scale-dependent phenological synchrony between songbirds and their caterpillar food source. Am. Nat. 186, 84–97 (2015).
Google Scholar 
41.Lambrechts, M. M., Blondel, J., Maistre, M. & Perret, P. A single response mechanism is responsible for evolutionary adaptive variation in a bird’s laying date. Proc. Natl Acad. Sci. USA 94, 5153–5155 (1997).CAS 

Google Scholar 
42.Dawson, A. Control of the annual cycle in birds: endocrine constraints and plasticity in response to ecological variability. Philos. Trans. R. Soc. Lond. B Biol. Sci. 363, 1621–1633 (2008).
Google Scholar 
43.Visser, M. E. et al. Phenology, seasonal timing and circannual rhythms: towards a unified framework. Philos. Trans. R. Soc. Lond. B Biol. Sci. 365, 3113–3127 (2010).CAS 

Google Scholar 
44.Caro, S. P., Schaper, S. V., Hut, R. A., Ball, G. F. & Visser, M. E. The case of the missing mechanism: how does temperature influence seasonal timing in endotherms? PLoS Biol. 11, e1001517 (2013).CAS 

Google Scholar 
45.Bourgault, P., Thomas, D., Perret, P. & Blondel, J. Spring vegetation phenology is a robust predictor of breeding date across broad landscapes: a multi-site approach using the Corsican blue tit (Cyanistes caeruleus). Oecologia 162, 885–892 (2010).
Google Scholar 
46.Bison, M. et al. Best environmental predictors of breeding phenology differ with elevation in a common woodland bird species. Ecol. Evol. 10, 10219–10229 (2020).
Google Scholar 
47.Bernhardt, J. R., O’Connor, M. I., Sunday, J. M. & Gonzalez, A. Life in fluctuating environments. Philos. Trans. R. Soc. Lond. B Biol. Sci. 375, 20190454 (2020).
Google Scholar 
48.Gienapp, P., Reed, T. E. & Visser, M. E. Why climate change will invariably alter selection pressures on phenology. Proc. R. Soc. Lond. B Biol. Sci. 281, 20141611 (2014).
Google Scholar 
49.Lönnstedt, O. M., McCormick, M. I., Chivers, D. P. & Ferrari, M. C. O. Habitat degradation is threatening reef replenishment by making fish fearless. J. Anim. Ecol. 83, 1178–1185 (2014).
Google Scholar 
50.Pellerin, F., Cote, J., Bestion, E. & Aguilée, R. Matching habitat choice promotes species persistence under climate change. Oikos 128, 221–234 (2019).
Google Scholar 
51.Firth, J. A., Verhelst, B. L., Crates, R. A., Garroway, C. J. & Sheldon, B. C. Spatial, temporal and individual-based differences in nest-site visits and subsequent reproductive success in wild great tits. J. Avian Biol. 49, e01740 (2018).
Google Scholar 
52.Naef-Daenzer, B. & Keller, L. F. The foraging performance of great and blue tits (Parus major and P. caeruleus) in relation to caterpillar development, and its consequences for nestling growth and fledging weight. J. Anim. Ecol. 68, 708–718 (1999).
Google Scholar 
53.Naef-Daenzer, B. Patch time allocation and patch sampling by foraging great and blue tits. Anim. Behav. 59, 989–999 (2000).CAS 

Google Scholar 
54.Bouwhuis, S., Sheldon, B. C., Verhulst, S. & Charmantier, A. Great tits growing old: selective disappearance and the partitioning of senescence to stages within the breeding cycle. Proc. R. Soc. Lond. B Biol. Sci. 276, 2769–2777 (2009).CAS 

Google Scholar 
55.Cole, E. F. & Sheldon, B. C. The shifting phenological landscape: within- and between-species variation in leaf emergence in a mixed-deciduous woodland. Ecol. Evol. 7, 1135–1147 (2017).
Google Scholar 
56.Wint, W. The role of alternative host-plant species in the life of a polyphagous moth, Operophtera brumata (Lepidoptera: Geometridae). J. Anim. Ecol. 52, 439–450 (1983).
Google Scholar 
57.Keller, L. F. & van Noordwijk, A. J. Effects of local environmental conditions on nestling growth in the great tit Parus major L. Ardea 82, 349–362 (1994).
Google Scholar 
58.Wilkin, T. A., Garant, D., Gosler, A. G. & Sheldon, B. C. Density effects on life-history traits in a wild population of the great tit Parus major: analyses of long-term data with GIS techniques. J. Anim. Ecol. 75, 604–615 (2006).
Google Scholar 
59.Wilkin, T. A. & Sheldon, B. C. Sex differences in the persistence of natal environmental effects on life histories. Curr. Biol. 19, 1998–2002 (2009).CAS 

Google Scholar 
60.Gagen, M. et al. The tree ring growth histories of UK native oaks as a tool for investigating chronic oak decline: an example from the Forest of Dean. Dendrochronologia 55, 50–59 (2019).
Google Scholar 
61.Sturrock, R. N. et al. Climate change and forest diseases. Plant Pathol. 60, 133–149 (2011).
Google Scholar 
62.MacColl, A. D. C. The ecological causes of evolution. Trends Ecol. Evol. 26, 514–522 (2011).
Google Scholar 
63.Grant, P. R. & Price, T. D. Population variation in continuously varying traits as an ecological genetics problem. Integr. Comp. Biol. 21, 795–811 (1981).
Google Scholar 
64.Hereford, J. A quantitative survey of local adaptation and fitness trade-offs. Am. Nat. 173, 579–588 (2009).
Google Scholar 
65.Hadfield, J. D. The spatial scale of local adaptation in a stochastic environment. Ecol. Lett. 19, 780–788 (2016).
Google Scholar 
66.Porlier, M. et al. Variation in phenotypic plasticity and selection patterns in blue tit breeding time: between- and within-population comparisons. J. Anim. Ecol. 81, 1041–1051 (2012).
Google Scholar 
67.Hidalgo Aranzamendi, N., Hall, M. L., Kingma, S. A., van de Pol, M. & Peters, A. Rapid plastic breeding response to rain matches peak prey abundance in a tropical savanna bird. J. Anim. Ecol. 88, 1799–1811 (2019).
Google Scholar 
68.Caro, S. P., Lambrechts, M. M., Balthazart, J. & Perret, P. Non-photoperiodic factors and timing of breeding in blue tits: impact of environmental and social influences in semi-natural conditions. Behav. Process. 75, 1–7 (2007).CAS 

Google Scholar 
69.Bourret, A., Bélisle, M., Pelletier, F. & Garant, D. Multidimensional environmental influences on timing of breeding in a tree swallow population facing climate change. Evol. Appl. 8, 933–944 (2015).
Google Scholar 
70.Nussey, D. H., Wilson, A. J. & Brommer, J. E. The evolutionary ecology of individual phenotypic plasticity in wild populations. J. Evol. Biol. 20, 831–844 (2007).CAS 

Google Scholar 
71.Morris, D. W. Toward an ecological synthesis: a case for habitat selection. Oecologia 136, 1–13 (2003).
Google Scholar 
72.Long, R. A. et al. Linking habitat selection to fitness-related traits in herbivores: the role of the energy landscape. Oecologia 181, 709–720 (2016).
Google Scholar 
73.Morris, D. W. Spatial scale and the cost of density-dependent habitat selection. Evol. Ecol. 1, 379–388 (1987).
Google Scholar 
74.Patten, M. A. & Kelly, J. F. Habitat selection and the perceptual trap. Ecol. Appl. 20, 2148–2156 (2010).
Google Scholar 
75.Ponchon, A., Garnier, R., Grémillet, D. & Boulinier, T. Predicting population responses to environmental change: the importance of considering informed dispersal strategies in spatially structured population models. Divers. Distrib. 21, 88–100 (2015).
Google Scholar 
76.Nilsson, A. L. K. et al. Hydrology influences breeding time in the white-throated dipper. BMC Ecol. 20, 70 (2020).
Google Scholar 
77.Nilsson, A. L. K. et al. Location is everything, but climate gets a share: analyzing small-scale environmental influences on breeding success in the white-throated dipper. Front. Ecol. Evol. 8, 542846 (2020).
Google Scholar 
78.Martin, R. O., Cunningham, S. J. & Hockey, P. A. R. Elevated temperatures drive fine-scale patterns of habitat use in a savanna bird community. Ostrich 86, 127–135 (2015).
Google Scholar 
79.Bailey, L. D. et al. Habitat selection can reduce effects of extreme climatic events in a long-lived shorebird. J. Anim. Ecol. 88, 1474–1485 (2019).
Google Scholar 
80.Kirby, K. J. et al. Changes in the tree and shrub layer of Wytham Woods (southern England) 1974–2012: local and national trends compared. Forestry 87, 663–673 (2014).
Google Scholar 
81.Perrins, C. & McCleery, R. Laying dates and clutch size in the great tit. Wilson Bull. 101, 236–253 (1989).
Google Scholar 
82.Wilkin, T. A., Perrins, C. M. & Sheldon, B. C. The use of GIS in estimating spatial variation in habitat quality: a case study of lay-date in the great tit Parus major. Ibis 149, 110–118 (2007).
Google Scholar 
83.Perrins, C. M. Population fluctuations and clutch-size in the great tit, Parus major L. J. Anim. Ecol. 34, 601–647 (1965).
Google Scholar 
84.Wesołowski, T. & Rowiński, P. Timing of bud burst and tree-leaf development in a multispecies temperate forest. For. Ecol. Manage. 237, 387–393 (2006).
Google Scholar 
85.Gibson, C. W. D. in Woodland Conservation and Research in the Clay Vale of Oxfordshire and Buckinghamshire (eds Kirby, K. J. & Write, F. J.) 32–40 (JNCC, 1988).86.Dawkin, H. C. & Field, D. R. B. A Long-Term Surveillance System for British Woodland Vegetation. Commonwealth Forestry Institute, Oxford, Occasional Paper No. 1. (1978).87.Horsfall, A. S. & Kirby, K. J. The Use of Permanent Quadrats to Record Changes in the Structure and Composition of Wytham Woods, Oxfordshire Research and Survey in Nature Conservation No. 1 (JNCC, 1992).88.Wilkin, T. A., King, L. E. & Sheldon, B. C. Habitat quality, nestling diet, and provisioning behaviour in great tits Parus major. J. Avian Biol. 40, 135–145 (2009).
Google Scholar 
89.Van Noordwijk, M. & Purnomosidhi, P. Root architecture in relation to tree–soil–crop interactions and shoot pruning in agroforestry. Agrofor. Syst. 30, 161–173 (1995).
Google Scholar 
90.Bailey, L. D. & van de Pol, M. climwin: an R toolbox for climate window analysis. PLoS ONE 11, e0167980 (2016).
Google Scholar 
91.van de Pol, M. et al. Identifying the best climatic predictors in ecology and evolution. Methods Ecol. Evol. 7, 1246–1257 (2016).
Google Scholar 
92.Simmonds, E. G., Cole, E. F. & Sheldon, B. C. Cue identification in phenology: a case study of the predictive performance of current statistical tools. J. Anim. Ecol. 88, 1428–1440 (2019).
Google Scholar 
93.Oksanen, J. et al. vegan: Community Ecology Package: R Package v.2.5-6 (2019); https://CRAN.R-project.org/package=vegan94.Sturges, H. A. The choice of a class interval. J. Am. Stat. Assoc. 21, 65–66 (1926).
Google Scholar 
95.Hadfield, J. D. MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package. J. Stat. Softw. https://doi.org/10.18637/jss.v033.i02 (2010).96.R Core Team R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2020); http://www.R-project.org/ More

Summary of sequencing dataWe obtained 168 million and 180 million 250 bp reads from Asian Badger and Northern Hog Badger, respectively. After removing transcripts and unigenes below 200 bp, we obtained 335,772 transcripts and 285,159 unigenes belonging to Asian Badger and 413,917 transcripts and 362,075 unigenes belonging to Northern Hog Badger (Table 1). Next, we analysed the length distribution of the unigenes and transcripts in these two species (Fig. 1). Their N50 of transcript length is longer than 1000 bp, and their N50 of unigene length is longer than 600 bp. The average GC content of the transcriptome data of Asian Badge was 52.71%, a value slightly higher than that of the Northern Hog Badger, which was 52.12% (Table 1).Table 1 Summary of the transcriptome of Asian Badgers and Northern Hog Badgers.Full size tableFigure 1Length and quantity distribution of transcripts and unigenes.Full size imageFunctional annotation and classification of the assembled unigenesThe success rate of annotation of these research data in the seven databases is shown in Table 2. In total, 34,150 (ZH) and 31,632 (GH) unigenes had GO terms (Table 2). Among them, there were three GO items related to digestion: positive regulation of the digestive system process (GH and ZH both have one gene), digestive tract development (GH and ZH both have four genes), and digestion (GH has five genes, ZH has three). Next, we compared the GO terms of Asian Badger and Northern Hog Badger transcriptomes and found that the distributions pattern of gene functions from these two species were particularly similar (Fig. 2). This predictable result indicates that there is no bias in the construction of the libraries from the Asian Badger and Northern Hog Badger. For both species, in the three main partitions (cellular component, molecular function, and biological process) of the GO classification, ‘Cellular process’, ‘Binding’ and ‘Metabolic process’, terms were principal individually (Fig. 2). In total, 8915 (ZH) and 10,203 (GH) unigenes had KOG terms (Table 2). In addition, 15,667 (ZH) and 17,823 (GH) were mapped to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Table 2) and grouped into 32 subclasses. Interestingly, the digestive system subcategory contains 695 and 611 unigenes in Asian Badger and Northern Hog Badger, respectively, involving 9 pathways, namely, bile, gastric acid, pancreatic, salivary secretion, carbohydrate, protein, vitamin, fat digestion and absorption, and mineral absorption.Table 2 Gene annotation success rate statistics.Full size tableFigure 2GO term Top20 for GH and ZH.Full size imageAnalysis of orthologous genesThe transcriptome evolution of different species can be understood by comparing transcriptome data. We analysed the possible orthologous genes between the transcriptome of Asian Badger and Northern Hog Badger obtained in this study. We selected a total of 5227 homologous gene pairs from these four species. After 5227 pairs of homologous genes were optimized and screened, 943 orthologous gene pairs were obtained (Supplementary Table S1).To explore whether the genes related to small intestinal digestion in Asian Badger and Northern Hog Badger have undergone adaptive evolution. We can predict the genes that affect the evolution of the two species through selection pressure on orthologous genes12. We selected 473 orthologous gene pairs with Ka/Ks  > 1 called divergent orthologous genes from the Ka/ks analysis results. We obtained 1263 orthologous gene pairs with Ka/Ks  More

1.Eitzel, M. V. et al. Citizen science terminology matters: Exploring key terms. Citiz. Sci. 2(1), 1 (2017).Article 

Google Scholar 
2.European Citizen Science Association (ECSA). ECSA’s Characteristics of Citizen Science. Apr 2020 [cited 3 July 2020]. In ECSA Our Documents [Internet]. https://ecsa.citizen-science.net/wp-content/uploads/2020/05/ecsa_characteristics_of_citizen_science_-_v1_final.pdf.3.Shirk, J. L. et al. Public participation in scientific research: A framework for deliberate design. Ecol. Soc. 17(2), 29 (2012).Article 

Google Scholar 
4.Follett, R. & Strezov, V. An analysis of citizen science based research: Usage and publication patterns. PLoS One 10(11), e0143687 (2015).Article 

Google Scholar 
5.Tauginienė, L. et al. Citizen science in the social sciences and humanities: The power of interdisciplinarity. Palgrave Commun. 6, 89 (2020).Article 

Google Scholar 
6.European Citizen Science Association (ECSA). Ten Principles of Citizen Science. Sept 2015 [cited 1 Dec 2020]. In ECSA Our Documents [Internet]. https://ecsa.citizen-science.net/wp-content/uploads/2020/02/ecsa_ten_principles_of_citizen_science.pdf.7.Vayena, E. & Tasioulas, J. Adapting standards: Ethical oversight of participant-led health research. PLoS Med. 10(3), e1001402 (2013).Article 

Google Scholar 
8.Resnik, D. B., Elliott, K. C. & Miller, A. K. A framework for addressing ethical issues in citizen science. Environ. Sci. Policy 54, 475–481 (2015).Article 

Google Scholar 
9.Riesch, H. & Potter, C. Citizen science as seen by scientists: Methodological, epistemological and ethical dimensions. Public Underst. Sci. 23(1), 107–120 (2014).Article 

Google Scholar 
10.Rothstein, M. A., Wilbanks, J. T. & Brothers, K. B. Citizen science on your smartphone: An ELSI research agenda. J. Law Med. Ethics 43(4), 897–903 (2015).Article 

Google Scholar 
11.Guerrini, C. J., Majumder, M. A., Lewellyn, M. J. & McGuire, A. L. Citizen science, public policy. Science 361(6398), 134–136 (2018).ADS 
CAS 
Article 

Google Scholar 
12.Resnik, D. B. Citizen scientists as human subjects: Ethical issues. Citiz. Sci. 4(1), 11 (2019).
Google Scholar 
13.Woolley, J. P. et al. Citizen science or scientific citizenship? Disentangling the uses of public engagement rhetoric in national research initiatives. BMC Med. Ethics 17, 33 (2016).Article 

Google Scholar 
14.Rasmussen, L. M. Research ethics in citizen science. In The Oxford Handbook of Research Ethics (eds Iltis, A. S. & MacKay, D.) (Oxford University Press, 2021).
Google Scholar 
15.Chesser, S., Poster, M. M. & Tuckett, A. G. Cultivating citizen science for all: Ethical considerations for research projects involving diverse and marginalized populations. Int. J. Soc. Res. Methodol. 23(5), 497–508 (2020).Article 

Google Scholar 
16.Vayena, E. & Tasioulas, J. “We the scientists”: A human right to citizen science. Philos. Technol. 28, 479–485 (2015).Article 

Google Scholar 
17.Rasmussen, L. M. “Filling the ‘Ethics Gap’ in Citizen Science Research”: A Workshop Report. 2017 [cited 15 July 2020]. In NIEHS Partnerships for Environmental Public Health [Internet]. https://www.niehs.nih.gov/research/supported/translational/peph/webinars/ethics/rasmussen_508.pdf.18.Louviere, J. J., Flynn, T. N. & Marley, A. A. J. Best-Worst Scaling: Theory, Methods and Applications (Cambridge University Press, 2015).Book 

Google Scholar 
19.Mühlbacher, A. C., Kaczynski, A., Zweifel, P. & Johnson, F. R. Experimental measurement of preferences in health and healthcare using best-worst scaling: An overview. Health Econ. Rev. 6, 2 (2016).Article 

Google Scholar 
20.Marti, J. A best-worst scaling survey of adolescents’ level of concern for health and non-health consequences of smoking. Soc. Sci. Med. 75, 87–97 (2012).Article 

Google Scholar 
21.Erdem, S. & Rigby, D. Investigating heterogeneity in the characterization of risks using best worst scaling. Risk Anal. 33(9), 1728–1748 (2013).Article 

Google Scholar 
22.Peay, H. L., Hollin, I. L. & Bridges, J. F. P. Prioritizing parental worry associated with Duchenne muscular dystrophy using best-worst scaling. J. Genet. Couns. 25, 305–313 (2016).Article 

Google Scholar 
23.Flynn, T. N., Louviere, J. J., Peters, T. J. & Coast, J. Best-worst scaling: What it can do for health care research and how to do it. J. Health Econ. 26, 171–189 (2007).Article 

Google Scholar 
24.Finn, A. & Louviere, J. J. Determining the appropriate response to evidence of public concern: The case of food safety. J. Public Policy Mark. 11(2), 12–25 (1992).Article 

Google Scholar 
25.Louviere, J. J. & Flynn, T. N. Using best-worst scaling choice experiments to measure public perceptions and preferences for healthcare reform in Australia. Patient 3(4), 275–283 (2010).Article 

Google Scholar 
26.Janssen, E. M., Benz, H. L., Tsai, J.-H. & Bridges, J. F. P. Identifying and prioritizing concerns associated with prosthetic devices for use in a benefit-risk assessment: A mixed methods approach. Expert Rev. Med. Devices 15(5), 385–398 (2018).CAS 
Article 

Google Scholar 
27.Auger, P., Devinney, T. M. & Louvier, J. J. Using best-worst scaling methodology to investigate consumer ethical beliefs across countries. J. Bus. Ethics 70, 299–326 (2007).Article 

Google Scholar 
28.Bridges, J. F. P. et al. Conjoint analysis applications in health—A checklist: A report of the ISPOR Good Research Practices for Conjoint Analysis Task Force. Value Health 14, 403–413 (2011).Article 

Google Scholar 
29.Coast, J. et al. Using qualitative methods for attribute development for discrete choice experiments: Issues and recommendations. Health Econ. 21, 730–741 (2012).Article 

Google Scholar 
30.Johnson, F. R. et al. Constructing experimental designs for discrete-choice experiments: Report of the ISPOR Conjoint Analysis Experimental Design Good Research Practices Task Force. Value Health 16, 3–13 (2013).Article 

Google Scholar 
31.Hauber, A. B. et al. Statistical methods for the analysis of discrete choice experiments: A report of the ISPOR Conjoint Analysis Good Research Practices Task Force. Value Health 19, 300–315 (2016).Article 

Google Scholar 
32.Flynn, T. N. Valuing citizen and patient preferences in health: Recent developments in three types of best-worst scaling. Expert Rev. Pharmacoecon. Outcomes Res. 10(3), 259–267 (2010).Article 

Google Scholar 
33.Louviere, J., Lings, I., Islam, T., Gudergan, S. & Flynn, T. An introduction to the application of (case 1) best-worst scaling in marketing research. Int. J. Res. Mark. 30, 292–303 (2013).Article 

Google Scholar 
34.Jannsen, E. M., Marshall, D. A., Hauber, A. B. & Bridges, J. F. P. Improving the quality of discrete-choice experiments in health: How can we assess validity and reliability?. Expert Rev. Pharmacoecon. Outcomes Res. 17(6), 531–542 (2017).Article 

Google Scholar 
35.Gallego, G., Bridges, J. F. P. & Flynn, T. Using best-worst scaling in horizon scanning for hepatocellular carcinoma technologies. Int. J. Technol. Assess. Health Care 28(3), 339–346 (2012).Article 

Google Scholar 
36.Bowser, A., Shilton, K., Preece, J. & Warrick, E. Accounting for privacy in citizen science: Ethical research in a context of openness. In CSCW ’17: Proceedings of the 2017 ACM Conference on Computer Supported Cooperative Work and Social Computing 2124–2136 (Association for Computing Machinery, 2017).37.Pandya, R. E. A framework for engaging diverse communities in citizen science in the US. Front Ecol. Environ. 10(6), 314–317 (2012).Article 

Google Scholar 
38.Citizen Science Association (CSA). Who We Are; Who We Serve. n.d. [cited 1 Sept 2021]. In CSA About [Internet]. https://www.citizenscience.org/about/.39.Sawtooth Software. The MaxDiff System Technical Paper v. 9. Oct. 2020 [cited 25 Jan 2021] In Technical Papers [Internet]. https://sawtoothsoftware.com/resources/technical-papers/maxdiff-technical-paper.40.English, P. B., Richardson, M. J. & Garzon-Galvis, C. From crowdsourcing to extreme citizen science: Participatory research for environmental health. Annu. Rev. Public Health 39, 335–350 (2018).CAS 
Article 

Google Scholar 
41.Conrad, C. C. & Hilchey, K. G. A review of citizen science and community-based environmental monitoring: Issues and opportunities. Environ. Monit. Assess. 176, 273–291 (2011).Article 

Google Scholar 
42.Cashman, S. B. et al. The power and the promise: Working with communities to analyze data, interpret findings, and get to outcomes. Am. J. Public Health 98(8), 1407–1417 (2008).Article 

Google Scholar 
43.Citizen Science Association (CSA). Trustworthy Data Practices. n.d. [cited 5 Sept 2020]. In CSA Events [Internet]. https://www.citizenscience.org/data-ethics-study/.44.Cheung, K. L. et al. Using best-worst scaling to investigate preferences in health care. Pharmacoecon. 34, 1195–1209 (2016).Article 

Google Scholar 
45.Hastings, J. J. A. When citizens do science. Narrat. Inq. Bioeth. 9(1), 33–34 (2019).Article 

Google Scholar 
46.Guerrini, C. J., Trejo, M., Canfield, I. & McGuire, A. L. Core values of genomic citizen science: Results from a qualitative interview study. BioSocieties https://doi.org/10.1057/s41292-020-00208-2 (2020).Article 

Google Scholar 
47.Sharon, T. Self-tracking for health and the Quantified Self: Re-articulating autonomy, solidarity, and authenticity in the age of personalized medicine. Philos. Technol. 30, 93–121 (2017).Article 

Google Scholar  More

Wet and dry weight of aerial partsDry weight of aerial parts was significantly affected by drought stress and genotype treatments and their interactions (Table 2). With increasing drought stress the amount of dry weight of aerial parts in all genotypes was decreased. Dry functions in I2, I3, and I4 levels in H genotype (Lavandula gngustifolia cv. Hidecot) were 15.68%, 40.35% and 48.15%, respectively. In S genotype (Lavandula stricta) these amounts were 0.78%, 48.58% and 51.72%, respectively; and in M genotype (Lavandula angustifolia cv. Muneasted) they were 22.29%, 49.38% and 52.63%, respectively. Compared to the control group, the most reduction in dry weight of aerial parts was in M genotype. The highest amount of dry weight (11.40 g in plant) was observed in H genotype in drought stress of 90–100% of field capacity. The lowest amount of dry weight of aerial parts (3.07 g) was seen in S genotype in drought stress of 30–40% field capacity (Fig. 2).Table 2 Variance analysis of the effect of drought stress on enzymatic activity of antioxidant enzymes, and quantity of essential oil from different lavender genotypes.Full size tableFigure 2The effect of drought stress on dry weight of aerial parts in different lavender genotypes.Full size imageIn this study drought stress had a negative effect on biomass of lavender plants. This effect can be due to water shortage. Because drought stress cause reduction in swelling, total water potential in cell and withering, it also results in closing stomata, reduction in cell division, and cell enlargement47,48. Reduction in cell division and cell enlargement as a result of drought, reduce the leaf surface, photosynthesis and growth function of the plant. In other words, reduction in photosynthesis products, cause reduction in leaf’s surface; and reduction in transfer of assimilated materials to aerial part, as a result of drought, cause decrease in aerial yield of the plant49. In this regard, Abbaszadeh et al. (2020) reported that due to drought stress of 30% and 60% of field capacity, dry weight of aerial parts in Rosmarinus officinalis L. has decreased. While contrary to our results Rhizopoulou and Diamantoglou (1991) observed that dry weight of leaves from Marjoram plant (Origanum majorana) was increased with increased soil moisture deficiency; which can be due to differences in plant species and ecological conditions50,51.Proline content of leavesThe results of variance analysis showed that drought stress, genotype and their interactions have significantly affected proline content of leaves (Table 2). With increasing drought stress the proline content was increased. The highest amount of proline content (4.96 mg per g) was observed in H genotype in I4 drought level (30–40% of field capacity). While the lowest amount (1.08 mg per g) was observed in S genotype in irrigation of 90–100% of field capacity (Fig. 3). In each genotype separately, in I2 to I3 drought levels the amount of proline was equal, but in H and M genotypes with increasing drought stress, the amount of proline was increased, While in S genotype with increasing water deficit proline did not show a significant increase. This indicates that two genotypes (H and M) have a similar function for using these types of osmolyte to deal with this level of drought. Which this result may be exist another osmolite production as a resistance mechanism in S genotype52.Figure 3The effect of drought stress on proline content in different genotypes of lavender.Full size imageOne change that happens in biological and non-biological stresses is increasing the amount of osmolytes in plant. To prevent negative effects of drought stress, the plant increases the amount of its osmolytes including proline53. Proline is an amino acid which in addition to act as an osmolyte, plays an important role in maintaining and stabilizing membranes by adding membrane phospholipids and changing the hydrated layer around macromolecules. Proline is also recognized as a stabilizer for cellular homeostasis under stressful conditions. This is due to high ability of proline to stabilize sub-cellular structures such as proteins and cell membranes and its ability to eliminate free radicals54. In present study, increasing proline content in different genotypes of lavender as a result of drought, can be for the same reason. It is proved that in some plants, changes in amount of proline is related to their ability to tolerate and adapt with drought stress; so, the proline content can be used as an indicator to select drought-resistant plants. Hosseinpour et al. (2020), reported that in response to drought stress, accumulation of compatible metabolites such as proline can participate in water absorption. In accordance with our results, an increase in proline content in different genotypes of Calendula officinalis plant due to drought stress has been reported as well55,56. However, in some plant species, other osmolites are produced under biological stress, the most important of them is glycine betaine. So that it is probable that the relationship between glycine betaine accumulation and stress tolerance, such as drought stress, is species- or even genotype specific57. As a results, the S genotype likely produced glycine betaine under drought stress, obviously, completed studies are needed to confirm this hypothesis.Relative water content of leavesThe relative water content (RWC) of leaves was significantly affected by drought stress, genotype and their interaction (Table 2). The highest amount of RWC (87.43%) was observed in H genotype in no drought stress condition. The lowest RWC (19.60%) was observed in S genotype in 30–40% of field capacity (Fig. 4). The results of comparing average data showed that in highest level of drought stress RWC in H, S and M genotypes is 57.25%, 65.19% and 58.88%, respectively; which compared to the control group, it is decreased in all genotypes. This suggests higher resistance of H genotype to maintain RWC of leaves (Fig. 4). In all evaluated genotypes, with increasing drought, RWC was decreased.Figure 4The effect of drought stress on RWC of leaves in different lavender genotypes.Full size imageRWC is a suitable indicator for water stress in plants. Drought stress by reducing RWC and total water potential of cell, result in reduction in growth of plants. The osmoregulation mechanisms in drought-resistant plants, maintains high RWC in them. Reduction in RWC of leaves as a result of water deficiency stress, is due to reduction in amount of water in tissue, reduction in amount of water in soil, and the negative soil water potential58. Alinejad et al. (2020), reported that RWC of leaves in Datura stramonium L. plant was decreased due to drought, in a way that the highest amount of RWC (80.22%) was seen in 55% of field capacity, compared to 35% and 15% of field capacity59. Also Mohammadi et al. (2018) suggested that RWC of leaves in Thymu vulgari L. was decreased to 18.41%, after being exposed to drought60.Total phenolic and flavonoids contents in leavesDrought, genotype and their interaction had a significant effect on total phenolic content of leaves (Table 2). The results suggest that in different levels of drought, total phenolic content was different in lavender genotypes. In the highest level of drought, total phenolic content in H, S, and M genotypes was respectively increased 18.64%, 28.57% and 98.07% in comparison with the control group. The highest difference in total phenolic content compared to control group was observed in M genotype (Fig. 5).Figure 5The effect of total phenolic content of leaves in different genotypes of lavender.Full size imageTotal flavonoids content of leaves was significantly (p ≤ 0.01) affected by drought and genotype (Table 2). The results of comparing averages showed that the highest amount of total flavonoids (1.12 mg quercetin per g of fresh weight) was in H genotype, and the lowest amount (0.95 mg quercetin per g of fresh weight) was in M genotype (Table 3). Moreover, our results showed that drought level from I2 to I4 caused an increase of 12.74%, 14.61% and 15.38% in total flavonoid content of leaves, respectively. Which indicates an increase in flavonoid amount with increasing drought level (Table 3). Table 3 Comparing simple effects of genotype and drought stress on traits of lavender plant.Full size tableTotal phenolic content is related to stress-resistance, indirectly by helping cell protection, and directly as an antioxidant61. Phenolic compounds due to their reductive properties, act as a free radical remover62. Our findings are similar to those of a study on growth of Mentha piperita in drought stress54.Total antioxidant activityTotal antioxidant activity was significantly affected by drought stress and genotype (Table 2). With increasing drought, antioxidant activity in H and S genotypes was increased. The results of comparing average data showed that compared to the control group, in drought levels of I2, I3 and I4, antioxidant activity in H genotype was increased by 98.43%, 98.36% and 118.78%, respectively; and in S genotype this amounts were increased by 89.85%, 111.78%, and 131.90% respectively (Table 5). In M genotype the antioxidant activity has reached its highest amount (49.38 mg/g) in I3 level of drought, and then with increasing drought stress the antioxidant activity was decreased, in a way that in highest drought level it had the lowest antioxidant activity (23.18 mg/g). M genotype was used as control (Fig. 6). Our results indicate that in highest drought level, antioxidant activity of S genotype was more than others. Figure 6The effect of drought stress on antioxidant activity in different lavender genotypes.Full size imageAntioxidant enzymesEnzymatic activity of antioxidant enzymes in lavender leaves was significantly affected by genotype and drought stress (Table 2). Our results showed that the highest activity of SOD (304.75 μmol min−1 mg−1 protein) was observed in interaction of H genotype and I4 drought level, and the lowest activity of SOD (144.52 μmol min−1 mg−1 protein) was observed in S genotype with no drought (Fig. 7). Moreover, our observations showed that in I2 and I3 drought levels, the highest amount of SOD enzymatic activity was related to M genotype (Fig. 7). In the highest drought level, enzymatic activity of SOD was increased in H and S genotypes, and it decreased in M genotype.Figure 7The effect of drought stress on enzymatic activity of SOD, POX and CAT in different lavender genotypes.Full size imageEnzymatic activity of peroxidase (POX) enzyme was increased in all three genotypes, with increasing drought. In all drought levels, H genotype had the highest amount of POX activity, compared to other genotypes. There was no significant difference in POX activity in S and H genotypes. The results showed that the highest amount of POX activity (274.48 μmol min−1 mg−1 protein) was observed in interaction of H genotype and 30–40% field capacity, and the lowest amount (117.66 μmol min−1 mg−1 protein) was observed in interaction of S genotype and no drought condition (control) (Fig. 7).Catalase (CAT) enzyme was affected by drought, genotype and their interaction (Table 2). The results of catalase enzyme activity assessment showed that with increasing drought, catalase activity is different in H, M and S genotypes. The most different reaction in production of CAT was related to H genotype, which with increasing drought stress up to I3 level, the enzyme activity was increased. But regarding M and S genotypes, with increasing drought level, CAT activity was increased in both genotypes. In this study the highest amount of CAT (460.51 μmol min−1 mg−1 protein) was observed in interaction of S genotype with 30–40% of field capacity; and the lowest amount (157.06 μmol min−1 mg−1 protein) was observed in interaction of H genotype with 90–100% of field capacity (Fig. 7).No significant effect was observed for APX enzyme in interaction of genotype and drought (Table 2). The results of comparing average data, suggest that the highest amount of APX activity (284.96 μmol min−1 mg−1 protein) was observed in H genotype (Table 3). Also the results showed that I2, I3 and I4 drought level resulted in an increase in APX enzyme activity by 32.38%, 49.16%, and 65.53% respectively. This indicates that APX enzymatic activity increases with increasing drought level (Table 3).Using physiological and biochemical mechanisms to reduce effects of stress shows that to overcome drought, oxidative stress and to eliminate ROS, plants will increase the amount of antioxidant content55. One of major mechanisms to cope with oxidative stress in plants, is activation of antioxidant enzymes61. Findings of the present study indicates that different lavender genotypes showed partial resistance against drought. In this research, increased activity of antioxidant enzymes in lavender genotypes under drought condition, was considered as an important drought-resistance factor. Among all antioxidant enzymes, SOD can have a good response against drought stress. In a way that H, S, and M genotypes of lavender in the highest level of drought stress (I4), showed an increased amount of SOD, by 57.42%, 35.85% and 60.69% compared to normal conditions (Fig. 7).In this study, the minimum enzymatic changes were related to the POX enzyme and the highest enzymatic changes were related to the CAT enzyme. Moreover, it was observed that the highest amount of catalase enzymatic activity was in H genotype. In a way that in plants under drought stress CAT activity was increased up to I3 drought level; but, after this level with increasing drought (I4 drought level) CAT enzymatic activity was decreased. CAT and POX are among important plants enzymes which can protect plant cells against free radicals63. In this study, in drought period, enzymatic activity of CAT and POX was increased, this means that lavender genotypes, in the face of stress produce antioxidant enzymes to protect themselves. While in H genotype compared to other genotypes, in high drought stress, CAT activity was decreased which this response indicates the different function of this genotype in dealing with ROS. Enzymatic response to drought condition was different in various lavender genotypes. Generally, the negative effect of drought is shown by production of reactive oxygen species (ROS). Increased enzymatic activity of antioxidant enzymes, particularly CAT and POX can reduce the negative effects of drought64, 65. In this regard, increased activity of antioxidant enzymes in different genotypes of Calendula officinalis plant was reported to56.Malondialdehyde (MDA) contentReaction of different lavender genotypes under drought stress was different in terms of malondialdehyde (MDA) production and accumulation (Table 2). With increasing drought, MDA content was significantly increased in M and H genotypes. The highest amount of MDA in these genotypes was 14.34 and 9.50 nmolg − 1 FM respectively, which was observed in drought level of 30–40% of field capacity. This indicates a significant increase in MDA content with increasing drought (Fig. 8). While the process of production and accumulation of MDA in S genotype was different at various drought levels. For S genotype, in first level of drought (I2), MDA content was increased which showed the vulnerability of the cell membrane at this drought level. But with increasing drought, gradually, the S genotype plants adapted to the dry environment, which in this level cell membrane damage was not obvious. Then, increasing in drought stress resulted in increased MDA content. Generally, in I2 and I3 drought levels, lavender genotypes underwent varying degrees of damage, which in M and H genotypes followed by increasing enzymatic activity, and in S genotype it resulted in decreased enzymatic activity. But in the highest level of drought (I4), the cell membrane was seriously damaged and in all three genotypes and MDA content was significantly increased (Fig. 8).Figure 8The effect of drought stress on MDA content in different lavender genotypes.Full size imageMembrane lipid peroxidation due to the accumulation of active oxygen species leads to cell damage and death. In plants this lipid peroxidation happens under drought stress66. MDA is the final product of membrane peroxidation and membrane processes. Simultaneously with peroxidation, the MDA content increases significantly67. So the MDA content can be considered as an indicator of drought-resistance in plants. Among lavender genotypes, in the highest level of drought, MDA content in M genotype was significantly increased compared to others genotypes; whish suggests that M genotype is more vulnerable in comparison with the two other genotypes. An increase in MDA content under drought stress, was reported in Thymus species as well66.Quantity and quality of essential oilMutual interaction between drought stress and genotype had a significant effect on percentage and yield of essential oil in lavender plants (Table 2). Our findings suggested a different essential oil percentage for each genotype in various levels of drought stress. With increasing drought to I3 level, the essential oil percentage was increased in M and H genotypes, but after that with increasing drought to a higher level (I4), essential oil percentage in these genotypes was decreased. While in S genotype, increasing essential oil percentage totally had an upward trend (Fig. 9).Figure 9The effect of drought stress on essential oil percent in different lavender genotypes.Full size imageEvaluation of essential oil percentage in different levels of drought, showed that in I2 drought level, the highest amount of essential oil (0.81%) was observed in H genotype; and in I3 and I4 drought levels, the highest amounts of essential oil were 1.29% and 1.68% respectively, which were observed in S genotype. Moreover, our results suggest that the highest difference in essential oil percentage in the studied genotypes compared to the control, was related to S genotype (Fig. 9). Totally, the highest percentage of essential oil was observed in S genotype in I4 drought level. This shows the high capacity of this genotype to produce essential oil under drought stress.Essential oil yield was significantly affected by genotype and drought. The results showed that the essential oil yield in S genotype was different from the others. So that the highest yield of essential oil (0.055 g per plant) was observed in this genotype in I3 drought level. While in H and M genotypes the highest amounts were 0.068 g and 0.065 g respectively, which were gained in I2 drought level (Fig. 10). Results of comparing average data showed that the highest yield of essential oil at I2 and I3 levels was obtained with 151/85% and 122.22% difference compared to the control, respectively, and they gained from H genotype. This indicates the high potential of H genotype to maintain biomass and produce essential oil in drought stress. Also our results suggest that in the highest drought level (I4), the highest essential oil yield (0.046 g per plant) was observed in M genotype (Fig. 10).Figure 10The effect of drought stress on essential oil yield in different lavender genotypes.Full size imagePrincipal component analysis (PCA)PCA analysis was performed to identify susceptibility of genotypes in irrigation regimes. According to physiological traits in the PCA analysis (Fig. 11a, b), the first factor (PC1) explains about 90% of the total variance of variables, and the second factor (PC2) about 8%.Figure 11Principal component analysis (PCA) for genotypes (a) and physiological traits (b) based on water status calculated for physiological traits. (R 4.0.4 packages, https://rstudio.com/products/rstudio/).Full size imageThe results of PCA analysis of different irrigation regimes showed that in the first component, which shows 89.91% of changes, the best traits are antioxidant enzymes CAT, SOD, APX, while in the second component, with 8.10% changes, only the trait Catalase is the best trait. Also, in total, the first and second components, which show 98.01% of the changes, show CAT as the most effective trait (Fig. 11a).The results of PCA analysis in lavender genotypes showed that the first and second main components could explain 98.91% of the existing changes. So that the first main component with 91.13% and the second component with 7.78% had a share in the total variation. Therefore, using these two components and ignoring other components will only cause the loss of a small part of about 1.09% of the data changes (Fig. 11b). These two principal components include peroxidase, ascorbate peroxidase, and superoxide. Physiological responses of Lavandula genotypes (L. angustifolia cv. Hidcote, L. angustifolia cv. Munstead, and L. stricta) submitted to drought stress were evaluated through principal component analysis (PCA), and the results are illustrated in Fig. 11a. Lavandula stricta presents higher levels of CAT activity than L. angustifolia cv. Hidcote and L. angustifolia cv. Munstead. In addition, APX and CAT increase in stress-treated in 30–40% FC. This result shows that L. stricta exhibits the most affected physiological changes while trying to adjust to changes in the water status of the environment, under the imposed conditions and shows the highest resistance.The results of analysis of essentials oils from H, S and M genotypes is shown is Tables 4, 5 and 6. The trend of changes in essential oils composition is described in all three genotypes. By studying the mass spectra and the Kovats retention index, 23 compounds were identified in the H genotype’s essential oil (Table 4). The yield of H genotype essential oil from I1 to I4 drought levels was 99.89%, 82.78%, 81.09% and 82.85%, respectively. The main components of H genotype essential oil in I1 to I4 drought levels, include 1.8-Cineol compounds (5.94%, 7.73%, 4.24% and 3.50%), Linalool (23.20%, 16.30%, 11.90% and 10.57%), Camphor (3.41%, 4.65%, 2.32% and 2.87%), Borneol (4.89%, 3.34%, 3.65% and 3.01%), Bornyl formate (27.32%, 16.04%, 19.45% and 20.03%), Lavandulyl acetate (1.40%, 4.21%, 6 and 8.35%), Caryophyllene oxide (10.92%, 11.77%, 12.16% and 19.91%), α-Muurolene (4.38%, 3.20%, 1.20% and 0%) (Table 4). The results of grouping the essential oil compounds showed that the amount of hydrocarbon monoterpenes from I1 to I4 drought level were 12.88%, 8.86%, 8.53% and 6.06%, respectively. The amount of oxygen monoterpenes was 64.76%, 50.70%, 43.32% and 42.45%; and hydrocarbon sesquiterpene compounds were 13.12%, 11.45%, 13.03% and 13.96%. The amount of oxygen sesquiterpene compounds were 10.92%, 11.77%, 16.21%, and 19.91%; which shows that increasing drought level, result in decreasing monoterpene compounds, and increasing sesquiterpene compounds.Table 4 Chemical composition of essential oils extracted from Lavandula angustifolia cv. Hidcote plants under different irrigation regime.Full size tableTable 5 Chemical composition of essential oils extracted from Lavandula stricta plants under different irrigation regime.Full size tableTable 6 Chemical composition of essential oils extracted from Lavandula angustifolia cv. Munstead plants different irrigation regime.Full size tableHeat map for the essential oil profile in Lavandula angustifolia cv. Hidcote corresponding to the different irrigation regime The similar discrimination was also supported by the heatmap constructed for essential compounds. Accordingly, 22 rows and 4 columns were achieved. α- pinene, β-Pinene, δ-3-Carene, type of Cymene, 1,8-Cineol, Camphor and Linalool from the main compounds, peaked at control. Moreover, lavandulyl acetate, Myrtenyl acetate, caryophyllene oxide, camphene and γ-Cadinene revealed highest percentage at 30–40% FC, Some compounds, such as Camphor and Linalyl acetate, are at the levels of the intermediate irrigation regime (Fig. 12). It is remarkable that as the water limit increases, the amount of monoterpene compounds decreases and the amount of sesquiterpene compounds increases.Figure 12Heatmap for the essential oil profile in aerial parts of Lavandula angustifolia cv. Hidcote corresponding to irrigation regimes (CIMminer, https://discover.nci.nih.gov/cimminer/oneMatrix.do).Full size imageWith evaluation of the essential oil from S genotype, 18 compounds were identified (Table 5). The amount of essential oil in I1 to I4 drought levels was 99.41%, 98.48%, 99.53% and 99.93% respectively (Table 5). Among identified compounds in S genotype the followings were accounted for the highest amount of components in the essential oil in I1 to I4 levels respectively; Linalool (32.60%, 28.45%, 20.12% and 19.12%), decanal (10.26%, 15.21%, 18.56% and 19.27%), 1-Decanol (8.01%, 10.31%, 17.88% and 21.34%), Kessane (2.44%, 4.43%, 9.99% and 11.50%), Hexadecane (1.26%, 5.77%, 6.10% and 11.9%), 2-methyl-1-hexadecanol (11.1%, 9.32%, 8.15% and 2.37%) and Hexahydrofarnesyl acetone (6.8%, 6.34%, 3.78% and 1.26%) (Table 5). The most obvious point was the high percentage of Linalool, decanal and 1-Decanol in the S genotype. With increasing drought, Linalool compounds were decreased and decanal and 1-Decanol compounds were increased. The grouping of essential oil components also showed that among the 18 compounds identified, the following were the highest in I1 to I4 drought levels, respectively; 3 hydrocarbon monoterpenes with total of (5.34%, 5.44%, 4.57% and 4.34%), 6 oxygen monoterpenes with total of (60.49%, 61.03%, 59.57% and 60.45%), 3 hydrogen sesquiterpenes with total of (5.69%, 10.27%, 11.85% and 15.24%) and 6 oxygen sesquiterpenes with total of (27.89%, 28.09%, 29.44% and 18.32%). With increasing drought, the amounts of hydrocarbon monoterpenes and oxygen sesquiterpenes were decreased; while the amount of hydrocarbon sesquiterpenes was increased. Also the highest amount of oxygen monoterpenes, by 61.03%, was seen in I2 drought level.Heat map for the essential oil profile in Lavandula stricta corresponding to the different irrigation regime The parallel discrimination was also supported by the heatmap constructed for essential compounds. Accordingly, 18 rows and 4 columns were achieved. α- pinene, Amyl isovalerate, Citronellol, β-Ionone and Linalool from the main compounds, peaked at control. Moreover, α-Thujene, decanal, 1-Decanol, Sesquiphellandrene, Kessane and Hexadecane revealed highest percentage at 30–40% FC (Fig. 13). These results confirm the results obtained from the Lavandula angustifolia cv. Hidcote so that as the water limit increases, the amount of monoterpene compounds decreases and the amount of Sesquiterpene compounds increases.Figure 13Heatmap for the essential oil profile in aerial parts of Lavandula stricta corresponding to irrigation regimes (CIMminer, https://discover.nci.nih.gov/cimminer/oneMatrix.do).Full size imageEssential oil yield in M genotype from I1 to I4 drought levels was obtained 99.90%, 98.38%, 93.08% and 87.04% (Table 6). As it is shown in Table 6, analysis of the essential oil from M genotype included 27 compounds which its major part was consisted of Camphor (16.82%, 16.32%, 17.11% and 18.30%), Borneol (44.96%, 42.80%, 37.54% and 30.99%) and Caryophyllene oxide (14.68%, 15.21%, 15.90% and 17.21%) from I1 to I4 drought levels, respectively. comparison of essential oil components (Table 6) showed that from 27 identified compounds in M genotype, the followings were the most prevalent from I1 to I4 levels respectively, including hydrocarbon monoterpene with total of (17.82%, 17.45%, 13.91% and 9.96%), 12 total oxygen monoterpene compounds with total of (65.95%, 62.05%, 56.96% and 50.42%), 4 hydrocarbon sesquiterpenes with total of (1.58%, 23.23%, 5.42% and 8.09%) and 2 oxygen sesquiterpenes with total of (14.91%, 15.65%, 16.79% and 18.37%). The highest drought level resulted in 31.76% and 17.23% increase in Camphor and Caryophyllene oxide. It also caused 31.07% decrease in Borneol compared to the control (Table 6). Totally, with increasing drought level, monoterpene compounds were decreased and sesquiterpene compounds were increased in lavender genotypes.The major components of essential oil were different in various lavender genotypes in the highest level of drought (I4). In this study in H genotype, the compounds Linalool, Bornyl formate and Caryophyllene oxide; in S genotype the compounds Linalool, decanal, 1-Decanol, Kessane and Hexadecane; and in M genotype the compounds Camphor, Borneol and Caryophyllene oxide, were the most prevalent components of essential oil. In this study, Borneol compound was not observed in S genotype. regarding the fact that essential oil extraction was performed on flowering branches in all three genotypes, and they were studied under similar drought conditions; and also comparing the results of this study with finding of other studies shows that the difference in types and percentage of essential oil’s components can be due to the effect of genetic differences; and to some extent, environmental factors on essential oil in different genotypes.A total comparison of essential oil analysis results for different lavender genotypes under drought stress showed that oxygen monoterpenes are the most prevalent components of the essential oil, which will decrease with increasing drought level. Sarker et al. (2012) reported that the essential oil of lavender (Lavandula angustifolia) contains high amounts of linalool and linalool acetate, along with scares amount of other monoterpenes68. A study by Hassan et al. (2014) showed that the compounds carvacrol, phenol-2-amino-4, 6-bis, trans-2-caren-4-ol, and n-hexadecanoic acid are the main constituents of Lavandula stricta plants which were collected from the Shaza Mountains in southern Saudi Arabia69. Total results from essential oil analysis in this study showed that Linalool was the main ingredient of essential oils in H and S genotypes. This compound is an oxygen monoterpene with a density of 0.85 and a pleasant smell, and is the main component of the essential oil from lavender plant. While in M genotype, Borneol was the main component of the essential oil, which is a circular monoterpene compound with density of Mohammadnejad ganji et al. (2017) suggested that the difference in natural quality of the essential oil from lavender plants is related to intrinsic factors (genetic or heredity capabilities and maturity), and external factors including sunlight, water, heat, pressure, latitude, and soil which affect plant growth and essential oil production70.Heat map for the essential oil profile in Lavandula angustifolia cv. Munstead corresponding to the different irrigation regime The parallel discrimination was also supported by the heatmap constructed for essential compounds. Accordingly, 18 rows and 4 columns were achieved. α- pinene, Tricycle, Camphene, Thuja-2,4(10)-diene, δ-3-Carene, ρ-Cymene, Borneol and limonene from the main compounds, peaked at control. Moreover, Camphor, α-Santalene, γ-Cadinene, δ-Cadinene, Caryophyllene oxide, α-Muurolene and Ledene oxide-(II) revealed highest percentage at 30–40% FC (Fig. 14). The results showed that the composition of the compounds was similar to the previous two genotypes and the water limit increases, the amount of monoterpene compounds decreases and the amount of Sesquiterpene compounds increases.Figure 14Heatmap for the essential oil profile in aerial parts of Lavandula angustifolia cv. Munstead corresponding to irrigation regimes (CIMminer, https://discover.nci.nih.gov/cimminer/oneMatrix.do).Full size imageEssential oils are generally in the group of terpenoids and The structure of terpenoids consists of two main precursors, isopentenyl pyrophosphate (IPP) and its isomer, dimethylallyl pyrophosphate (DMAPP). These compounds are synthesized via the cytosolic pathway of mevalonic acid (MVA) or plasticity of methylerythritol phosphate (MEP)71. The MVA pathway is primarily responsible for the synthesis of Sesquiterpenoids and triterpenoids, while the MEP pathway is used for the biosynthesis of monoterpenoids, diterpenoids and tetraterpenoids72. Monoterpenes and Sesquiterpenes are the main constituents of essential oils that play a role in aroma, flavor, photosynthetic pigments and antioxidant activities73.In drought conditions, the amount of these isoprenes does not decrease in relation to the mediators of the MEP pathway and in contrast sometimes increases. Therefore, sesquiterpene compounds increase in drought conditions because most of these compounds are synthesized through the MVA pathway74. Another reason for the decrease in MEP path flux is the location of this path, which has a significant impact in drought conditions. In this case, plastids are not able to provide the required IPP of this path, so most monoterpene compounds are reduced75.Also, since the quality of the essential oil is due to the presence of linalool and linalyl acetate76. According to the results obtained from heatmaps related to essential oils, three genotypes are identified, the highest amount of linalool amount in S genotype was remained under mind- (I2) till severe-drought (I4) condition. This indicates more compatibility with maintaining the desired quality of drought conditions in this plant than the other two commercial genotypes. And then the H genotype is in the second stage due to the presence of important compounds.Comparing the grouping created in the heat maps related to the essential oil of 3 genotypes, it is clear that the two genotypes S and H were divided into two groups I1, I2 and I3, I4 in the genotype. But in the genotype M, the results were divided into I4 and I3 groups I2 were divided into genotypes. This can be due to differences in the resistance mechanism of plants in different genotypes, so in genotypes S and H of the plant through increasing sesquiterpene compounds showed resistance to drought stress, while in genotype M increased resistance to drought levels through higher monoterpene compounds. Another conclusion that can be drawn from these heat maps is that in genotypes S and H, the rate of drought resistance in the first and second levels of drought with the third and fourth levels has shown more changes in the type of essential oil compounds, while in the third genotype (M) these changes in the last level drought has been most evident.At a glance, it seems Genotype S has a different mechanism in reducing the negative effects of drought compared to genotypes M and H, So that, among the enzymatic and non-enzymatic mechanisms, it tends to use the enzymatic pathway more. In association with the production of “proline”, drought stress index osmolyte, genotype S has a different trend from genotypes H and M and this osmolyte in this plant has a lower production flux compared to other genotypes. Also, due to the fact that the production of soluble sugars in this plant has been moderate compared to other genotypes, it is expected this genotype replace proline with another osmolyte or uses an enzymatic mechanism to deal with drought, as the results of antioxidant enzyme “catalase” related to genotype S had the highest value with a significant difference under drought stress, while, in the H and M genotypes, the SOD enzyme was responsive to drought.On the other hand, the high resistance of genotype S can be attributed to the greater activation of the pathway of essential oil compounds. Because by examining the constituents of the essential oil (monoterpene and sesquiterpene), it can be concluded that genotype H and then M at high drought levels still retain the ability to produce monoterpene compounds, while in genotype S with increasing drought, the amount of semi-heavy compounds (sesquiterpene) has increased significantly (Fig. 15), this can confirm the existence of a different resistance mechanism in the S genotype. Because some structural compounds of the membrane, such as sterols, are made from the mevalonic acid (MVA) pathway of acetyl coenzyme A origin. For this reason it seems that S genotype by setting up terpenoid pathways involved in the production of steroids another solution to drought is by preserving its plasma membrane. Steroids are derivatives of triterpenes that, along with phospholipids, are major components of plasma membranes70. Also, the study of MDA content as the final product of membrane lipid peroxidation in genotypes at the fourth level of drought (the most severe drought) showed the M genotype is most sensitive to drought. In this way, the two genotypes S and H have almost equal MDA content, so that it can be said that with a small difference from genotype S, genotype H has less composition.Figure 15The amount of monoterpenes and sesquiterpene compounds in different genotypes under irrigation regimes.Full size imageContinuous production of isoprene under drought conditions shows that despite the reduction in the synthesis of osmolyte and relative increasing of MDA (with very little difference from genotype H) that occurs under these conditions, the function of this pathway is essential for the S genotype. Isoprene has long been used to protect plants from drought, high temperatures and oxidative stress are recommended77. Of course, it was showed which is possible with increasing drought, sufficient isoprene is not produced to counteract and launch defense pathways and instead used as a general signal to increase drought tolerance78,79.Reasons such as further activation of terpenoid skeletal pathways towards the production of semi-heavy (sesquiterpene) compounds, production of steroids via the MVA pathway could be a reason for lower susceptibility of S genotype and high resistance of this genotype through these mechanisms compared to other genotypes. In contrast, on the one hand, H genotype using proline production, soluble sugar levels and decreased MDA in response to stress caused by drought and on the other hand, the ability to produce substances important monoterpenes, such as Linalool and Linalyl acetate, with the aim of using medicine and aromatherapy76, It (H genotype) can be considered as a cultivar with high commercial value and significant resistance to M genotype. More