Philippa Kaur

Nash, R. D. M., Valencia, A. H. & Geffen, A. J. The origin of Fulton’s condition factor—setting the record straight. Fisheries 31(5), 236–238 (2006).
Google Scholar 
Tarkan, A. S., Gaygusuz, Ö., Acıpınar, H., Gürsoy, Ç. & Özuluğ, M. Length–weight relationships of fishes from the Marmara region (NW-Turkey). J. Appl. Ichthyol. 22(4), 271–273 (2006).Article 

Google Scholar 
Al Nahdi, A., de Leaniz, C. G. & King, A. J. Spatio-temporal variation in length-weight relationships and condition of ribbonfish Trichiurus lepturus (Linnaeus, 1758): Implications for fisheries. PLoS One 11(8), e0161989 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Froese, R., Tsikliras, A. C. & Stergiou, K. I. Editorial note on weight–length relations of fishes. Acta Ichthyol. Piscat. 41(4), 261–263 (2011).Article 

Google Scholar 
Torres, M. A. et al. Length–weight relationships for 22 crustecans and cephalopods from the Gulf of Cadiz (SW Spain). Aquat. Liv. Resour. 30, 12 (2017).Article 

Google Scholar 
Rocha, F., Guerra, A. & Gonzalez, A. F. A review of reproductive strategies in cephalopods. Biol. Rev. 76, 291–304 (2001).CAS 
PubMed 
Article 

Google Scholar 
Laptikhovsky, V. & Salman, A. On reproductive strategies of the epipelagic octopods of the superfamily Argonautoidea (Cephalopoda: Octopoda). Mar. Biol. 142, 321–326 (2003).Article 

Google Scholar 
Forsythe, J. W. & van Heukelem, W. F. Growth. In Cephalopod Life Cycles (ed. Boyle, P. R.) 135–156 (Academic Press, 1987).
Google Scholar 
Jereb, P., et al. (eds) 2015. Cephalopod biology and fisheries in Europe: II. Species Accounts. ICES Cooperative Research Report No. 325, p. 360.Salman, A. Cephalopod research in the eastern Mediterranean (East of 23°E): A review. Boll. Malacol. 45, 47–59 (2009).
Google Scholar 
Salman, A. & İzmirli, C. Ege Üniversitesi Su Ürünleri Fakültesi Müzesi (ESFM)’nin cephalopod envanteri. EgeJFAS 37(4), 357–361. https://doi.org/10.12714/egejfas.37.4.06) (2020) (in Turkish).Article 

Google Scholar 
Önsoy, B. & Salman, A. Reproductive biology of the common cuttlefish Sepia officinalis L. (Sepiida: Cephalpoda) in the Aegean Sea. Turk. J. Vet. Anim. Sci. 29, 613–619 (2005).
Google Scholar 
Salman, A. Reproductive biology of the elegant cuttlefish (Sepia elegans) in the Eastern Mediterranean. Turk. J. Fish. Aquat. Sci. 15(2), 265–272 (2015).Article 

Google Scholar 
Dursun, D., Eronat, E. G. T., Akalın, M. & Salman, M. A. Reproductive biology of pink cuttlefish Sepia orbignyana in the Aegean Sea (eastern Mediterranean). Turk. J. Zool. 37, 576–581 (2013).Article 

Google Scholar 
Salman, A. Reproductive biology of Sepietta oweniana (Pfeffer, 1908) (Sepiolidae: Cephalopoda) in the Aegean Sea. Sci. Mar. 62(4), 379–383 (1998).Article 

Google Scholar 
Salman, A. & Önsoy, B. Reproductive biology of the bobtail squid Rossia macrosoma (Cephalopoda: Sepiolidea) from the eastern Mediterranean. Turk. J. Fish. Aquat. Sci. 10, 81–86 (2010).Article 

Google Scholar 
Salman, A. Fecundity and spawning strategy of shortfin squid Illex coindetii (Oegopsida: Ommastrephidae), in the eastern Mediterranean. Turk. J. Fish. Aquat. Sci. 17, 841–849 (2017).
Google Scholar 
Mangold-Wirz, K. Biologie des céphalopodes benthiques et nectoniques de la Mer Catalane. Vie Millieu suppl. 13, 1–285 (1963).
Google Scholar 
Salman, A. Fecundity, spawning strategy and oocyte development of shortfin squid Alloteuthis media (Myopsida: Loliginidae) in the eastern Mediterranean. Cah. Biol. Mar. 55, 163–171 (2014).
Google Scholar 
Önsoy, B. & Salman, A. Reproduction patterns of the Mediterranean endemic, Eledone moschata (Lamarck, 1798) (Octopoda: Cephalopoda) in the eastern Mediterranean. (In Turkish) 1st National Malacology Congress, 1–3 September 2004, Izmir-Turkey (Bilal Öztürk & Alp Salman, eds). Turk. J. Aquat. Life 2(2), 55–60 (2004).
Google Scholar 
Tesch, F. W. Age and growth. In Methods for Assessment of Fish Production in Fresh Waters (ed. Ricker, W. E.) 99–130 (Blackwell Scientific Publications, 1971).
Google Scholar 
Merella, P., Quetglas, A., Alemany, F. & Carbonell, A. Length–weight relationship of fishes and cephalopods from the Balearic Islands (western Mediterranean). Naga ICLARM Q. 20(3–4), 66–68 (1997).
Google Scholar 
Manfrin Piccinetti, G. & Giovanardi, O. Données sur la biologie de Sepia officinalis L. dans l’Adriatique obtenues lors de expéditions pipeta. FAO Fish. Rep. 290, 135–138 (1984).
Google Scholar 
Bello, G. Length–weight relationship in males and females of Sepia orbignyana and Sepia elegans (Cephalopoda: Sepiidae). Rapp. Comm. Int. Mer. Médit. 31(2), 254 (1988).
Google Scholar 
Ragonese, S. & Jereb, P. Length-weight relationship and growth of the pink and elegant cuttlefish Sepia orbignyana and Sepia elegans in the Sicilian Channel. In Acta of the 1st International Symposium on the Cuttlefish (ed. Boucaud-Camou, E.) 31–47 (SEPIA. Centre de Publications de l’Universite de Caen, 1991).
Google Scholar 
Akyol, O. & Metin, G. An investigation on determination of some morphological characteristics of Cephalopods in Izmir Bay (Aegean Sea). EU J. Fish. Aquat. Sci. 18(3–4), 357–365 (2001).
Google Scholar 
Lefkaditou, E., Verriopoulos, G. & Valavanis, V. VII9. Research on Cephalopod resources in Hellas. In State of Hellenic Fisheries (eds Papaconstantinou, C. et al.) 440–451 (HCMR Publications, 2007).
Google Scholar 
Duysak, Ö., Sendão, J., Borges, T., Türeli, C. & Erdem, Ü. Cephalopod distribution in Iskenderun bay (eastern Mediterranean–Turkey). J. Fish. Sci. 2, 118–125 (2008).
Google Scholar 
Giordano, D. et al. Distribution and biology of Sepietta oweniana (Pfeffer, 1908) (Cephalopoda: Sepiolidae) in the southern Tyrrhenian Sea (central Mediterranean Sea). Cah. Biol. Mar. 50, 1–10 (2009).
Google Scholar 
Andriguetto, J. M. Jr. & Haimovici, M. Effects of fixation and preservation methods on the morphology of a Loliginid squid (Cephalopoda: Myopsida). Am. Malac Bull. 6(2), 213–217 (1988).
Google Scholar 
Sanchez, P. Donnés preliminaires sur la biologie de trois species de cephalopods de la Mer Catalan. Rapp. Comm. Int. Mer. Médit. 30(2), 247 (1986).
Google Scholar 
Belcari, P., Sartor, P., Nannini, N. & De Ranieri, S. Length-weight relationship of Toda- ropsis eblanae (Cephalopoda: Ommastrephidae) of the northern Tyrrhenian Sea in relation to sexual maturation. Biol. Mar. Mediter. 6, 524–528 (1999).
Google Scholar 
Belcari, P. Length–weight relationship in relation to sexual maturation of Illex coindetii (Cephalopoda: Ommastrephidae) in the northern Tyrrhenian Sea (western Mediterranean). Sci. Mar. 60, 379–384 (1996).
Google Scholar 
Petric, M., Ferri, J., Skeljo, F. & Krstulovic Sifner, S. Body and beak measures of Illex coindetii (Cephalopoda: Ommastrephidae) and their relation to growth and maturity. Cah. Biol. Mar. 51, 275–287 (2010).
Google Scholar 
Ceriola, L., Ungaro, N. & Toteda, F. Some information on the biology of Illex coindetii Verany, 1839 (Cephalopoda, Ommastrephidae) in the south-western Adriatic Sea (central Mediterranean). Fish. Res. 82, 41–49 (2006).Article 

Google Scholar 
Arvanitidis, C. et al. A comparison of the fishery biology of three Illex coindetii Verany, 1839 (Cephalopoda: Ommastrephidae) populations from the European Atlantic and Mediterranean Waters. Bull. Mar. Sci. 71, 129–146 (2002).
Google Scholar 
Quetglas, A., Alemany, F., Carbonell, A., Merella, P. & Sanchez, P. Some aspects of the biology of Todarodes sagittatus (Cephalopoda: Ommastrephidae) from the Balearic Sea (western Mediterranean). Sci. Mar. 62, 73–82 (1998).Article 

Google Scholar 
Krstulovic Sifner, S. K. & Vrgoc, N. Population structure, maturation and reproduction of the European squid, Loligo vulgaris, in the central Adriatic Sea. Fish. Res. 69, 239–249 (2004).Article 

Google Scholar 
Moreno, A. et al. Biological variation of Loligo vulgaris (Cephalopoda: Loliginidae) in the eastern Atlantic and Mediterranean. Bull. Mar. Sci. 71(1), 515–534 (2002).
Google Scholar 
Guerra, A. & Manriquez, M. Parametros biometricos de Octopus vulgaris. Invest. Pesq. 44, 177–198 (1980).
Google Scholar 
Quetglas, A., Alemany, F., Carbonell, A., Merella, P. & Sanchez, P. Biology and fishery of Octopus vulgaris Cuvier, 1797, caught by trawlers in Mallorca (Balearic Sea, western Mediterranean). Fish. Res. 36, 237–249 (1998).Article 

Google Scholar 
Sanchez, P., & Obarti, R. 1993. The biology and fishery of Octopus vulgaris caught with clay pots on the Spanish Mediterranean coast. In: Jereb, P., Allcock, A. L., Lefkaditou, E., Piatkowski, U., Hastie, L. C., Pierce, G. J. (Eds.) 2015. Cephalopod biology and fisheries in Europe: II. Species Accounts. ICES Cooperative Research Report No. 325, p 360.Gonzalez, M., Barcala, E., Perez-Gil, J. L., Carrasco, M. N. & Garcia-Martinez, M. C. Fisheries and reproductive biology of Octopus vulgaris (Mollusca: Cephalopoda) in the Gulf of Alicante (Northwestern Mediterranean). Medit. Mar. Sci. 12, 369–389 (2011).Article 

Google Scholar 
Jabeur, C., Nouira, T., Khoufi, W., Mosbahi, D. S. & Ezzedine-Najai, S. Age and growth of Octopus vulgaris Cuvier, 1797 along the east coast of Tunisia. J. Shellf. Res. 31, 119–124 (2012).Article 

Google Scholar 
Quetglas, A., Ordines, F., Gonzalez, M. & Franco, I. Life history of the bathyal octopus Pteroctopus tetracirrhus (Mollusca, Cephalopoda) in the Mediterranean Sea. Deep Sea Res. Part I 56, 1379–1390 (2009).Article 

Google Scholar 
Quetglas, A., Gonzalez, M. & Franco, I. Biology of the upper-slope cephalopod Octopus salutii from the western Mediterranean Sea. Mar. Biol. 146, 1131–1138 (2005).Article 

Google Scholar 
Moriyasu, M. Etude biometrique de la croissance d’E. cirrhosa [LAM. 1798 (Cephalopoda, Octopoda)] du Golfe du Lion. Oceanol. Acta 6, 35–41 (1983).
Google Scholar 
Massi, D. Effetti del congelamento sull’accuratezza delle misure in Eledone cirrhosa (Lamarck, 1798). Biol. Mar. Suppl. al Notiziario S.I.B.M. 1, 379–380 (1993).
Google Scholar 
Agnesi, S., Belluscio, A. & Ardizzone, G. D. Biologia e dinamica di populazione di Eledone cirrhosa (Cephalopoda: Octopoda) nel Tirreno Centrale. Biol. Mar. Mediterr. 5, 336–348 (1998).
Google Scholar 
Giordano, D. et al. Population dynamics and distribution of Eledone cirrhosa (Lamarck, 1798) in the Southern Tyrrhenian Sea (Central Mediterranean). Cah. Biol. Mar. 51, 213–227 (2010).
Google Scholar 
Krstulovic Sifner, S. K. & Vrgoc, N. Reproductive cycle and sexual maturation of the musky octopus Eledone moschata (Cephalopoda: Octopodidae) in the northern and central Adriatic Sea. Sci. Mar. 73, 439–447 (2009).Article 

Google Scholar 
Ikica, Z., Krstulovic Sifner, S. & Joksimovic, A. Some preliminary data on biological aspects of the musky octopus, Eledone moschata (Lamarck, 1798) (Cephalopoda: Octopodidae) in Montenegrin waters. Stud. Mar. 25, 21–36 (2011).
Google Scholar 
Akyol, O., Şen, H. & Kinacigil, H. T. Reproductive biology of Eledone moschata (Cephalopoda: Octopodidae) in the Aegean Sea (Izmir Bay, Turkey). J. Mar. Biol. Assoc. UK 87, 967–970 (2007).Article 

Google Scholar 
Quetglas, A., Gonzalez, M., Carbonell, A. & Sanchez, P. Biology of the deep-sea octopus Bathypolypus sponsalis (Cephalopoda: Octopodidae) from the western Mediterranean Sea. Mar. Biol. 138, 785–792 (2001).Article 

Google Scholar  More

For a schematic overview of the experimental design, see Fig. 2.Experimental evolutionProtocolThe stock population (Stock population below) was allowed to expand for one generation and from this we established eight replicate experimental evolution populations, four selected for fighter morphs (F-lines) and four selected for scrambler morphs (S-lines). Each population was founded by 1,000 recently eclosed adults (500 random females and 500 random males of the desired morph). The classification of the morphs was based on visual inspection using a stereoscopic microscope and was unambiguous due to the discontinuous distribution of the phenotypes (Classifying male morphs below). Adults were allowed to interact freely for 6 days, all surviving adults (with previously laid eggs discarded) were transferred to a new container for 24 h of egg laying, after which adults were removed. The resulting offspring were allowed to mature over 13 days and 1,000 individuals from the newly eclosed adults selected for founding the following generation, again 500 random females and 500 random males of the desired morph, with this protocol repeated every generation (Extended Data Fig. 1). The isolation of nymphs to use virgins was unfeasible with our experimental design and population sizes. However, the period of 6 days after selecting the founders of the next generation and collecting eggs for the next generation was probably enough to displace most sperm stored by females mated with any unselected males due to the high number of remating that will be occurring over this duration (females on average remate after 80 min, ref. 88) and last male sperm precedence89. The timing of generation was chosen to reflect maturation rates from our stock population to avoid indirect selection on this trait. Moreover, a previous study90 showed there was no difference between male morphs in maturation rates and that over similar lengths of time to the protocol here the fertility of both morphs remains similar. Therefore, our protocol was not likely to impose strong differential selection on morph life histories.Tracking morph proportionWe assayed the proportion of male morph in each population every 6–7 generations, by isolating 200 larvae (ten per vial) from the container, allowing maturation within vials and recording the morph of all males that eclosed (mean n = 86 per population, per generation, range 71–109). Our selection protocol was highly effective in driving an increase in the frequency of the desired male morph to >90% after 20 generations in both treatments, with this effect considerably faster within F-lines indicted by a significant two-way interaction between proportion of the desired morph and generation (χ2 = 39.9, d.f. = 6, P 90% is probably a consequence of a longer interaction period (3 versus 6 days) in which the stored sperm of males before selection was able to be displaced and/or because selection was acting more efficiently in our larger populations. The difference between rate of changes in morph proportion between F- and S-lines in the current study, and also found by Plesnar-Bielak et al.39, may be associated with the genetic architecture of morph expression. Alternatively, selection could be less effective in scrambler lines if they are less efficient than fighters in displacing sperm of previous females’ partners, but this is unlikely as R. robini male morphs have previously been demonstrated to not differ in their sperm competiveness89.Stock populationWe established a stock population by mixing three laboratory populations that were collected from three sites in Poland (Krakόw, collected in 1998 and 2008, Kwiejce, collected in 2017 and Mosina, collected in 2017; Extended Data Fig. 1), where the line derived for material used in creating the reference genome (below) was also established from the same collections at Mosina in 2017. All populations were maintained in cultures with several hundred individuals per generation before mixing and establishment of the stock population. The mixing of distinct populations increased the genetic variance in the stock population, which otherwise would probably have been limited due to founder events and the limited population size of each of the contributing populations73, thus decreasing our power to detect the effects of SSTs on genetic variation. The newly mixed stock population was maintained with several hundred individuals per generation for roughly 12 generations before the onset of this experiment. This time period is probably enough to break linkage disequilibria that could have arisen due to mixing (for unlinked loci, linkage disequilibrium should decay by half each generation91).One generation before establishing experimental evolution populations the proportion of male morphs was determined from 176 random males, indicating a roughly equal morph ratio (95 fighters, 81 scramblers) of the stock population (Extended Data Fig. 3).General housing and husbandryThe stock population and experimental evolution populations were maintained in plastic containers (approximate length, 9.5 cm; width, 7 cm; height, 4.5 cm), filled with roughly 1 cm of plaster-of-Paris. The same containers were used when sampling mites for sequencing for the reference genome or resequencing from experimental evolution populations, but either replaced the plaster-of-Paris with 5% agarose gel or added a thin layer of 5% agarose gel above the plaster-of-Paris, respectively. The agarose gel was used to reduce the number of contaminates within our samples and on the basis of preliminary extractions that indicated that small pieces of plaster-of-Paris may reduce the quality of DNA during extractions. Individuals, pairs and small groups of ten mites were housed in glass vials (approximate height, 2 cm; diameter, 0.8 cm) and large groups of 60 or 150 mites in plastic containers (approximate height, 1.5 cm; diameter, 2 cm diameter or height, 1.5 cm; diameter, 3.5 cm diameter, respectively) all with an approximate 1 cm base of plaster-of-Paris. All plaster-of-Paris bases were completely soaked in water before mites were transferred into them. All mites were reared at a constant 23 °C, at high humidity ( >90%) and were provided an excess of powdered yeast ad libitum.Classifying male morphsTo illustrate the discontinuous distribution of the weapon and to demonstrate that this classification based on visual inspection is non-subjective, we performed phenotypic measurements from male mites from a population collected near Krakόw, Poland, that had previously been fixed onto microscope slides for a separate study66. The measurements taken were idiosoma (body without mouthparts) length and width of third proximal segment of the third right leg (genu). Measurements were preformed using Lecia DM5500B microscope and Lecia Application Suite v.4.6.1. We then performed an analysis to, first, determine whether the allometric relationship between idiosoma length and width of third pair of legs is best described as discontinuous and, second, to verify that classification by simple visual inspection matches the same classification from allometric analysis. One researcher performed all the measurements and classified each male as a fighter (n = 50) or scrambler (n = 50), a separate researcher was then given the measurements but not the classification of the male morph.Broadly, guidelines for the analysis of non-linear allometries92 were followed. The log–log scatterplots of idiosoma length against leg width were visualized, which showed there was clear evidence for non-linear scaling relationships. Next histograms of idiosoma length, leg width and relative leg width (leg width/idiosoma length) were visualized (Extended Data Fig. 2a–c). Where a normal distribution of idiosoma length, and a binomial distribution in leg width and relative leg width are further indications of a discontinuous relationship. On the basis of the lowest point between the two peaks of the density plot of relative leg width (Extended Data Fig. 2c) males were classified as scramblers (relative leg width 0.125). Replotting the log–log scatterplot of idiosoma length and leg width, and using the classification of morph described above clearly demonstrates the discontinuous allometric relationship of idiosoma length and leg width in R. robini (Extended Data Fig. 2d). Moreover, on the basis of the Akaike information criterion (AIC), the discontinuous model where males were assigned a morph (AIC = 646.5) clearly has a substantially better fit than a simple linear and quadratic models (AIC = 918.5 and 920.2, respectively). Further models were omitted from comparison (for example, breakpoint or sigmoidal) due to the clear discontinuous allometry observed. Finally, all 100 males were assigned the same morph by visual inspection and blind allometric analysis, demonstrating that the former is effective and accurate in classifying male morph.Phenotypic assaysFecundity assays were performed using experimental evolution females at F20 and F32. Eggs laid by females between days 4–8 were counted, encompassing the window of time of most evolutionary relevance for female fitness during maintenance of selection lines (that is, egg laying period in selection lines was between days 6–7) and also likely to capture variation in lifetime fecundity that remains largely consistent throughout the first 3 weeks of life93. Nymphs were individually isolated to gain virgin females, which on maturation females from each experimental evolution population (n = 30) were paired with a male from the stock population (15 with fighters and 15 with scramblers). Pairs were transferred to a new vial on day 4, with the pair being removed from the second vial after a further 4 days and all eggs in the second vial counted. If the male had died in the first vial, they were replaced with a stock male of the same morph. Any female deaths in the first or second vials were recorded.Longevity assays were also performed at F20 and F32. At F20, females used in fecundity assays, including the stock male they were paired with (replaced if dead), were transferred to a new vial at day 8. After this point, vials were then checked every 2 days for female deaths and pairs were moved to new vials every 4 days. Males were replaced with stock males of the same morph if found dead. Similarly, at F20, on maturation males from experimental evolution populations (n = 30) were paired with stock females, vials were checked every 2 days and changed every 4 days, with females being replaced if dead. At F32, only female longevity was determined and was performed in groups; 30 experimental evolution females and 30 stock males (15 of each morph) were placed in plastic containers, two per experimental population. This logistically easier estimate of longevity was done due to local restrictions during the SARS-CoV-2 pandemic and the imposed limitations on people working closely together. Groups were checked for dead females every other day and all remaining live mites transferred to a new container every 4 days. When mites were transferred to a new container the sex and morph ratio were balanced to that of the remaining females, by either removing or adding males of the desired morph from the stock population.To determine whether the survival of mites differed between F- and S-lines when competition between males was allowed, at F45 we created small colonies from each population and survival of males and females recorded over 6 days, the same period as used between selecting founders of the next generation and subsequent egg laying period. Colonies were at a 50:50 sex ratio, established with 150 newly eclosed mites placed into small plastic containers. This was approximately the same density after selection of the next generations founders during the maintenance of experimental evolution populations (150 mites in roughly 9.5 cm2 = 16 mites per 1 cm2; 1,000 mites in roughly 67 cm2 = 15 mites per 1 cm2). After 3 days, all colonies were checked and any dead mites identified by sex. After another 3 days, again dead mites were recorded and all surviving mites sexed and counted.Additionally, at F45 we performed further fecundity assays to obtain estimates of inbreeding depression within experimental evolution populations. To establish family groups, larvae were isolated and on maturation F0 males and females (n = 16) from within the same experimental evolution population were paired together. Pairs were allowed to produce eggs for 48 h, after which adults were removed from vials. After hatching from each pair, 12 F1 larvae were isolated into new vials. On their maturation, these F1 mites were either paired with a full sibling, that is, from the same family, or with an individual from a different family but from the same experimental evolution population. When possible, we made two inbred and two outbred pairs with same family lines used. Again, pairs were allowed to produce eggs for 48 h before their removal for the vial. After a further 5 days, vials were checked for larvae, if larvae were present in the first vial six were individually isolated and the second vial discarded, if no larvae were present in the first vial the second vial was checked for larvae and, if present, they were isolated. This protocol therefore produced inbred and outbred individuals from within the same experimental evolution population. Which, as above, on maturation F2 inbred and outbred females were paired with stock males (fighter males only) and number of eggs laid between days 4 and 8 counted. Only a single female from each unique inbred or outbred family was used. Either due to pairs failing to produce offspring or there being no F2 females, samples sizes were not exactly equal. In total, 59 outbred and 55 inbred females from F-lines, and 56 outbred and 54 inbred females from S-lines were paired with stock males.Phenotypic assay statistical analysesAll phenotypic analysis was conducted using R statistical software94 (v.3.5.2) and data were visualized using ggplot2 (ref. 95).Analysis of male morph proportion was performed using a generalized linear mixed model with binomial error structure, fitted using lme4 (ref. 96). Where the proportion of desired morph was compared in model with morph selection and generation (as a factor) including their two-way interaction as explanatory variables, and population included as a random effect.All fecundity data were analysed using generalized linear mixed models with Poisson error structures, fitted using lme4. Due to the differences in stock population males used between F45 and earlier generations, and slightly different rearing conditions between females in the fecundity assays from generations F20 and F32, they were analysed separately from data collected in F45. However, we noted that the fecundity of females in Fig. 5a was comparable to the outbred females in Fig. 5b. Explanatory variables fitted to fecundity data from F20 and F32 were, morph selection treatment, generation, including their two-way interaction term, and stock male morph. The explanatory variables fitted to fecundity data from inbreeding depression data were, morph selection treatment and status of female (that is, inbred or outbred), including their two-way interaction term. In both analyses, we included population as a random effect and an observation level random effect to account for overdispersion, we omitted fitting random slopes due to issues with increasing the complexity of random effects close to reaching a singular fit. Females that died before the end of the fecundity assay and those that laid zero eggs were removed from analysis. This excluded five females from F20 (three F-line and two S-line), 20 from F32 (13 F-line and seven S-line) and 16 from F45 (three inbred and three outbred F-line, and nine inbred and one outbred S-line).Longevities of females at F20 and F32, and males at F20, were analysed separately using mixed effects Cox models, fitted using coxme97. In all analyses, we included a random effect of population, with morph selection treatment as an explanatory variable and extra variable of male morph included in female longevity analysis at F20. Survival of mites over 6 days at F45 was analysed using a GLM with counts of dead and surviving mites fitted with a quasibinomial error structure, the model included morph selection treatment and sex, including their interaction term, as explanatory variables. If individuals were lost due to handling error (that is, killed or escaped) they were right-censored during analysis.Genome assemblySample originA line of R. robini originated from a wild-collected population from the Mosina region (Wielkopolska, Poland). In October 2017, onions were collected from the field and approximately 200 individuals of R. robini were identified under dissecting microscope. The line used for DNA isolation in the genome sequencing project was developed from full sib × sib mating for 14 generations (to maximize homozygosity) following and continuing the protocol described in ref. 67.DNA extractionFor DNA extraction we used only mite eggs, that were laid by 500 females, collected in a container (see above for a description) Females were kept in this container for 3 days. After that time, they were removed, and eggs were filtered using fine sieves and washed for 1 min in 0.3% sodium hypochlorite solution and in Milli-Q water for 2 × 2 min to remove any potential foreign DNA contamination. These eggs were collected in 1.5 ml Eppendorf tube and after short centrifugation, the remains of the water (supernatant) removed with a pipette. The sample was immediately transferred to ice and prepared for DNA extraction. DNA was extracted using Bionano Prep Animal Tissue DNA Kit for HMW DNA isolation according to the manufacturer’s instructions. Briefly, eggs were smashed with a sterile pellet pestle on ice in 500 μl homogenization buffer; the sample was fixed with 500 μl cold ethanol and incubated 60 min on ice, after that time the sample was centrifuged at 1500g for 5 min at 4 °C and the supernatant was discarded. Next, after resuspension in a homogenization buffer pellet, this was cast in four agarose plugs as described in the original protocol. Agarose plugs were incubated with Proteinase K and Lysis buffer solution for 2 h with intermittent mixing. After that time, the digestion solution was replaced with a freshly made one and incubated overnight with intermittent mixing. According to the original protocol, after RNase A digestion and plug washing, DNA was recovered by incubation of the plugs in TE buffer, followed by plug melting and addition of agarase. Recovered DNA was dialysed and homogenized on a membrane for 45 min at room temperature and transferred to a clean tube with a wide bore tip.SequencingSequencing was done using Oxford Nanopore Technologies (ONT, MinION). Isolated DNA purified using AMPure XP beads and resuspended in H2O before library preparation. Two separate libraries were prepared using ligation sequencing kit, SQK-LSK109 and Rapid Sequencing Kit SQK-RAD004, respectively, according to the manufacturer’s protocols and were sequenced on a FLO-MIN106 R9.4.1 SpotON flow cell on a MinION Mk 1B sequencer (ONT). The total yields from sequencing were 484,700 reads (2,417,068,187 nt) with a read-N50 of 10,044 nt (ranging from 216,403 to 100). Base calling of the raw reads was done using Guppy (v.3.3) resulting in a total sum of the reads 7,979,616,172, equivalent to 26× coverage aiming for a genome of 300 megabases (Mb). The reads N50/N90 were estimated at 7,958/1,719.Assembling reference genomeReads aligning with the Mitochondrion genome were identified using BLASTN and filtered from the raw reads before assembling the genome. The remaining ONT reads were assembled using the Flye software (v.2.6), with –min-overlap 3,000 to increase stringency at the initial overlay step, and default parameters including five rounds of polishing through consensus, contigs were additionally polished two times with Medaka (v.0.11.2). Illumina paired-end RNA dataset is assembled using CLC Assembler (CLC Assembly Cell). Both RNA assemblies and paired-end 10X genomic dataset (unpublished data) were mapped onto the contigs using minimap2 (v.2.16) and BWA mapper (v.0.7.17), respectively, and the assembly was further polished using PILON (v.1.20) to error correct potential low-quality regions. The resulting assembly yielded a genome of 307 Mb, assembled into 1,533 contigs ranging from 10,840,357 to 100 basepairs (bp) and an assembly-N50 of 1.670 Mb. Moreover, the BUSCO completeness analysis using the Arachnida (odb10) reference set confirmed our assembly represents the complete genome C:94.8%(S:89.1%,D:5.7%),F:0.9%,M:4.3%,n:2934 (=arachnida_odb10), only missing 126 genes from the whole reference set. Knowing that BUSCO only gives a rough estimation, we remain confident that this assembly represents well the bulb mite genome.Flow-cytometryWhole individual R. robini were homogenized in 500 μl of ice-cold LB01 detergent buffer along with the head of a male Drosophila melanogaster (1 C = 0.18 pg) as an internal standard. The homogenized tissue was filtered through a 30-μm nylon filter. Then 12 μl of propidium iodide with 2 μl of RNase was added, and stained for 1 h on ice in the dark. All samples were run on an FC500 flow cytometer (Beckman-Coulter) using a 488-nm blue laser, providing output as single-parameter histograms showing relative fluorescence between the standard nuclei and the R. robini nuclei. Six replicate samples were run to account for variation in fluorescence outputs. The genome size of R. robini was estimated at 0.30 pg, or about 293 Mb, and consistent with estimation of size from the genome assembly described above.Mitochondrial genomeONT reads aligned with R. robini mitochondrion genome were de novo assembled with Flye (v.2.6) assembler and polished with Racon. Mitochondrion genome is assembled in one single contig with a size of 15,335 bases.Gene predictionOn the polished final genome, protein coding genes have been predicted. For this, AUGUSTUS was used including hints coming from R. robini RNA-sequencing (RNA-seq) (samples SRR3934324, SRR3934325, SRR3934326, SRR3934327, SRR3934328, SRR3934329, SRR3934330, SRR3934331, SRR3934332, SRR3934333, SRR3934335, SRR3934337, SRR3934338 and SRR3934339 from the PRJNA330592 BioProject deposited at the National Center for Biotechnology Information (NCBI) Short Read Archive) and proteins coming from highly curated Tetranychus urticae (v.2020-03-20) as well as proteins from the previous version of the unpublished, Illumina-sequenced R. robini genome (https://public-docs.crg.eu/rguigo/Data/fcamara/bulbmite.v4a/). The PE RNA-seq reads were mapped on the genome using HISAT2 (-k 1 —no-unal) and further processed with Regtools to extract junction hints and filtered for junctions with a minimum coverage of 10. All the RNA-seq reads were also assembled with CLC Assembly Cell (v.5.2.0) software, setting the word size for the Bruijn graph at 50 and maximum bubble size at 31. The reads were assembled into 689,563 contigs (ranging from 10,675 to 180 bp), which were later mapped on the genome with GenomeTheader to generate complementary DNA hints. Protein hints were generated by using with Exonerate (v.2.2) with Protein2Genome model. To reduce the amount of overprediction due to repeated elements (transposable elements, simple sequence repeats) we de novo predicted high abundant repeats using RepeatModeler. The accompanying parameter file for extrinsic data for AUGUSTUS was adapted to include these hints as well as the softmasking of the genomic sequence. The resulting gene predictions from AUGUSTUS were further curated with EvidenceModeler using the same extrinsic data. The BUSCO analysis confirmed that our gene prediction indeed captured the expected genes well (C:94.6%(S:86.3%,D:8.3%),F:0.4%,M:5.0%,n:2934 (=arachnida_odb10)). The final predicted gene set was subsequently processed to be uploaded into ORCAE (https://bioinformatics.psb.ugent.be/orcae/overview/Rhrob)98.ResequencingGenomic sampling and mappingFor genomic analyses we sampled material from each of the morph selection lines (n = 8) at F1, F12 and F29. Following the experimental evolution protocols, after the first 24 h of egg laying all adults were transferred to a new container (described above) for a second 24 h to lay eggs and from these second dishes genomic material was sampled. On maturation, adults were transferred to and kept for 3 days in containers. Adults were then randomly selected and placed into Digestion Solution for MagJET gDNA Kit (F1&12) or ATL buffer (F29) before freezing at −20 °C. From each population two samples were collected consisting of 100 individuals (1 × 100 females and 1 × 100 males of random morph), the two samples separated by sex were used as technical replicates. The tissue from the 100 individuals within each sample was homogenized and DNA was extracted by Proteinase K digestion (24 h) followed by standard procedures using MagJET Genomic DNA Kit (ThermoScientific, F1&12) or DNeasy Blood and Tissue (Qiagen, F29). DNA concentration was controlled with the Qubit double-stranded DNA HS Assay Kit and DNA quality was assessed on agarose gels. The library preparation was performed using NEBNext Ultra II FS DNA Library Prep kit for Illumina.Whole-genome resequencing was carried out by National Genomics Infrastructure (Uppsala, Sweden) using the Illumina Nova-Seq 6000 platform with S4 flow cell to produce 2 × 150 bp reads (average 160.7 × 106: range 130.7 × 106 − 189.9 × 106). Adaptors were trimmed from reads using Trimmomatic99 software (v.0.39) and unpaired reads discarded. Fastq files were mapped to the assembled genome with bwa mem100 (v.0.7.17-r1188) using default settings. Sam files were converted to bam files, sorted, duplicates marked and ambiguously mapped reads removed using samtools101 (v.1.9). On average, 90% (range, 86–93%) of the reads from each sample were mapped successfully, of which an average of 17% (range, 15–19%) were marked as duplicates. This left us with an average of 117.7 × 106 pair end reads per sample, ranging between 99.6 × 106 and 145.9 × 106 (Supplementary Table 1).Genomic analysisFile preparation and filteringPreparation of files used in genomic analysis was done as follows: bam files were converted to a pile-up file using samtools, following which indels and surrounding windows (5 bp either side) were filtered, using identify-genomic-indel-regions.pl and filter-pileup-by-gtf.pl in PoPoolation102 (v.1.2.2) to avoid false SNPs, with the resulting filtered pile-up file converted to a sync file using mpileup2sync.pl in PoPoolation2 (ref. 103) (v.1.201). Using custom python scripts, the distribution of coverage from each sample (single sex) was determined by recording the coverage of positions every 10 kb across the genome from the sync file to give information on expected coverage (Supplementary Fig. 1). On the basis of this, we filtered the sync and pile-up files to contain only regions within a range of informative coverage, where the mean coverage of all samples at every position was between 50% of the expected coverage and 200% of the expected coverage (56×, range 23−112×). The pile-up and sync files containing individual male and female samples (48 in total) were then merged by sex to give files containing allele frequencies from 24 samples (eight populations across three generations), each consisting of allele frequencies of 200 individuals (100 males and 100 females, above) and used in all subsequent analysis (unless stated otherwise). Similarly, we drew coverage of a position every 10 kb from each sample in the sex-merged sync file to determine a distribution from which we decided to subsample to (Supplementary Fig. 1). We putatively identified X-linked contigs (below) and excluded them autosomal analysis. A similar, but, separate analysis on genes and SNPs from X-linked contigs was performed by using different parameters (below).Estimating nucleotide diversityUsing PoPoolation we determined various estimates of genetic diversity per sample (that is, 24 sex-merged samples). The pile-up file from each sample was subsampled using subsample-pileup.pl to a coverage of 63× (max coverage, 252×) to standardize estimations of genetic diversity across the genome, between populations and across generations. First, nucleotide diversity (Tajima’s Pi, π) and number of segregating sites (Watterson’s theta, ϴ) were estimated within genes. We performed analysis of exons using Syn-nonsyn-at-position.pl, in which genetic diversity of synonymous and non-synonymous positions were determined. Further analysis of overall genetic diversity within exons and introns were performed using Variance-at-position.pl, Tajima’s D (D) also estimated in the former. We used a minimum count of three (equal to a minor allele frequency of roughly 5%) for a SNP to be called, and a phred score >30 and a pool size of 400. Further analysis using 10 kb sliding windows (step size 10 kb) across the genome were performed using Variance-sliding.pl, and also included estimation of D. Estimates of D require the minimum count to be 2, but otherwise all the same parameters were used.We filtered genes to be included in our analysis (and all subsequent analysis) on the basis of a number of criteria. On the basis of extensive RNA-seq data from both males and females (Plesnar-Bielak, unpublished data with NCBI accession number PRJNA796800), we only included genes in our analyses that were expressed at a mean level of fragments per kilobase of transcript per million mapped reads >1 across 72 samples originated from both sexes and both morphs rearing in three different temperatures (18, 23 and 28 °C). A further filtering step was performed to remove genes with inconsistent mapping between samples, only genes with >60% exons mapped to (calculated from positions used to calculate parameters in the Syn-nonsyn-at-position.pl π outputs), with 63−252× coverage, in all 24 samples were included in the analysis. The final dataset contained 13,389 autosomal genes and subsequently used to filter other datasets to retain this set of genes only (see Supplementary Table 8 for a list of genes). Similarly, windows were discarded from outputs if 60% of genes being mapped to in all 24 samples) and reducing the final X-linked dataset to contain fewer than 200 genes. We therefore opted to reduce the target coverage further to 40×, in an attempt to retain more genes. This slight reduction of target coverage increased the number of genes in the final dataset substantially to 587 genes. We therefore opted to use a minimum coverage of 40× in all analysis of X-linked SNPs, genes and windows.Diverging SNPsTo determine divergent SNPs between F- and S-lines, we extracted the allele frequencies of all samples from the sex combined sync file. Samples from F29 were then used to filter the entire dataset to only contain SNPs on the basis of a number of criteria. First, positions within all samples were required to have a coverage >63× and 5% (that is, the average of all samples but not necessarily above >5% in all samples). Thus, our dataset contained only positions with the target coverage in all F29 samples and in which polymorphisms were unlikely to be a consequence of sequencing errors. After this filtering we were left with roughly 6 million SNPs used in further analysis. We performed a GLM, at each position by comparing the count of the major allele against counts of minor alleles at F29, to determine consistent allele frequency changes between treatments70. If any population had minor or major allele count of 0, +1 was added to minor and major alleles from all samples. To correct for multiple testing, we converted P values to q values using the qvalues R package (v.2.14.1)104 and applied a FDR with a q 900,000), GLMs were performed (identical to above) on the simulated major and minor allele counts. Using a FDR with a q  More

It’s surprisingly hard to catch a sloth. Although they’re slow — very, very slow — if you climb a tree to catch one, it will move along to the next tree. Once you climb the new tree, it will move back again.My team does this regularly, as we conduct the Sloth Backpack Project, a data-logging initiative here in Costa Rica, where many sloths coexist with people. In 2017, I wanted to do more than research, so I started the Sloth Conservation Foundation.In this photograph, I’m fitting a backpack to a brown-throated three-fingered sloth (Bradypus variegatus) that we named Baguette, after a nearby bakery. The backpack will collect data on her location, movement and living patterns.We had found Baguette about 20 minutes earlier, balancing atop construction fencing as she attempted to escape two pit bulls. Baguette wasn’t all that grateful. She’s a feisty old girl. She’s old: she’s missing fingers, and she’s got scars on her face.I adore sloths, but I also envy them. They’re a powerful symbol of the slowness that our society needs more of. They don’t let anything stress them out unless it’s really important — they just get on with life.The backpack project will help us to understand sloth behaviour, so we can better protect them as the urban environment grows. This year, I received a €50,000 (US$52,220) Future For Nature award, which we will use to train a dog to detect sloth faeces. We can use faeces as a proxy for sloth numbers and locations in the region, and ultimately work out the boundaries of the species, how fast populations are declining and which conservation measures work.I’m happy I’ve moved away from academia — I can put all my energy into conservation as opposed to bashing out papers. That’s what I feel ecology should focus on — how we can use what we’re learning to give back to other species. More

This section is devoted to show results and to highlight eventual effects of the interplay between the nonlinearity characterizing the system dynamics and the presence of noisy fluctuations for the irradiance variable.Analysis of experimental dataThe need of taking into account noisy fluctuations of such an environmental variable is well demonstrated in Fig. 1. In the first panel (a) the experimental time behaviour of the irradiance is shown. This noisy curve is based on the experimental data (purple points) of the Boussole buoy located in the Gulf of Lion, collected over a period of nine years, precisely from 2004 to 2013. The time series of the experimental data presents quite a few gaps in time due to the malfunction of the buoy. This aspect has been remedied by merging the experimental data with those of the OASIM model validated for the Boussole site61 (yellow points). The latter is a multispectral atmospheric radiative transfer model that is in turn forced by experimental-model data based on ECMWF ERAINTERIM reanalyses which provide, for example, cloud cover data. The radiative model is partly stochastic since it considers the effects stemming from the presence of clouds, averaged along a single day (this explains why the yellow points are slightly less scattered). We see that the OASIM model accurately reproduces the profile which emerges from the experimental data. Further, we stress that the experimental data are only used in this initial analysis. In the biogeochemical simulations the irradiance signal is fully reconstructed starting from a realistic seasonal cycle combined with a range of different random fluctuations, and the information from OASIM is not used. In the second panel (b) the daily (black points) as well as the three-month (red points) running mean of the experimental series are plotted. Figure 1c shows the irradiance noisy fluctuations (INF) which have been obtained by subtracting the three-month running mean curve (3MRM, red curve in Fig. 1b) from the daily running mean one (DRM, black curve in Fig. 1b) and normalizing with respect to the mean of the 3MRM ((overline{3MRM})), namely (INF = (DRM – 3MRM) / overline{3MRM}). We see that a seasonal overall trend with higher oscillations during the winter time can be seen, implying that the characteristics of the noise may change over the year. Moreover, a slight imbalance between positive and negative values of the noisy fluctuations (that is, different values of the maximum fluctuation intensity) is present. The physical reason for the occurrence of such an aspect can be ascribed to the fact that the maximum value of solar irradiance corresponds to that measured during a sunny day. Conversely, the minimum level tends to zero corresponding to a dense darkness. While the former is close to the mean value of the solar irradiance (most of all in summer), the latter is much further away and then a natural asymmetry arises in the random fluctuations. However, it should be noted that, apart from the intense spikes, the asymmetry is not so pronounced, as proved by the mean value (red line in Fig. 1c) which is practically zero, namely (0.4%) of the (overline{3MRM}). Therefore, basing on this last observation, to model the noise affecting the irradiance dynamics, as a first approximation we consider a symmetric Gaussian autocorrelated noise as described in the next subsection.On the basis of such experimental results, we postulate the hypothesis that random fluctuations of light cannot be neglected, most of all in the study of ecological systems where light profoundly determines the system dynamics, governing fundamental processes at the basis of of the food web.Figure 1(a) Experimental data (purple points) of the stochastic solar irradiance collected by the Boussole buoy in a time-window of 9 years (2004-2013); the yellow points are the data generated by the OASIM model used to fill the gaps present in the experimental time-series due to malfunctioning of the buoy. (b) Daily (black points) and three-month (red points) running mean of the light curve in panel (a). (c) Irradiance noisy fluctuations (INF), obtained by subtracting the three-month running mean curve (3MRM) from the daily running mean one (DRM) and normalizing with respect to the mean value of 3MRM ((overline{3MRM})), namely (INF = (DRM – 3MRM) / overline{3MRM}); the red line represents the mean value of such fluctuations. Data already presented and validated in61.Full size imageSolar irradianceThe solar irradiance forcing is derived considering a deterministic seasonal oscillation combined with an Ornstein-Uhlenbeck process. The coefficient of variation (CV) of simulated light forcing, Fig. 2, (CV=sigma / mu) ((mu) and (sigma) being mean value and standard deviation calculated over both time and numerical realizations), is shown for 231 (D-tau) pairs. D and (tau) represent the intensity of a Gaussian noise source and the auto-correlation time of the fluctuations, respectively (see Eqs. (2) and (3)).Each pixel represents the mean value on time of CV calculated with respect to 1000 different stochastic realizations. Figure 2Coefficient of variation ((CV=sigma / mu)) of irradiance resulting from numerical integration of model equations for 231 (D-tau) different scenarios.Full size imageIt is easy to see the agreement between the results obtained from the numerical integration and the theoretical ones derivable from Eq. (5) by putting (text {var}{F_L(0)}=0) and (t gg 1), getting (sigma ^2_L=D / 2tau). In Fig. 2, indeed, the maximum values of (sigma) lie in the upper left part of the plot corresponding to small (high) values of (tau) (D). As it is clear the values of D have been chosen in order to obtain a relative standard deviation ranging from (5%mu) to (60%mu). We underline that, in this case, it is possible to interchangeably consider (sigma) and CV since the dependence of CV on D and (tau) does not differ from that of (sigma) (meaning that the dependence of (sigma) is not altered by dividing by (mu)) (results not shown).Effects on population dynamicsIn this section the noise-induced effects on the population dynamics are examined. The nine planktonic populations present a different qualitative behaviour of the CV, compared to that of the irradiance. In this case, the CV is characterized by a strong non-monotonic dependence on the parameter (tau). This aspect can be appreciated in Fig. 3 where different curves of CV versus the time correlation parameter are shown for different fixed values of D.Figure 3Coefficient of variation ((CV=sigma / mu)) of the nine planktonic populations resulting from numerical integration of model equations plotted versus the considered values of (tau); the different curves are related to different values of the noise intensity D.Full size imageThe existence of a maximum value for CV can be appreciated for each species. Although the qualitative behaviour is the same for all strains, particular attention has to be payed on diatoms and nanoflagellates. All the other species, indeed, present a percent variation of standard deviation between (2%) and (15%). In the case of nanoflagellates, instead, the D-dependent range is (20-90%), while diatoms reach values over the (100%) for the highest values of D. Therefore, these two species, in particular, and the whole system, in general, are extremely sensitive to the auto-correlation time which characterizes the noise.We note that the different curves related to the different selected values of D approach the horizontal axis, tending asymptotically to vanish as (tau) increases. Such a behaviour can be explained by the fact that high values of (tau) give rise to a more correlated dynamics, so that (tau rightarrow infty) implies fully correlated time-behaviours corresponding to the deterministic case. In this instance, then, all the different realizations give the same results, making the standard deviation vanish. The same happens, independently of the value of (tau), for low values of noise intensity for which the corresponding curves approach the same almost vanishing value (see orange, gray and yellow lines). Differently from the previous case, when (tau rightarrow 0) the noise tends to a delta-correlated noise, that is a white noise; for (tau ne 0), instead, the noise spectrum is not flat, being characterized by a Cauchy-Lorentz distribution. The strong nonmonotonicity of CV with respect to (tau), emerging when there are relatively high values of CV, implies a greater variability of the system biomass. Lower values of CV indicate that the system dynamics is less influenced by the presence of noise where very little or no differences with respect to the deterministic case are present. Conversely, high values of CV clearly demonstrate the remarkable signature of the presence of an impacting noise source. It is interesting to note that the noise influence on the ecosystem strongly depends on both (tau) and D, that is, just an intense noise is not enough to generate a greater response of the ecosystem. In particular, experimental data are characterized by a CV approximately equal to 0.361, which corresponds to values of D and (tau) lying on the diagonal strip in Fig. 2 ranging from ((tau ,D)=(0.5,10^4)) to ((tau ,D)=(365,10^7)). Finally we note the presence of a noise suppression effect. High values of D, indeed, can generate slight effects when the correlation time (tau) does not take on suitable values.The results shown here are an extension of the previous work by Benincà et al.56. There, the authors analyse a simpler, less realistic model of two interacting populations, whose dynamics is affected by a randomly fluctuating temperature. In that case, moreover, the deterministic oscillations of the temperature are suppressed, and the system exhibits intrinsic Lotka-Volterra oscillations whose frequency match with the characteristic one(s) of the noise. On the contrary, here, the observed maximum response (see Fig. 3) cannot be interpreted as a synchronization effect, since our model does not present intrinsic Lotka-Volterra-like oscillations and the periodic population variability is only due to the deterministic forcing(s).The nonmonotonic behaviour of the CV can be then interpreted as the signature of the intimate interplay between the ecological system and the noise. This interplay, indeed, has a pivotal role in both determining the dynamics of the populations and defining the characteristics of the ecosystem.In Fig. 3 it can be observed that the value of (tau) for which CV is maximum strongly depends on the noise intensity D. In particular, it is possible to note that the peaks in Fig. 3 move towards higher values of (tau) as the noise intensity increases. Thus, Fig. 3 demonstrates that the maximum-response effect to the random fluctuations is sensitive to the noise intensity D.However, it is important to underline that the response of the system to the noisy signal does not depend on the yearly oscillations induced by the deterministic forcings. Indeed, by considering constant the deterministic part of all external forcings (temperature, irradiance, wind and salinity), the non monotonic behaviour of CV with respect to both (tau) and D is still present, provided that the populations are not extinct (plot not shown). In this scenario indeed, besides dinoflagellates, diatoms and nanoflagellates are practically extinct as well, exhibiting thus a constant vanishing variance. All the other strains, instead, present qualitatively the same nonmonotonicity with only slight differences (shift of the peaks and different mean values of the CV curves), probably due to the extinction of diatoms and nanoflagellates which causes relevant differences in the system dynamics. More specifically, the system’s response seems to depend on both the noise intensity and the correlation time (see Fig. 3).In this scenario (absence of seasonal driving) we have studied the dependence on both parameters D and (tau) of the probability density functions (PDFs) of the non-vanishing populations. In Fig. 4, the PDFs of bacteria (B1), picophytoplankton (P3), microzooplankton (Z5) and etherotrophic nanoflagellates (Z6) are plotted for (tau =0.5) and eight different values of the parameter D.Figure 4Dependence of the probability density functions of non-vanishing populations on the parameter D for (tau =0.5). The curves are normalized within the interval taken into account. For this reason the relative peaks of the curves in the bottom panels have different values compared to those of the top panels. However, the figure aims at showing the existence of the value of the noise intensity for which the system is more sensitive as well as the generation of a stationary out-of-equilibrium state induced by the noise.Full size imageWe see that the mean value and the variance of these populations are strongly affected by the presence of random fluctuations in the irradiance. Specifically, as the noise intensity increases the mean values of picophytoplankton and bacteria concentrations exhibit a shift. In particular, the results indicate that picophytoplankton is disavantaged by the presence of a noisy component in the irradiance, which indeed tends to inhibit its ability to absorbe the solar light, slowing down its growth. As a consequence, since phytoplankton and bacteria compete for the same resources, as the former declines the latter are favoured, with a compensation mechanism which allows their predators (zooplankton populations) to be almost not affected by the noisy behaviour of the irradiance. Further, we note that for intermediate values of the noise intensity ((D = 10^4 – 10^5)) a maximum of the variance occurs (the PDFs are clearly spread on a wider range of values). Such an effect indicates that the noisy behaviour of irradiance strongly influences the whole ecosystem dynamics. Moreover, the nonmonotonic behaviour of the variance (its PDFs become larger and then tighter again as the noise intensity increases) indicates that the noise pushes the ecosystem away from equilibrium, driving it towards a non-equilibrium steady state. Finally, we note that the nonmonotonic behaviour of CV as a function of the noise intensity remains also in the presence of seasonal driving.Figure 5Coefficient of variation ((CV=sigma / mu)) of nine planktonic populations resulting from numerical integration of model equations plotted versus the considered values of D; different curves correspond to different values of the correlation time (tau).Full size imageFigure 5 shows indeed the nonmonotonic response of the ecosystem to the change of D when the deterministic seasonal cycling of the four environmental parameters (temperature, irradiance, wind and salinity) is present. It is easy to observe that also in this instance the major noise-induced effect appears in nanoflagellates and diatoms with a percent standard deviation of 50(%) and 100(%), respectively. The coalescence of different curves (related to different values of (tau)), as D decreases, is due to the fact that for (D rightarrow 0) the impact of the noise is negligible and the evolution of the system practically resembles the deterministic one. On the contrary, for higher values of D remarkable differences arise and clear peaks of CV appear in the considered range of variation.These plots show that, for a fixed value of (tau), there exists a value of the noise intensity for which the planktonic concentrations are maximally spread around their mean values (corresponding to the maximum value of CV and then of the variance). Moreover, such a nonmonotonic behaviour suggests the presence of a resonance, which can be interpreted as the effect of the interplay between the nonlinearity of the system and the environmental random fluctuations.Also in this case, the interplay between the two parameters D and (tau) in determining and characterizing the dynamics of the ecosystem transparently emerges. The value of D corresponding to the maximum value of CV, indeed, basically depends on the specific value of (tau).Finally, we point out that the different dynamic scenarios identified by the D-(tau) couples can be experienced by the system during the year, since the two parameters may seasonally vary depending on the different weather conditions. In other words, a seasonally varying noise (see Fig. 1c) may cause the nine populations explore different regions of the D-(tau) space during the year. Therefore, the results reported in this paper can highlight the detectable yearly variability of a marine ecosystem which does not stem from the deterministic seasonal variation of environmental parameters.Effects on the organic carbonIn this subsection the effects of the irradiance noise on the biogechemistry are analysed. In Fig. 6 the dependence on (tau) of both the CV [panel (a)] and the mean value concentration [panel (b)] of detritus, labile dissolved organic carbon (L-DOC), semi-labile dissolved organic carbon (SL-DOC) and gross primary production (GPP) are shown. All these biogeochemical properties are correlated with carbon cycling. Gross primary production is related to the amount of carbon entering in the ecosystem, and is related to the maximum energy available in the ecosystem progressively dissipated in the trophic web. Gross primary production is directly affected by light fluctuation and its CV shape is very similar to that of the irradiance, Fig. 2. We selected also detritus and DOC because they are important indicators for the carbon cycling dynamics and are related to the cycling of chemicals like heavy metals62. The different curves, related to different values of D, approach the same (vanishing) value for large (tau). As previously discussed for the CV [Fig. 6(a)] of biomass concentrations, this circumstance is due to the fact that, in this case, the system dynamics tends to the deterministic case, characterized by a unique possible realization implying a vanishing standard deviation. For high correlation times thus the system is insensitive to the noise intensity. On the contrary, for small values of (tau), different values of D lead to significant differences of the variance. In particular, detritus, L-DOC and SL-DOC exhibit a clear non-monotonic behaviour whose maximum value depends on the combined values of D-(tau). Only the GPP presents a decreasing monotonic behaviour.The dependence of the mean value concentration on (tau), instead, is qualitatively the same for all the four parameters. Also in this case we can note a diversification with respect to D occurring at small (tau) and a (deterministic) constant value arising for low (high) values of D ((tau)).These results manifest that not only the population dynamics, but also all the biogeochemical processes are profoundly affected by the presence of stochastic environmental variables. The values and the behaviour of the examined quantities are indeed determined by the intimate interplay between the intensity and the time correlation of the noise fluctuations.Figure 6(a) Coefficient of variation ((CV=sigma / mu)) and (b) mean value concentration ((mu)) of detritus, labile dissolved organic carbon (L-DOC), semi-labile dissolved organic carbon (SL-DOC) and gross primary production (GPP) resulting from numerical integration of model equations plotted versus the considered values of (tau); the different curves are related to different values of the correlation time D.Full size image More

Chapron, G. et al. Recovery of large carnivores in Europe’s modern human-dominated landscapes. Science 346, 1517–1519 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Lute, M. L., Carter, N. H., López-Bao, J. V. & Linnell, J. D. C. Conservation professionals’ views on governing for coexistence with large carnivores. Biol. Cons. 248, 108668 (2020).Article 

Google Scholar 
Gantchoff, M. G. & Belant, J. L. Regional connectivity for recolonizing American black bears (Ursus americanus) in southcentral USA. Biol. Cons. 214, 66–75 (2017).Article 

Google Scholar 
Ripple, W. J. et al. Status and ecological effects of the world’s largest carnivores. Science 343, 25 (2014).Article 
CAS 

Google Scholar 
Kays, R. et al. Does hunting or hiking affect wildlife communities in protected areas?. J. Appl. Ecol. 54, 242–252 (2017).Article 

Google Scholar 
Schipper, J. et al. The status of the world’s land and marine mammals: diversity, threat, and knowledge. Science 322, 225–230 (2008).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Smith, J. A., Wang, Y. & Wilmers, C. C. Top carnivores increase their kill rates on prey as a response to human-induced fear. Proc. R. Soc. B Biol. Sci. 282, 20142711 (2015).Article 

Google Scholar 
Stillfried, M., Belant, J. L., Svoboda, N. J., Beyer, D. E. & Kramer-Schadt, S. When top predators become prey: Black bears alter movement behaviour in response to hunting pressure. Behav. Proc. 120, 30–39 (2015).Article 

Google Scholar 
Støen, O.-G. et al. Physiological evidence for a human-induced landscape of fear in brown bears (Ursus arctos). Physiol. Behav. 152, 244–248 (2015).PubMed 
Article 
CAS 

Google Scholar 
Evans, M. J., Rittenhouse, T. A. G., Hawley, J. E. & Rego, P. W. Black bear recolonization patterns in a human-dominated landscape vary based on housing: New insights from spatially explicit density models. Landsc. Urban Plan. 162, 13–24 (2017).Article 

Google Scholar 
LaRue, M. A. et al. Cougars are recolonizing the midwest: Analysis of cougar confirmations during 1990–2008. J. Wildl. Manag. 76, 1364–1369 (2012).Article 

Google Scholar 
Cove, M. V., Fergus, C., Lacher, I., Akre, T. & McShea, W. J. Projecting mammal distributions in response to future alternative landscapes in a rapidly transitioning region. Remote Sens. 11, 2482 (2019).ADS 
Article 

Google Scholar 
Frid, A. & Dill, L. Human-caused disturbance stimuli as a form of predation risk. Conserv. Ecol. 6, 25 (2002).
Google Scholar 
Clinchy, M. et al. Fear of the human “super predator” far exceeds the fear of large carnivores in a model mesocarnivore. Behav. Ecol. 27, 1826–1832 (2016).
Google Scholar 
Suraci, J. P., Clinchy, M., Zanette, L. Y. & Wilmers, C. C. Fear of humans as apex predators has landscape-scale impacts from mountain lions to mice. Ecol. Lett. 22, 1578–1586 (2019).PubMed 
Article 

Google Scholar 
Gaynor, K. M., Hojnowski, C. E., Carter, N. H. & Brashares, J. S. The influence of human disturbance on wildlife nocturnality. Science 360, 1232–1235 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Smith, J. A., Thomas, A. C., Levi, T., Wang, Y. & Wilmers, C. C. Human activity reduces niche partitioning among three widespread mesocarnivores. Oikos 127, 890–901 (2018).Article 

Google Scholar 
Tucker, M. A. et al. Moving in the Anthropocene: Global reductions in terrestrial mammalian movements. Science 359, 466–469 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
Carter, N. H., Brown, D. G., Etter, D. R. & Visser, L. G. American black bear habitat selection in northern Lower Peninsula, Michigan, USA, using discrete-choice modeling. Ursus 21, 57–71 (2010).Article 

Google Scholar 
Naidoo, R. & Burton, A. C. Relative effects of recreational activities on a temperate terrestrial wildlife assemblage. Conserv. Sci. Pract. 2, e271 (2020).
Google Scholar 
Geffroy, B., Samia, D. S. M., Bessa, E. & Blumstein, D. T. How nature-based tourism might increase prey vulnerability to predators. Trends Ecol. Evol. 30, 755–765 (2015).PubMed 
Article 

Google Scholar 
Geffroy, B. et al. Evolutionary dynamics in the Anthropocene: Life history and intensity of human contact shape antipredator responses. PLoS Biol. 18, e3000818 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Beeco, J. A., Hallo, J. C. & Brownlee, M. T. J. GPS visitor tracking and recreation suitability mapping: tools for understanding and managing visitor use. Landsc. Urban Plan. 127, 136–145 (2014).Article 

Google Scholar 
Thorsen, N. H. et al. Smartphone app reveals that lynx avoid human recreationists on local scale, but not home range scale. Sci. Rep. 12, 1–13 (2022).Article 
CAS 

Google Scholar 
Evans, M. J., Hawley, J. E., Rego, P. W. & Rittenhouse, T. A. G. Hourly movement decisions indicate how a large carnivore inhabits developed landscapes. Oecologia 190, 11–23 (2019).ADS 
PubMed 
Article 

Google Scholar 
Carlos, A. W. D., Bright, A. D., Teel, T. L. & Vaske, J. J. Human-black bear conflict in urban areas: an integrated approach to management response. Hum. Dimens. Wildl. 14, 174–184 (2009).Article 

Google Scholar 
Johnson, H. E. et al. Human development and climate affect hibernation in a large carnivore with implications for human–carnivore conflicts. J. Appl. Ecol. 55, 663–672 (2018).Article 

Google Scholar 
Gould, N. P., Powell, R., Olfenbuttel, C. & DePerno, C. S. Growth and reproduction by young urban and rural black bears. J. Mammal. 102, 1165–1173 (2021).Article 

Google Scholar 
Ditmer, M. A., Noyce, K. V., Fieberg, J. R. & Garshelis, D. L. Delineating the ecological and geographic edge of an opportunist: The American black bear exploiting an agricultural landscape. Ecol. Model. 387, 205–219 (2018).Article 

Google Scholar 
McFadden-Hiller, J. E. Jr. & Belant, J. L. Spatial distribution of black bear incident reports in michigan. PLoS One 11, e0154474 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ladle, A., Steenweg, R., Shepherd, B. & Boyce, M. S. The role of human outdoor recreation in shaping patterns of grizzly bear-black bear co-occurrence. PLoS One 13, e0191730 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Wilbur, R. C., Lischka, S. A., Young, J. R. & Johnson, H. E. Experience, attitudes, and demographic factors influence the probability of reporting human–black bear interactions. Wildl. Soc. Bull. 42, 22–31 (2018).Article 

Google Scholar 
Lustig, E. J., Lyda, S. B., Leslie, D. M., Luttbeg, B. & Fairbanks, W. S. Resource selection by recolonizing American Black Bears. J. Wildl. Manage. 85, 531–542 (2021).Article 

Google Scholar 
Sun, C. C., Fuller, A. K., Hare, M. P. & Hurst, J. E. Evaluating population expansion of black bears using spatial capture-recapture. J. Wildl. Manage. 81, 814–823 (2017).Article 

Google Scholar 
Kautz, T. M. et al. Large carnivore response to human road use suggests a landscape of coexistence. Glob. Ecol. Conserv. 30, e01772 (2021).Article 

Google Scholar 
Michigan Department of Natural Resources (MIDNR) (2021).Blount, J. D., Chynoweth, M. W., Green, A. M. & Şekercioğlu, Ç. H. Review: COVID-19 highlights the importance of camera traps for wildlife conservation research and management. Biol. Cons. 256, 108984 (2021).Article 

Google Scholar 
Weather Atlas. https://www.weather-atlas.com/enEvans, J. S. Spatial Analysis and Modelling Utilities. Package ‘spatialEco’. https://cran.r-project.org/web/packages/spatialEco/spatialEco.pdf (2021).Díaz-Ruiz, F., Caro, J., Delibes-Mateos, M., Arroyo, B. & Ferreras, P. Drivers of red fox (Vulpes vulpes) daily activity: prey availability, human disturbance or habitat structure?. J. Zool. 298, 128–138 (2016).Article 

Google Scholar 
Moore, J. F. et al. Comparison of species richness and detection between line transects, ground camera traps, and arboreal camera traps. Anim. Conserv. 23, 561–572 (2020).Article 

Google Scholar 
Parsons, A. W. et al. Urbanization focuses carnivore activity in remaining natural habitats, increasing species interactions. J. Appl. Ecol. 56, 1894–1904 (2019).Article 

Google Scholar 
Allen, M. L., Sibarani, M. C., Utoyo, L. & Krofel, M. Terrestrial mammal community richness and temporal overlap between tigers and other carnivores in Bukit Barisan Selatan National Park, Sumatra. Anim. Biodiv. Conserv. 43(1), 97–107 (2020).Article 

Google Scholar 
Tian, C. et al. Temporal niche patterns of large mammals in Wanglang National Nature Reserve, China. Glob. Ecol. Conserv. 22, e01015 (2020).Article 

Google Scholar 
Meredith, M. & Ridout, M. Estimates of coefficient of overlapping for animal activity patterns. Package ‘overlap’. https://cran.r-project.org/web/packages/overlap/overlap.pdf (2020).RStudio Team. RStudio: Integrated Development for R. RStudio, PBC, Boston, MA. http://www.rstudio.com/ (2021).Ridout, M. S. & Linkie, M. Estimating overlap of daily activity patterns from camera trap data. JABES 14, 322–337 (2009).MathSciNet 
MATH 
Article 

Google Scholar 
Lashley, M. A. et al. Estimating wildlife activity curves: comparison of methods and sample size. Sci. Rep. 8, 1–11 (2018).CAS 
Article 

Google Scholar 
Rowcliffe, M. Animal Activity Statistics. Package ‘activity’. https://cran.r-project.org/web/packages/activity/activity.pdf (2021).MacKenzie, D. I., Nichols, J. D., Hines, J. E., Knutson, M. G. & Franklin, A. B. Estimating site occupancy, colonization, and local extinction when a species is detected imperfectly. Ecology 84, 2200–2207 (2003).Article 

Google Scholar 
Wei, T., & Simko, V. Visualization of a Correlation Matrix. Package ‘corrplot’. https://cran.r-project.org/web/packages/corrplot/corrplot.pdf (2017).Norton, D. C. et al. Female American black bears do not alter space use or movements to reduce infanticide risk. PLoS One 13, e0203651 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ditmer, M. A. et al. Behavioral and physiological responses of American black bears to landscape features within an agricultural region. Ecosphere 6, 1–21 (2015).Article 

Google Scholar 
Clark, D. et al. Using machine learning methods to predict the movement trajectories of the Louisiana black bear. SMU Data Sci. Rev. 5, 25 (2021).
Google Scholar  More

Growth conditions and plant materialsTwo experiments were conducted concurrently (sites A and B) in the same screen house in 2019 between the months of May and July at the Landmark University Greenhouse and Hydroponic Technology Center, a section of the Teaching and Research Farm of the University in Omu-Aran, Kwara State Nigeria. Experiment at site B was conducted simultaneously as A so as to validate the results of experiment A. Landmark University lies within Latitude 8° 7′ 26.21388″ and 5° 5′ 0.1788″. Both experiments (A & B) involved tomato (Solanum lycopersicum L. variety cherry) and cucumber (Cucumis sativus L. variety marketer) crops. For each crop, seeds were sown into a separate seed tray filled with coco peat (Coco peat, SRIMATHI EXPORT, INDIA). Cocopeat is the mesocarp tissue or husk after the grinding of coconut fruit. It has a lightweight and high water and nutrient holding capacities, it has an acceptable pH, electrical conductivity, and other chemical attributes27. Rice husk is the by-product of rice after milling. The rice husk used was collected from the rice processing mill of Landmark University. Rice husk is a highly porous and light weighted material with a very high specific area28.Two sets of seed trays (one for organic and another for inorganic fertilizers) were used each for tomato and cucumber crops in the nursery. Both were raised in the nursery for two weeks before transplanting. Black grow bags (30 × 17 cm) filled with a coco peat/rice husk (1:4 ratio by volume) mixture with a weight of about 10 kg were arranged in a screen house. Both the nursery and establishment of crop proper take place in a screen house. The screen house has a galvanized iron as the frame, a UV covering on top, side net for screening insect pests the floor fairly covered with granite. Temperature and relative humidity within the screen house during the period of the experiment was monitored using a Thermograph and a Barograph, and they were at an average of 31 °C and 75%, respectively.The grow bags were randomly placed in the screen house for the unbiased application of amendments. For both tomato and cucumber crops, the treatment comprised of six (6) levels of liquid organic fertilizer (5, 15, 25, 35, 45, 55 mL), in-organic fertilizer, and a control (ordinary borehole water). Levels of organic fertilizers were selected based on the recommendation of 20 mL of liquid organic fertilizer by29. The eight (8) treatments both for tomato and cucumber were arranged in a Completely Randomized Design replicated three times. One healthy plant was maintained per grow bag and four grow bags represent a treatment and there were 32 plants per block each for tomato and cucumber. For both crops, the experiment lasted for 90 days.Organic and in-organic nutrient solutionsThe liquid organic fertilizer used was obtained from the biomass of Mexican sunflower (Tithonia diversifolia). Fresh biomass (mainly leaves and stems) of the plant was collected from the Teaching and Research Farms of Landmark University, Nigeria. After rinsing, they were cut with a sterile knife into pieces of ≤ 1 cm size. A sample was taken for initial physicochemical analyses by grinding in a sterile mortal, diluted with sterile water and analyzed. The biomass was then soaked in sterile water inside a clean container, and allowed to ferment spontaneously for a period of 14 days. During the fermentation, samples were taken every 4 days for microbial analyses of the major players during the fermentation. At the end of fermentation, the mixture was separated using a sieve of mesh size ≤ 2 mm. The liquid portion was then refrigerated prior to the planting regime while another sample was taken to ascertain the physicochemical and microbial qualities of the produced liquid fertilizer. The chemical analysis is presented in Table 4. For inorganic fertilizer, Water soluble fertilizers employed in hydroponics were used (Hydroponics fertilizer, Anmol chemicals, India); calcium nitrate 650 mg L−1, potassium nitrate 450 mg L−1, magnesium 400 mg L−1, chelate 20 mg L−1, mono-ammonium phosphate 400 mg L−1. The electrical conductivity (EC) of the solution was 1.9 dS m-1.Irrigation and fertigationThe tomato and cucumber plants were fertigated morning and evening daily for one hour on each occasion according to the treatments. Preparation of the nutrient solution was with borehole water and was supplied to plants by an online pressure drip irrigation system set at 2.0 L h-1 using an arrowhead on each tomato and cucumber plant. Different tanks (250 L) were installed according to the various treatments making a total of 8 tanks. The organic fertilizer was diluted according to the various treatments equivalent to 1.25, 3.75, 6.25, 8.75, 11.25, and 13.75 L per 250 L of water respectively for 5, 15, 25, 35, 45, and 55 mL treatments. The nutrient solutions were refilled when the consumption is less than 20% of the initial volume (250 L) in the tank. One day per week, crops were irrigated with ordinary water to wash out pipes and prevent deposits of salts. The same concentration of nutrient was used from transplanting to the termination of the study for both tomato and cucumber crops, however, at the flowering of the crops, the volume of fertigation was increased to 3.0 L h-1 to be able to cope with the size of the plants.Trellising, pest and diseases controlFor both tomato and cucumber crops, plant vines were supported by twisting them around a wire that is- attached to the roof of the screen house and 2 m from the ground. Lateral outgrowths were cut off every week to ensure a sturdy single stem. Pests and diseases were scouted every day. Whiteflies, aphids, and other insects were controlled with orizon (Producer, location of producer) (active ingredient, acetamiprid, and abamectin) using 0.133% v/v. Fungi were controlled using ridomil gold (Producer, Location of producer) at 2% w/v.Determination of growth and yield of tomato and cucumberThree tomato and cucumber plants were randomly selected for each treatment for the determination of growth parameters (plant height, leaf area, number of leaves per plant, and stem diameter) at mid the flowering stage of tomato and cucumber plants.The leaf area of tomato was calculated using the model (A = KL2) developed by Lyon30, where L = Length of tomato leaf, K = constant which is 0.1551, and A = leaf area of tomato. Similarly, the leaf area of cucumber was calculated using A = 0.88LW – 4.27, where L = cucumber leaf length and W = cucumber leaf width, A = leaf area of cucumber31.Tomato fruits were ready for harvest from 65 days after transplanting, harvestings were done twice every week (Mondays and Fridays) for up to 85 days after transplanting. Similarly, harvesting of cucumber fruits started 35 days after transplanting and harvestings were also done twice a week (Mondays and Fridays), harvesting was carried out till 60 days after transplanting. Tomato and cucumber fruit yields were counted and weighed at each harvest.Analysis of tomato and cucumber leaves and fruitsAt the 50% flowering stage of tomato and cucumber plants, ten leaf samples were collected from each treatment. The leaf samples were oven-dried at 75 °C for 24 h and thereafter grounded. The grounded samples were later analyzed for nitrogen (N), phosphorous (P), potassium (K), calcium (Ca), and magnesium (Mg) content using the method of described by32. At harvest, four matured tomato and cucumber fruits of uniform size were selected per treatment, and their nutrient compositions were determined using the method of33.Statistical analysisAll data collected on the growth, yield, leaf, and fruit nutrient contents of tomato and cucumber were subjected to analysis of variance (ANOVA). The SPSS V 21.0 (New York, USA) software was used to perform ANOVA and Duncan’s multiple range test (DMRT) was used to compare means at a 5% probability level.
Ethical approvalI confirm that all the research meets ethical guidelines and adheres to the legal requirements of the study country.Compliance with international, national and/or institutional guidelinesExperimental research (either cultivated or wild), comply with relevant institutional, national, and international guidelines and legislation. Experimental studies were carried out in accordance with relevant institutional, national or international guidelines or regulation. More

Armbrust, E. V. The life of diatoms in the world’s oceans. Nature 459, 185–192 (2009).CAS 
Article 

Google Scholar 
Mann, D. G. The species concept in diatoms. Phycologia 38, 437–495 (1999).Article 

Google Scholar 
Smetacek, V. Diatoms and the ocean carbon cycle. Protist 150, 25–32 (1999).CAS 
Article 

Google Scholar 
Rynearson, T. A. et al. Major contribution of diatom resting spores to vertical flux in the sub-polar North Atlantic. Deep. Res. Part I Oceanogr. Res. Pap. 82, 60–71 (2013).CAS 
Article 

Google Scholar 
Allen, J. T. et al. Diatom carbon export enhanced by silicate upwelling in the northeast Atlantic. Nature 437, 728–732 (2005).CAS 
Article 

Google Scholar 
Boyd, P. W., Strzepek, R., Fu, F. & Hutchins, D. A. Environmental control of open-ocean phytoplankton groups: Now and in the future. Limnol. Oceanogr. 55, 1353–1376 (2010).CAS 
Article 

Google Scholar 
Hátún, H., Somavilla, R., Rey, F., Johnson, C. & Mathis, M. The subpolar gyre regulates silicate concentrations in the North Atlantic. Sci. Rep. 1–9 https://doi.org/10.1038/s41598-017-14837-4 (2017).Bopp, L. Response of diatoms distribution to global warming and potential implications: a global model study. Geophys. Res. Lett. 32, 2–5 (2005).Article 
CAS 

Google Scholar 
Warner, A. J. & Hays, G. C. Sampling by the Continuous Plankton Recorder survey. Prog. Oceanogr. 6611, 237–256 (1994).Article 

Google Scholar 
Edwards, M., Beaugrand, G., Reid, P. C., Rowden, A. A. & Jones, M. B. Ocean climate anomalies and the ecology of the North Sea. Mar. Ecol. Prog. Ser. 239, 1–10 (2002).Article 

Google Scholar 
Allen, S. et al. Interannual stability of phytoplankton community composition in the North-East atlantic. Mar. Ecol. Prog. Ser. 655, 43–57 (2020).Article 

Google Scholar 
Edwards, M. & Richardson, A. J. Impact of climate change on marine pelagic phenology and trophic mismatch. Nature 430, 881–884 (2004).CAS 
Article 

Google Scholar 
Wihsgott, J. U. et al. Observations of vertical mixing in autumn and its effect on the autumn phytoplankton bloom. Prog. Oceanogr. 177, 1157–1165 (2019).Article 

Google Scholar 
Kamykowski, D. & Zentara, S. J. Predicting plant nutrient concentrations from temperature and sigma-T in the upper kilometer of the world ocean. Deep. Res. Part A-Oceanogr. Res. Pap. 33, 89–105 (1986).CAS 
Article 

Google Scholar 
Kamykowski, D. & Zentara, S. J. Changes in world ocean nitrate availability through the 20th century. Deep. Res. Part I-Oceanogr. Res. Pap. 52, 1719–1744 (2005).Article 

Google Scholar 
Behrenfeld, M. J. et al. Climate-driven trends in contemporary ocean productivity. Nature 444, 752–755 (2006).CAS 
Article 

Google Scholar 
López-Urrutia, A., San Martin, E., Harris, R. P. & Irigoien, X. Scaling the metabolic balance of the oceans. Proc. Natl. Acad. Sci. USA 103, 8739–8744 (2006).Article 
CAS 

Google Scholar 
Richardson, A. J. & Schoeman, D. S. Climate impact on plankton ecosystems in the Northeast Atlantic. Science 305, 1609–1612 (2004).CAS 
Article 

Google Scholar 
Edwards, M., Johns, D., Leterme, S., Svendsen, E. & Richardson, A. Regional climate change and harmful algal blooms in the northeast Atlantic. Limnol. Oceanogr. 51, 820–829 (2006).Article 

Google Scholar 
Hinder, S. L. et al. Changes in marine dinoflagellate and diatom abundance under climate change. Nat. Clim. Chang. 2, 271–275 (2012).Article 

Google Scholar 
Batten, S. et al. CPR sampling: the technical background, materials and methods, consistency and comparability. Prog. Oceanogr. 58, 193–215 (2003).Article 

Google Scholar 
Reid, P. C. et al. The Continuous Plankton Recorder: concepts and history, from plankton indicator to undulating recorders. Prog. Oceanogr. 58, 117–173 (2003).Article 

Google Scholar 
Edwards, M., Beaugrand, G., Hays, G. C., Koslow, J. A. & Richardson, A. J. Multi-decadal oceanic ecological datasets and their application in marine policy and management. Trends Ecol. Evol. 25, 602–610 (2010).Article 

Google Scholar 
Edwards, M. et al. North Atlantic warming over six decades drives decreases in krill abundance with no associated range shift. Commun. Biol. 4, 1–10 (2021).Article 

Google Scholar 
Richardson, A. J. et al. Using continuous plankton recorder data. Prog. Oceanogr. 68, 27–74 (2006).Article 

Google Scholar 
Hélaouët, P., Beaugrand, G. & Reygondeau, G. Reliability of spatial and temporal patterns of C. finmarchicus inferred from the CPR survey. J. Mar. Syst. 153, 18–24 (2016).Article 

Google Scholar 
Owens, N. J. P. et al. All plankton sampling systems underestimate abundance: response to “Continuous plankton recorder underestimates zooplankton abundance” by J.W. Dippner and M. Krause. J. Mar. Syst. 128, 240–242 (2013).Article 

Google Scholar 
Jonas, T. D., Walne, A., Beaugrand, G., Gregory, L. & Hays, G. C. The volume of water filtered by a Continuous Plankton Recorder sample: the effect of ship speed. J. Plankton Res. 26, 1499–1506 (2004).Article 

Google Scholar 
O’Reilly, C. H., Zanna, L. & Woollings, T. Assessing external and internal sources of Atlantic multidecadal variability using models, proxy data, and early instrumental indices. J. Clim 32, 7727–7745 (2019).Article 

Google Scholar 
Qin, M., Dai, A. & Hua, W. Quantifying contributions of internal variability and external forcing to atlantic multidecadal variability since 1870. Geophys. Res. Lett. 47, 1–11 (2020).
Google Scholar 
Mann, M. E., Steinman, B. A. & Miller, S. K. Absence of internal multidecadal and interdecadal oscillations in climate model simulations. Nat. Commun. 11, 1–9 (2020).Article 
CAS 

Google Scholar 
Enfield, D. B., Mestas-Nuñez, A. M. & Trimble, P. J. The Atlantic Multidecadal Oscillation and its relation to rainfall and river flows in the continental U.S. Geophys. Res. Lett. 28, 2077–2080 (2001).Article 

Google Scholar 
Gray, S. T., Graumlich, L. J., Betancourt, J. L. & Pederson, G. T. A tree-ring based reconstruction of the Atlantic Multidecadal Oscillation since 1567 A.D. Geophys. Res. Lett. 31, 2–5 (2004).Article 

Google Scholar 
Edwards, M., Beaugrand, G., Helaouët, P., Alheit, J. & Coombs, S. Marine ecosystem response to the Atlantic Multidecadal Oscillation. PLoS One 8, e57212 (2013).CAS 
Article 

Google Scholar 
Jones, P. D., New, M., Parker, D. E. & Martin, S. & Rigor, I. G. Surface air temperature and its changes over the past 150 years. Rev. Geophys. 37, 173–199 (1999).Article 

Google Scholar 
Beaugrand, G. et al. Reorganization of North Atlantic marine copepod biodiversity and climate. Science296, 1692–1694 (2002).CAS 
Article 

Google Scholar 
Alvain, S., Moulin, C., Dandonneau, Y. & Bréon, F. M. Remote sensing of phytoplankton groups in case 1 waters from global SeaWiFS imagery. Deep Sea Res. Part I Oceanogr. Res. Pap. 52, 1989–2004 (2005).Article 

Google Scholar 
Kalnay, E. et al. The NCEP/NCAR 40-year reanalysis project. Bull. Am. Meteorol. Soc. 77, 437–471 (1996).Article 

Google Scholar 
Huang, B. et al. Extended reconstructed sea surface temperature, version 5 (ERSSTv5): upgrades, validations, and intercomparisons. J. Clim. 30, 8179–8205 (2017).Article 

Google Scholar 
Lam, N. S. N. Spatial interpolation methods: a review. Am. Cartogr. 10, 129–150 (1983).Article 

Google Scholar 
Beaugrand, G., McQuatters-Gollop, A., Edwards, M. & Goberville, E. Long-term responses of North Atlantic calcifying plankton to climate change. Nat. Clim. Chang. 3, 263–267 (2012).Article 
CAS 

Google Scholar 
Beaugrand, G. et al. Prediction of unprecedented biological shifts in the global ocean. Nat. Clim. Chang. 9, 237–243 (2019).Article 

Google Scholar  More

Overall characteristics of air pollutantsThe results of previous studies indicated that local pollution is highly important in determining the emissions of air pollutant. Therefore, in this study, we estimated the changes in pollution and the AQI between the pre-COVID and COVID lockdown periods and among the different regions in Ji’nan. A comparison of the different pollutant concentrations analysed in this study shows that the concentrations of almost all pollutants decreased during the COVID lockdown period; only the concentration of O3 increased continuously as the COVID lockdown period progressed (Fig. 1).Figure 1Spatial distributions of the different observation sites and industrial enterprises above a designated size threshold in Ji’nan city. JCE, machine tool factory No. 2; LSX, technical college; JNS, Ji’nan fourth building group; KFQ, economic development zone; KGS, Kegansuo; LWZ, Laiwu memorial hall; NKS, Agricultural Scientific Institute; SZZ, Seed warehouse of Shandong Province; SJC, Ji’nan monitoring station; TXG, Taixing company; CQD, Changqing school. Red circles, red triangles and red squares represent stations in urban, urban-industrial and suburban regions, respectively. The map of Observation site was completed by the geostatistical analysis module of ArcGIS (version 10.3, https://developers.arcgis.com/).Full size imageDuring the observation period, the daily average mass concentrations of PM10, PM2.5, SO2, NO2, CO, and O3 in Ji’nan were 137.09 µg/m3, 101.35 µg/m3, 22.70 µg/m3, 39.77 µg/m3, 1.28 mg/m3, and 71.84 µg/m3, respectively (Fig. 2). The mass concentrations of PM10 and PM2.5 exceeded the daily average Grade I values (50 µg/m3 and 35 µg/m3) of the Ambient Air Quality Standard of China (CAAQS, GB 3095-2012) during the whole observation period. In contrast, the mass concentrations of NO2, SO2, CO and O3 were substantially lower than the daily average Grade I values (80 µg/m3, 50 µg/m3, 4 mg/m3 and 100 µg/m3, respectively) of the CAAQS each day. During the pre-COVID period, the daily average mass concentrations of PM10, PM2.5, SO2, NO2, CO, and O3 in Ji’nan were 177.03 µg/m3, 125.94 µg/m3, 26.39 µg/m3, 54.52 µg/m3, 1.59 mg/m3, and 60.72 µg/m3, respectively. The mass concentrations of all these pollutants, except NO2, CO and O3, exceeded the daily average Grade I values of the CAAQS. The mass concentration trends during the COVID lockdown period were consistent with those during the pre-COVID period, but there were significant differences in the concentrations between the periods. In summary, the air quality in Ji’nan was generally good from January 24 to February 7, 2020, mainly due to the strict prevention and control measures for COVID-19.Figure 2Temporal variations in the mass concentrations of air pollutants (PM10, PM2.5, NO2, SO2, CO and O3) at the urban site in Ji’nan during the observation period.Full size imageEffects of regional differences and lockdown on air pollutantsOur results reveal that the PM10, PM2.5, NO2, SO2, CO and O3 concentrations in the urban, suburban and urban-industrial regions differed significantly between the COVID lockdown and pre-COVID periods (Figs. 3, 4).Figure 3Mean concentrations (± SD, mg/m3) of PM10, PM2.5, NO2, SO2, CO and O3 during the pre-COVID and COVID lockdown periods in 2020; the values were determined by combining the urban, suburban and urban-industrial areas at the regional scale. *, ** and *** represent significant differences between the pre-COVID and COVID lockdown periods in the same region (Duncan test, *p = 0.05; **p = 0.01; ***p = 0.001), with nonsignificant results being excluded.Full size imageFigure 4General reductions in the concentrations of major air pollutants.Full size imageNOx, one of the most important pollutants and a major health hazard, was studied in different countries across the world during COVID-19-related lockdowns. In all three regions studied herein, the highest rate of reduction in NO2 concentrations was observed during the COVID lockdown period (Fig. 4), with the NO2 levels in the COVID lockdown period being 54.02% on average lower than those during the pre-COVID period (53.07% in urban area, 48.31% in the suburban areas and 55.74% in the urban-industrial area) (Fig. 4); this reduction is greater than that reported at other sites by 26–42%11 and 14–38%18 but lower than that (50–62%) in Barcelona and Madrid in Spain33. As shown in Fig. 3E, the NO2 concentrations in the urban, suburban and urban-industrial areas were significantly higher in the pre-COVID period than in the COVID lockdown period, with the pre-COVID the NO2 levels in the urban area being 13.46% and 27.63% higher than those in the suburban and urban-industrial areas, respectively. During the COVID-19 lockdown period, the NO2 levels in urban areas were 4.69% and 31.75% higher than those in the suburban and urban-industrial areas, respectively. Blocking and controlling the air pollution associated with COVID-19 has helped reduce ground NO2 levels34 and this effect might be correlated with the tropospheric NO2 column density27. Among all sources of NO2, automobile emissions and power generation are the most important5. A systematic review confirmed that a short-term increase in the NO2 concentration in urban areas correlates to an increase in the number of pneumonia hospitalizations5,35.The trends in the CO concentration were similar to those in the NO2 level. During the COVID-lockdown period, the average CO mass concentrations in the urban, suburban and urban industrial areas were 1.08 mg/m3, 1.16 mg/m3 and 1.14 mg/m3, respectively, which decreased by 27.78%, 29.46% and 36.61%, respectively, compared with those during the pre-COVID period. The highest levels of PM10 were also observed during the pre-COVID period in the urban, suburban, and urban-industrial areas in Ji’nan (Fig. 4). The reductions in PM2.5 and CO emissions in urban and urban-industrial areas are generally higher than those in suburban areas25, supporting our findings. Notably, PM2.5 and CO are generated mainly by construction activities and from road dust, natural soil dust and dust from urban-industrial activities36. In contrast, the differences in the PM10 concentrations among the three regions were not significant during either the pre-COVID period or the COVID-lockdown period (Fig. 3A), which suggests that particles in Ji’nan are strongly diffused. However, the COVID lockdown period had a significant effect on the PM10 concentrations, with 42.86%, 44.26% and 50.60% differences in the PM10 concentration between the pre-COVID and COVID lockdown periods in the urban, suburban and urban-industrial areas, respectively (average of 44.92%, Fig. 4). The main reasons for the decreases in the concentration of PM were the severe restrictions on vehicle traffic, the cessation of industrial activities, and the stopping of construction projects, which are important sources of floating dust in the urban air37. Despite the overall consistency among the observed changes in all regions for the different air pollutants (except O3), at the regional level, some differences were statistically significant, while others were not due to the variability among stations, with the differences being more pronounced at the urban, suburban and urban-industrial stations.O3 is a secondary pollutant involved in different atmospheric reaction mechanisms and acts as both a source and sink. Generally, the impact of lockdowns on O3 was mixed, with its levels generally falling within ± 20%38, but total O3 levels remained relatively stable18. In this study, by comparing the regional mean concentrations throughout the COVID-19 period, we found that O3 concentrations were higher during the COVID lockdown period than during the pre-COVID period, especially in the urban regions (Fig. 3). Furthermore, the mean O3 concentration at all stations during the COVID lockdown period was 37.42% higher than that during the pre-COVID period (46.84% in the urban areas, 18.27% in the suburban area, and 19.84% in the urban-industrial areas) (Fig. 4); this finding is consistent with the outcomes of other studies, which reported that O3 concentrations increased by (on average) 20% during lockdowns39, potentially due, in part, to atmospheric reactivity37. The higher lockdown O3 concentrations can be attributed to the following three reasons: (1) low PM concentrations can result in more sunlight passing through the atmosphere, encouraging increased photochemical activities and thus higher O3 production40; (2) a reduction in NOx emissions increases O3 formation41; and (3) lower PM2.5 concentrations means their role as a sink for hydroperoxy radicals (HO2) is less effective, which would increase peroxy radical-mediated O3 production42. During the pre-COVID period, the O3 levels were not significantly different among the region, and the same results were observed during the COVID lockdown period. However, in the urban and urban-industrial areas, the O3 levels during the COVID lockdown period were significantly higher than those in the pre-COVID period (p  More