Ecology

Seasonal changes in Photosystem II- and Photosystem I-related functional activities
To monitor photosynthetic performance, we recorded several PSII and PSI parameters during three consecutive growth seasons (2015–2016, 2016–2017, and 2017–2018) along with concomitant changes in daily air temperature and solar radiation (Fig. 1a–c). For simplicity, in the main figures we only present (here onwards) data from 2017 to 2018, which we divided into five distinct seasons, based on weather parameters: Summer (S, June–Aug), autumn (A, Sept–mid Nov), winter (W, mid Nov–mid Feb), early spring (ES, mid Feb–mid Apr), and late spring (LS, mid Apr–June). Data from the other two growth seasons are provided in Supplementary Figs. 1–3.
Fig. 1: Seasonal dynamics of weather and photochemical performance of PSII measured by chlorophyll fluorescence in Scots pine needles.

Changes in temperature (°C) (Left Y-axis) and solar radiation (watt m−2) (Right Y-axis) during 2015–2016 (a), 2016–2017 (b), 2017–2018 (c) measuring seasons. Seasonal dynamics of PSII photochemistry: (d) Changes in maximal (Fm) and basic (Fo) fluorescence. (e) Maximal quantum efficiency of PSII measured as Fv/Fm (f) Effective quantum yield of PSII (Φ(II)). Energy dissipation measured as regulated and non-regulated non-photochemical quenching: (g) Changes in NPQ with increasing PAR. (h) Induction of NPQ with constant actinic light. (i) Quantum yield of non-regulated non-photochemical quenching. Recovery of pine needles under artificial conditions at 80 µmol photons m−2 s−1 of light for 48 h with 18/6 photoperiod: (j) Changes in Fo. Fm and Fv/Fm. (k) Changes in Φ(II). (l) Induction of NPQ with constant actinic light. Quantum yields were calculated at actinic light illumination of 300 μmol m−2 s−1 in the light response curve and NPQ induction was measured either at 300/800 μmol m−2 s−1 of actinic light illumination. All measurements were taken after 30 min of dark adaptation at 4 °C from winter to late spring and room temperature in summer. Error bars indicate the mean ± SD (n = 3) for (d–g), and (i); ± SD (n = 4) for (j, k). Error bars indicate the mean ± SE (n = 3) for (h, l). The statistically significant mean differences (t-test) are marked by asterisks indicating the following confidence intervals: *≤95%; **≤99%.; ***≤99.9%. Exact p values are provided in the source data file.

Full size image

We found characteristic seasonal patterns in maximum PSII fluorescence (Fm) and maximum quantum efficiency of PSII (expressed as Fv/Fm), in accordance with earlier reports8,10. Fv/Fm was highest in S, fell with reductions in ambient temperatures during A and W (Fig. 1d, e), and was lowest in ES (63% lower than in S), when low temperatures coincided with rises in solar irradiance (Fig. 1a–c; see also Supplementary Fig. 2a). During LS, Fv/Fm gradually increased with increasing temperatures and peaked in S. The seasonal changes recorded in Fv/Fm were found to be mostly affected by changes in Fm rather than Fo (basic intrinsic fluorescence) (Fig. 1d). For a deeper understanding of PSII performance, the fraction of absorbed light energy utilized by PSII photochemistry15 (ΦPSII) was measured (Fig. 1f). During W ΦPSII decreased significantly and reached minimal values (19% of S values) in ES (Supplementary Fig. 2b).
Furthermore, we tried to quantify the amount of absorbed light energy thermally dissipated by non-photochemical quenching (NPQ)3. The component of NPQ that play a crucial role under fluctuating light conditions is the fast component (ΔpH-PsbS dependent-qE3 or zeaxanthin dependent-qZ4), which increased with increasing light intensities in all the seasons. In S and A, the fast component did not reach a stationary phase even at 1500 µmol m−2 s−1 illumination, while in W and LS samples a stationary phase was reached at 500 µmol m−2 s−1 and in ES already at 300 µmol m−2 s−1 (Fig. 1g). Most importantly, in ES inducible steady-state NPQ was much smaller and overall slower due to the smaller amplitude of the fast component (~50% less than in S) in ES (Fig. 1h). This is due to the fact that ES needles have already developed static NPQ. Instead, the quantum yield of non-regulated and/or constitutive loss (ΦNO) of energy was high during ES (Fig. 1i) (Supplementary Fig. 2d). This strongly suggests that this static NPQ is the fraction of absorbed light energy neither going to drive photochemistry (ΦPSII) nor thermally dissipated by rapid regulated NPQ processes (qE/qZ).
To confirm that this static quenching is a manifestation of a “sustained mode of quenching” we artificially relaxed the ES needles (hereafter, ER) for 48 h in low light (80 µmol m−2 s−1) with 18/6 h photoperiod. During this recovery, we observed no significant changes in Fv/Fm for 6–8 h (less than a 10% increase), a modest increase (30%) after 24 h (Fig. 1j), and almost complete (95%) recovery to S levels after 48 h. Fo did not change significantly throughout the recovery period (Fig. 1j), but changes in Fm followed the same recovery dynamics as Fv/Fm (Fig. 1j), although maximal fluorescence (Fm) was much higher in S than in 48 h ER samples. ΦPSII did not change significantly until 12 h of recovery but recovered to 50% and 90% of S levels after 24 and 48 h, respectively (Fig. 1k). qE/qZ-dependent fast NPQ kinetics followed similar patterns to ΦPSII (Fig. 1l). It should be noted that light-dependent induction of the fast component of NPQ (qE/qZ) was slower in the ER samples than in S samples, even after 48 h of recovery, which appears to be mainly due to photoinhibition.
Seasonal changes in the partitioning of absorbed light energy within PSI showed that Y(I) was highest in S, declined in A and W, and reached minimal values in ES (Fig. 2a). Although it decreased, Y(I) was much less affected than the PSII activity (ΦPSII) during the cold periods (Fig. 1f) (Supplementary Fig. 3a). Y(NA) followed the same seasonal pattern, with minimum values registered in ES and a steady recovery until S (Fig. 2b) (Supplementary Fig. 3b). In contrast to Y(NA), Y(ND) considerably increased during the cold periods and peaked (at 69% higher than S values) in ES, then gradually declined to minimum values in S (Fig. 2c, Supplementary Fig. 3c). This suggests that during early spring most of the absorbed light energy is dissipated by oxidized P700 (P700+) to prevent photoinhibition caused by over reduction of the acceptor side as iron-sulfur clusters get damaged16. During the first 3 h of the recovery period, ΦPSI did not significantly change, but it gradually returned to S levels over 48 h (Fig. 2d). Y(NA) did not change much during recovery (Fig. 2e). Y(ND) increased slightly during the first 3 h, but then gradually decreased, and ER samples did not quite reach S values. Based on these results, we conclude that the restoration of PSI light use efficiency takes longer than 48 h (Fig. 2f).
Fig. 2: Seasonal changes of PSI photochemistry in Scots pine needles.

Energy partitioning in PSI considering Y(I) + Y(ND) + Y(NA) = 1, where Y(I) (a), Y(NA) (b), Y(ND) (c) are the photochemical quantum yield of PSI (when P700 is reduced and A is oxidized), energy dissipation in PSI (a measure of acceptor side limitation, when P700 and A both are reduced) and energy dissipation in PSI (a measure of donor side limitation, when P700 and A both are oxidized), respectively. Recovery of PSI photochemistry under artificial conditions: (d) Y(I), (e) Y(NA), (f) Y(ND) of PSI during the recovery period. Quantum yields were calculated from 300 μmol m−2 s−1 light illumination period of a light response curve. All measurements were taken after 30 min of dark adaptation at 4 °C from winter to late spring and room temperature in summer. Error bars indicate the mean ± SD (n = 3) for (a–c); ± SD (n = 4) for (d–f). The statistically significant mean differences (t-test) are marked by asterisks indicating the following confidence intervals: * ≤95%; ** ≤99%.; *** ≤99.9%. Exact p values are provided in the source data file.

Full size image

Sustained non-photochemical quenching dominates in early spring needles
To elucidate the mechanism of “sustained quenching” observed in pine during the early spring periods, ultrafast time-resolved fluorescence measurements on intact pine needles were performed using the setup and approach described previously17,18. In addition, a special chamber was designed to keep the rotational cuvette with the sample at −20 °C (Supplementary Fig. 4). This chamber was used in order to maintain the “native” sustained quenched state in pine needles. Fluorescence decay traces were collected from pine needles in the same states as reported in the previous subsection: S, ES, ER: Summer needles (S) do not contain any “sustained quenching”, while E.spring needles (ES) contain the highest levels of “sustained quenching”. To resolve the mechanisms of this quenching, fluorescence decay traces of S needles in the original state (dark-adapted, Fig. 3a, Summer dark) were compared to the needles collected in ES when sustained quenching was present (Fig. 3a, Early spring) and after the quenching in the needles had been relaxed at room temperature (ER) (Fig. 3a, E. spring recovered). To understand whether sustained quenching observed in ES occurs via a similar quenching mechanism to light-induced quenching, S needles exposed to HL for 30 min (SQ) were also measured (Fig. 3a, Summer quenched). As shown by the direct comparison of the decay curves (Fig. 3a), the fluorescence in SQ needles is much shorter-lived as compared to that of S needles. However, the fluorescence decays of ES needles are still pronouncedly shorter-lived than that of the SQ samples, thus characterizing the “sustained quenching” in the ES samples as the most pronounced quenching at all detection wavelengths and under all conditions. In contrast, ER and S samples showed very similar fluorescence decays, indicating that the ES sample recovered quite (although not completely, vide infra) well from sustained quenching within 48 h of recovery treatment.
Fig. 3: Time-resolved fluorescence of intact pine needles measured using TCSPC.

(a) Fluorescence decay traces measured at −20 °C and at four characteristic wavelengths: 686 nm (mainly PSII, LHCII contributions), 698 nm (PSII, PSI contributions), 723 nm (mainly PSI contribution), and 741 nm (mainly PSI contribution). (b) Global analysis of pine needles in three states: Summer dark (S, dark-adapted summer needles, left row), ES (E.spring needles with “sustained quenching” present, middle row), Summer quenched (SQ, right row) (c) Kinetic target analysis of pine needles in the three states. See the main text for the explanations. The kinetic target analysis (SAS top, a kinetic model with rate constants in ns-1, bottom) shows the results of the detailed target modeling of the fluorescence kinetics of pine needles. The rate constants (ns−1) and Species-associated emission spectra (SAS) resulted were determined from global target analysis. Species-associated emission spectra (SAS) resulted from the fit of the target kinetic model in the corresponding state. Note that fluorescence decay measurements below 680 nm to detect the decreasing short-wavelength part of the spectra were not possible due to the extremely high scattering of the pine needles. This has no effect however on the ability to distinguish the various lifetime components kinetically and spectrally.

Full size image

Global analysis of the fluorescence kinetics of pine needles
To get the first hints on the mechanism(s) underlying sustained quenching in early spring global analysis19 was performed on all datasets: ES, S, SQ (Fig. 3b), and ER (Supplementary Figs. 5–8) needles. Six decay-associated spectra (DAS) were required to fit the fluorescence kinetics in all four cases. Three DAS were tentatively assigned to PSI: 11–16 ps, 50–60 ps, 150–180 ps (Fig. 3b, Supplementary Figs. 5, 6). Their spectra and lifetimes are reminiscent of previously reported PSI-related DAS in other plant species in vivo17,18,20. The fastest component represents energy equilibration between Chlred and Chlblue21,22. The second DAS in S needles peaks around 685–690 nm with a broad emission in the 700–720 nm region. Therefore, it mainly represents PSI core decay in combination with the 710–720 Chls of Lhcas. The third DAS (185 ps) peaks at 680 nm with some contribution at 730 nm, and therefore, it should be primarily assigned to LHCII and some Chlred in PSI. In quenched states (SQ and ES, in particular) these two DAS have substantially higher contribution around 680 nm and, as a result, the reconstructed steady-state PSI spectra have higher emission around 680–685 nm as compared to the unquenched states (S and ER see Supplementary Fig. 5). It indicates that upon quenching PSII and/or LHCII kinetics also contribute to these PSI-related components. In S and ER states, the remaining three DAS peak at 682 nm and therefore, are assigned to PSII (Fig. 3b and Supplementary Fig. 9). Their lifetimes are 0.5–0.6 ns, 1.5–2.0 ns, and 4 ns, similar to what was resolved previously for PSII in the closed state17,18. However, in ES two considerably shorter-lived PSII-related DAS were observed with lifetimes of 0.3 ns and 0.8 ns. This demonstrates that PSII kinetics is indeed quenched strongly in ES needles. Additionally, we resolved a DAS with 8 ps lifetime in the ES sample. The “oscillations” in the spectra giving rise to the sharp bands are due to the mixing of the contributing components, representing excitation energy transfer (EET) from Chl a pools emitting at ~680 nm to those emitting at 690 nm and 700–730 nm, respectively. This suggests that there occurs efficient EET between PSII and/or LHCII and PSI in the ES needles. However, these processes could be unambiguously resolved only by target modeling.
In a nutshell, global analysis results suggest that sustained quenching in the ES state involves both PSII/LHCII and PSI through a direct energy transfer. This mechanism provides much stronger quenching than the high-light-induced quenching in the summer needles (SQ).
Target modeling of the fluorescence kinetics of unquenched pine needles
To resolve the mechanism(s) responsible for the “sustained quenching” in pine needles in detail we performed target analysis23 of all the kinetic data. In target analysis, a specific kinetic model is tested using a linear time-invariant compartmental model. Each compartment in Fig. 3 is presented by a square box. Each box belongs to any of the three different components: PSI, PSII, LHCII quenched. Fluorescence emitting compartments are presented by colored boxes. The red flashes with numbers directing to that compartments represent the so-called excitation vector, providing the percentage of excitations directed to each compartment and should roughly illustrate the relative number of Chl a (i.e., relative Chl a antenna size) in it. Non-colored boxes represent non-emitting compartments, like radical pairs (RP) in RCs and CT states, which are contributing to the overall observable kinetics and thus have to be included. Target modeling is required to describe the real time-dependent concentrations of each compartment, rates between them (arrows in Fig. 3, ns−1), and spectra (species associated spectra (SAS)) of the interacting components. Note, that the global analysis, presented in the previous subsection, is equivalent to a number of noninteracting, parallelly decaying compartments and gives only the first indications of the processes involved. It is a mere mathematical analysis, not considering any physical description or information on the studied system. By testing kinetic models that had been validated previously on isolated subcomplexes24,25,26 it was possible to separate the PSI from the PSII kinetics. Therefore, our assignment of PSI and PSII is based on a combination of the used kinetic model, resolved rates as well as SAS. The target model and the excitation vector presented in Fig. 3c resulted in the best fitting quality (Supplementary Fig. 7). In S needles (Fig. 3c, Summer dark) the PSI kinetics was described with two emitting compartments (PSI red, Ant/core) and one non-emitting radical pair (RP, Fig. 3c). The resolved rates and SAS are strongly reminiscent of the ones reported for the PSI-LHCI complex26. Kinetics from the PSII antenna/RC complex was described by one emitting compartment (Ant/core) and two non-emitting radical pairs (RP1 and RP2, Fig. 3c) as reported previously17,24,25,27. According to their SAS and rates, the emitting species represents the equilibrated excited state of the PSII-LHCII super-complex. In our analysis, to achieve a very good fit, two pools of PSII complexes with different decay rates were required (pool 1, pool 2, Fig. 3c, see the detailed comparison of the residuals plots in Supplementary Fig. 7) for the S sample. The two pools differ in the rate constants of charge separation in a manner that is typical for PSII particles with different antenna sizes (see18,19 for discussion and explanation), which result in different average lifetimes (Supplementary Table 3). In the S sample, PSI with its antenna accounts for 30% of the total absorption cross-section at the excitation wavelength (662 nm), while the two PSII pools account for 40% and 30% of the total absorption cross-section. No internal quenching related to any NPQ process of these PSII compartments was observed.
The target model for the ER needles is rather similar to that of S needles (Supplementary Fig. 5). The main difference is that the smaller PSII pool (with a very small absorption cross-section of ca. 10%) shows an unusually small charge separation rate. The most likely explanation of this phenomenon is a small percentage of PSII RCs has been photoinhibited during the 48 h relaxation process from the sustained quenching state, which agrees with NPQ measurements.
Direct energetic transfer between PSII and PSI is the dominating component of sustained quenching
As for the S samples, also the description of the ES sample (Fig. 3c) required a three-compartmental target model for PSI and the presence of two PSII pools, each described by a three-compartmental model. The two PSII pools differ in the rate constants of charge separation in a similar manner as in S samples. The kinetics are also very similar to that of the S samples, however with a slightly increased charge separation rate. This can be interpreted as a consequence of the partial reduction in the physical antenna size. This goes in line with an increase in the Chl a/b ratio in ES as compared to S samples (Supplementary Fig. 9, Supplementary Table 2). Resolved antenna SAS (Ant/Core compartment of PSI, Fig. 3) of PSI in ES samples contains stronger emission in the 680 nm region, as compared to S needles. That is why we suggest that some LHCII strongly coupled to PSI contribute to it, like recently reported20.
The striking feature is the presence of a very strong contribution of a direct energy transfer from PSII to PSI. Without allowing for this direct energy transfer, the data could not at all be fitted adequately (Supplementary Figs. 6, 7). This process has a transfer rate of 4.4 ns−1. A small compartment representing quenched and functionally detached LHCII was also needed for a satisfactory description of the ES samples (See more details on this component in the next section). Overall, the target analysis shows that direct energy transfer from PSII to PSI complexes provides the dominant quenching mechanism of PSII complexes in ES needles causing the pronounced “sustained quenching”.
Comparison of sustained quenching and light-induced quenching
To model the fluorescence kinetics of SQ needles again three-compartmental models for PSI and PSII were required. However, unlike in previous cases, one single PSII pool was sufficient to describe the PSII kinetics.
Besides PSII and PSI, one additional component was needed to satisfactorily describe the fluorescence kinetics in the SQ samples. The rates, spectra, and average lifetime (0.4 ns, Supplementary table 3) of these compartments are strongly reminiscent of what was reported previously for highly quenched LHCII aggregates28, and the quenched LHCII complexes under NPQ conditions in wild type Arabidopsis17 and in Arabidopsis mutants depleted of the reaction centers29. We, therefore, assign this highly quenched component to the quenched LHCII antenna energetically disconnected from both PSII and PSI. However, the addition of this quenched LHCII compartment was not sufficient to fully describe the fluorescence kinetics (Supplementary Fig. 7). The species-associated spectra (SAS) and kinetics of PSII still contained a substantial amount of 690–700 and 730 nm emission, entirely uncharacteristic for PSII super-complexes, but characteristic of PSI complexes22,26. Interestingly, also the SAS of the PSI antenna/core complex contains a 680 nm emission, which is characteristic of PSII. This strongly suggests—as can be directly deduced from the spectral properties—the presence of direct energetic contact between PSI and PSII (Fig. 3c). Introducing such a connection in the target model indeed resulted in a considerably better fit (χ2 decreased from 1.198 to 1.093, Supplementary Fig. 6). Also, the PSII and PSI kinetics became well-separated in such a model. The resolved direct energy transfer rate from PSII to PSI is quite high (2.3 ns−1), implying that a part of energy absorbed by PSII is actually funneled to PSI. However, this rate of direct energy transfer in SQ needles is only about half of that in ES needles, suggesting that unlike in cold-induced sustained quenching, upon light illumination, spillover is not the major process. The contribution of the quenched antenna, on the other hand, plays a major role.
Efficiency of direct energy transfer quenching in different conditions
By evaluating the kinetic information provided in Fig. 3c one can get a quantitative estimation of the degrees of light use efficiencies (for photosynthesis) on the one hand, and for photoprotection/quenching on the other hand, present in both PSII and PSI in the different situations (c.f. Supplementary Fig. 8, Supplementary Table 4). In S needles 70% of the total absorbed energy drives PSII charge separation and only 30% PSI charge separation. In the SQ needles the percentage of energy flowing into the PSII reaction center drops to 7.1%, while it increases to 67% for PSI. Upon strong light illumination, most of the PSII reaction centers will be closed, and a fraction of that total energy flowing through the RC can produce harmful reactive species like ROS30,31. Since the percentage of energy flow into PSII decreases 10 times (70% in S and 7.1% in SQ) in SQ as compared to S, this means that the potential of oxidative damage is reduced by a factor of 10 by the quenching process. This factor is increased by a further factor of approx. 2 due to the large pool of functionally detached and quenched LHCII. So, in SQ needles the overall quenching provides a protection factor of 20 for PSII compared to S needles. In the sustained quenched ES samples this effect is much more extreme: Only 1.5% of the total absorbed energy in the system flows through the PSII reaction center potentially causing oxidative damage (assuming all other factors are equal). Thus, in total, “sustained quenching” in ES samples provides a protection factor of approx. 40 to PSII (double of SQ) compared to S needles. This huge protection factor explains well why “sustained quenching” is so effective enabling pine needles to survive the harsh winter and early spring conditions.
Massive destacking of grana membranes in the winter
In accordance with early reports32,33, we observed strong seasonal changes in chloroplast ultrastructure, including loss of grana stacks during early spring (Fig. 4a). Ultrastructural changes derived from morphometric analysis of electron micrographs (Fig. 4a) indicate that average numbers of grana per chloroplast steadily decreased from autumn to winter and reached minimal values in early spring, corresponding to just 68% of S values (Supplementary Table 1). More strikingly, the number of appressed thylakoid layers/granum dramatically declined from A (4.97) to ES (2.72) and rose back in S (6.50) (Supplementary Table 1). These changes corroborate the shift in grana stacking illustrated in Fig. 4b. In S and A, 3–6 layered grana stacks accounted for 60–70% of total stacks. In contrast in ES samples, 70–75% and 15–20% of the grana had only two and three layers, respectively (Fig. 4b). These changes in ES samples were accompanied by a transient doubling of the number of lipid globules (plastoglobuli) per cell (Supplementary Table 1) when the chloroplast structure deviated most strongly from the S state.
Fig. 4: Transmission electron microphotographs depicting seasonal variations in chloroplast ultrastructure in pine needles.

a Chloroplast structure in Autumn, Winter, E. Spring, L. Spring, and Summer. b Histograms of frequency distributions of numbers of thylakoids per granum during the five distinct seasonal periods. The histograms were calculated from 80–100 electron micrographs per season, representing 2–3 chloroplasts per image [Error bars indicate mean ± SD calculated from n = 1124 (Autumn), 1235 (Winter), 521 (Early Spring), 578 (Late Spring) and 1048 (Summer) of grana stacks].

Full size image

To obtain a deeper understanding of the thylakoid plasticity in ES needles we also subjected needles collected at specific early spring and summer dates to EM analysis after artificially induced recovery under different light conditions (Fig. 5a). When both S and ER samples were exposed to HL for 30 (SQ1 for S samples, ERQ1 for ER samples) and 60 min (SQ2 for S samples, ERQ2 for ER samples) similar destacking of grana (as of ES samples) have been observed (Fig. 5b). These changes in grana structures recorded in ER and S samples (Fig. 5c, d) following short term HL exposure strongly corroborate the high dynamic plasticity of the thylakoid membrane, which appears to be essential for acclimation to the harsh boreal winters.
Fig. 5: Artificial induction of changes in chloroplast ultrastructure of pine needles.

a Changes in chloroplast ultrastructure in E. spring (ES), E. spring samples recovered (ER) at 18oC for 48 h with a photoperiod of 18 h at 80 µmol m−2 s−1, ER samples treated with 800 µmol  m−2 s−1 high light for 30 min (ERQ1), for 60 min (ERQ2). Summer (S), Summer samples treated with 1200 µmol m−2 s−1 high light for 30 min (SQ1), for 60 min (SQ2). b The number of grana per chloroplasts (Error bars indicate mean ± SD (n = 75); c Histograms of frequency distributions of numbers of thylakoids per granum in different E. spring treated [n = 220 (ES), 250 (ER), 272 (ERQ1), 246 (ERQ2)]. d Summer treated [n = 576 (S), 498 (SQ1), 415 (SQ2)] samples. Error bars indicate the mean ± SD obtained from the analysis of grana stacks.

Full size image

Changes in pigment and protein composition of ES, ER, SQ, and S samples
We also analyzed pigments of the pine needles used for spectroscopic measurements and found substantial differences in pigment contents between S and ES needles (Supplementary Table 2). In S and ES samples the Chl a/b ratios were 2.8 and 3.4, respectively, indicating that S needles had higher relative amounts of Lhcbs and other Chl b-binding proteins. The Chl/carotenoid ratio was also substantially lower in ES samples than in S samples, indicating that they had higher amounts of carotenoids (mainly lutein and violaxanthin/zeaxanthin). We observed no significant difference in total Chl pigment contents between S and SQ samples.
For quantifying the core and antenna proteins, isolated thylakoids were subjected to SDS-PAGE and were immunoblotted against specific antibodies (Supplementary Fig. 9). This revealed that ES samples contain a lower abundance of core (PsaD, PsbD) and antenna (Lhcb2, Lhca4) proteins than in S, in line with earlier reports34 (Supplementary Fig. 9). ER samples contained lower amounts of PsbD and Lhcb2 compared to S, which might reflect on the minor amount of photoinhibition seen in fluorescence data. More

1.
Barrat, A., Barthélemy, M., Pastor-Satorras, R. & Vespignani, A. The architecture of complex weighted networks. Proc. Natl. Acad. Sci. USA 101, 3747–3752 (2004).
ADS  CAS  PubMed  MATH  Article  Google Scholar 
2.
Colizza, V., Barthélemy, M., Barrat, A. & Vespignani, A. Epidemic modeling in complex realities. Comptes Rendus Biol. 330, 364–374 (2007).
Article  Google Scholar 

3.
Craft, M. E., Volz, E., Packer, C. & Meyers, L. A. Distinguishing epidemic waves from disease spillover in a wildlife population. Proc. R. Soc. B Biol. Sci. 276, 1777–1785 (2009).
Article  Google Scholar 

4.
Vernon, M. C. & Keeling, M. J. Representing the UK’s cattle herd as static and dynamic networks. Proc. Biol. Sci. 276, 469–476 (2009).
PubMed  Google Scholar 

5.
Kao, R. R., Green, D. M., Johnson, J. & Kiss, I. Z. Disease dynamics over very different time-scales: foot-and-mouth disease and scrapie on the network of livestock movements in the UK. J. R. Soc. Interface 4, 907–916 (2007).
PubMed  PubMed Central  Article  Google Scholar 

6.
Craft, M. E. Infectious disease transmission and contact networks in wildlife and livestock. Philos. Trans. R. Soc. B Biol. Sci. 370, 20140107–20140107 (2015).
Article  Google Scholar 

7.
Chiner-Oms, Á. & Comas, I. Large genomics datasets shed light on the evolution of the Mycobacterium tuberculosis complex. Infect. Genet. Evol. https://doi.org/10.1016/j.meegid.2019.02.028 (2019).
Article  PubMed  Google Scholar 

8.
Köser, C. U. et al. Rapid whole-genome sequencing for investigation of a neonatal MRSA outbreak. N. Engl. J. Med. 366, 2267–2275 (2012).
PubMed  PubMed Central  Article  Google Scholar 

9.
Roetzer, A. et al. Whole genome sequencing versus traditional genotyping for investigation of a Mycobacterium tuberculosis outbreak: a longitudinal molecular epidemiological study. PLoS Med. 10, e1001387 (2013).
PubMed  PubMed Central  Article  Google Scholar 

10.
Kao, R. R., Haydon, D. T., Lycett, S. J. & Murcia, P. R. Supersize me: how whole-genome sequencing and big data are transforming epidemiology. Trends Microbiol. 22, 282–291 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

11.
Wymant, C. et al. PHYLOSCANNER: Inferring transmission from within- and between-host pathogen genetic diversity. Mol. Biol. Evol. 35, 719–733 (2018).
CAS  PubMed  Article  Google Scholar 

12.
Gutiérrez, S., Michalakis, Y. & Blanc, S. Virus population bottlenecks during within-host progression and host-to-host transmission. Curr. Opin. Virol. 2, 546–555 (2012).
PubMed  Article  CAS  Google Scholar 

13.
Buckee, C. O. F., Koelle, K., Mustard, M. J. & Gupta, S. The effects of host contact network structure on pathogen diversity and strain structure. Proc. Natl. Acad. Sci. USA 101, 10839–10844 (2004).
ADS  CAS  PubMed  Article  Google Scholar 

14.
Rasmussen, D. A., Ratmann, O. & Koelle, K. Inference for nonlinear epidemiological models using genealogies and time series. PLoS Comput. Biol. 7, e1002136 (2011).
ADS  MathSciNet  CAS  PubMed  PubMed Central  Article  Google Scholar 

15.
Rasmussen, D. A., Volz, E. M. & Koelle, K. Phylodynamic inference for structured epidemiological models. PLoS Comput. Biol. 10, e1003570 (2014).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

16.
Rasmussen, D. A., Kouyos, R., Günthard, H. F. & Stadler, T. Phylodynamics on local sexual contact networks. PLoS Comput. Biol. 13, 1–23 (2017).
Article  CAS  Google Scholar 

17.
Cottam, E. M. et al. Integrating genetic and epidemiological data to determine transmission pathways of foot-and-mouth disease virus. Proc. R. Soc. B Biol. Sci. 275, 887–895 (2008).
Article  Google Scholar 

18.
Ypma, R. J. F. et al. Unravelling transmission trees of infectious diseases by combining genetic and epidemiological data. Proc. R. Soc. B Biol. Sci. 279, 444–450 (2012).
CAS  Article  Google Scholar 

19.
Ypma, R. J. F., van Ballegooijen, W. M. & Wallinga, J. Relating phylogenetic trees to transmission trees of infectious disease outbreaks. Genetics 195, 1055–1062 (2013).
PubMed  PubMed Central  Article  Google Scholar 

20.
Morelli, M. J. et al. A Bayesian inference framework to reconstruct transmission trees using epidemiological and genetic data. PLoS Comput. Biol. 8, e1002768 (2012).
MathSciNet  CAS  PubMed  PubMed Central  Article  Google Scholar 

21.
Lau, M. S. Y., Marion, G., Streftaris, G. & Gibson, G. A systematic Bayesian integration of epidemiological and genetic data. PLoS Comput. Biol. 11, 1–27 (2015).
Article  CAS  Google Scholar 

22.
De Maio, N., Wu, C. H. & Wilson, D. J. SCOTTI: efficient reconstruction of transmission within outbreaks with the structured coalescent. PLoS Comput. Biol. 12, 1–23 (2016).
Google Scholar 

23.
Li, L. M., Grassly, N. C. & Fraser, C. Quantifying transmission heterogeneity using both pathogen phylogenies and incidence time series. Mol. Biol. Evol. 34, 2982–2995 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

24.
Romero-Severson, E. O., Bulla, I. & Leitner, T. Phylogenetically resolving epidemiologic linkage. Proc. Natl. Acad. Sci. USA 113, 2690–2695 (2016).
ADS  CAS  PubMed  Article  Google Scholar 

25.
Firestone, S. M. et al. Reconstructing foot-and-mouth disease outbreaks: a methods comparison of transmission network models. Sci. Rep. 9, 4809 (2019).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

26.
Campbell, F., Strang, C., Ferguson, N., Cori, A. & Jombart, T. When are pathogen genome sequences informative of transmission events?. PLoS Pathog. 14, 1–21 (2018).
Google Scholar 

27.
Allen, A. R. One bacillus to rule them all? Investigating broad range host adaptation in Mycobacterium bovis. Infect. Genet. Evol. 53, 68–76 (2017).
PubMed  Article  Google Scholar 

28.
Biek, R. et al. Whole genome sequencing reveals local transmission patterns of Mycobacterium bovis in sympatric cattle and badger populations. PLoS Pathog. 8, e1003008 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

29.
Glaser, L. et al. Descriptive epidemiology and whole genome sequencing analysis for an outbreak of bovine tuberculosis in beef cattle and white-tailed deer in northwestern Minnesota. PLoS ONE 11, e0145735 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

30.
Orloski, K., Robbe-Austerman, S., Stuber, T., Hench, B. & Schoenbaum, M. Whole genome sequencing of Mycobacterium bovis isolated from livestock in the United States, 1989–2018. Front. Vet. Sci. 5, 1–10 (2018).
Article  Google Scholar 

31.
Palmer, M. V. Mycobacterium bovis: characteristics of wildlife reservoir hosts. Transbound. Emerg. Dis. 60, 1–13 (2013).
PubMed  Article  Google Scholar 

32.
Keeling, M. J. & Ross, J. V. On methods for studying stochastic disease dynamics. J. R. Soc. Interface 5, 171–181 (2008).
CAS  PubMed  Article  Google Scholar 

33.
Sharkey, K. J. Deterministic epidemiological models at the individual level. J. Math. Biol. 57, 311–331 (2008).
MathSciNet  PubMed  MATH  Article  Google Scholar 

34.
Stollenwerk, N. & Jansen, V. A. A. Meningitis, pathogenicity near criticality: the epidemiology of meningococcal disease as a model for accidental pathogens. J. Theor. Biol. 222, 347–359 (2003).
PubMed  MATH  Article  Google Scholar 

35.
Delahay, R. J. et al. Long-term temporal trends and estimated transmission rates for Mycobacterium bovis infection in an undisturbed high-density badger (Meles meles) population. Epidemiol. Infect. 141, 1445–1456 (2013).
CAS  PubMed  Article  Google Scholar 

36.
Crispell, J. et al. Combining genomics and epidemiology to analyse bi-directional transmission of Mycobacterium bovis in a multi-host system. Elife https://doi.org/10.7554/eLife.45833.001 (2019).
Article  PubMed  PubMed Central  Google Scholar 

37.
Ford, C. B. et al. Use of whole genome sequencing to estimate the mutation rate of Mycobacterium tuberculosis during latent infection. Nat. Genet. 43, 482–488 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

38.
Colangeli, R. et al. Whole genome sequencing of Mycobacterium tuberculosis reveals slow growth and low mutation rates during latent infections in humans. PLoS ONE 9, 1–9 (2014).
MathSciNet  Article  CAS  Google Scholar 

39.
Didelot, X., Fraser, C., Gardy, J., Colijn, C. & Malik, H. Genomic infectious disease epidemiology in partially sampled and ongoing outbreaks. Mol. Biol. Evol. 34, 997–1007 (2017).
CAS  PubMed  PubMed Central  Google Scholar 

40.
Lycett, S. J. et al. Role for migratory wild birds in the global spread of avian influenza H5N8. Science 354, 213–217 (2016).
Article  CAS  Google Scholar 

41.
Saulnier, E., Gascuel, O. & Alizon, S. Inferring epidemiological parameters from phylogenies using regression-ABC: a comparative study. PLoS Comput. Biol. 13, 1–31 (2017).
Article  CAS  Google Scholar 

42.
De Maio, N., Worby, C. J., Wilson, D. J. & Stoesser, N. Bayesian reconstruction of transmission within outbreaks using genomic variants. PLoS Comput. Biol. 14, 1–23 (2018).
Google Scholar 

43.
Pybus, O. G. & Rambaut, A. Evolutionary analysis of the dynamics of viral infectious disease. Nat. Rev. Genet. 10, 540–550 (2009).
CAS  PubMed  PubMed Central  Article  Google Scholar 

44.
Biek, R., Pybus, O. G., Lloyd-Smith, J. O. & Didelot, X. Measurably evolving pathogens in the genomic era. Trends Ecol. Evol. 30, 306–313 (2015).
PubMed  PubMed Central  Article  Google Scholar 

45.
Brooks-Pollock, E., Roberts, G. O. & Keeling, M. J. A dynamic model of bovine tuberculosis spread and control in Great Britain. Nature 511, 228–231 (2014).
ADS  CAS  PubMed  Article  Google Scholar 

46.
Campbell, F., Cori, A., Ferguson, N. & Jombart, T. Bayesian inference of transmission chains using timing of symptoms, pathogen genomes and contact data. PLOS Comput. Biol. 15, e1006930 (2019).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

47.
Robertson, A., Palphramand, K. L., Carter, S. P. & Delahay, R. J. Group size correlates with territory size in European badgers: implications for the resource dispersion hypothesis?. Oikos 124, 507–514 (2015).
Article  Google Scholar 

48.
Roper, T. Badger (Collins, London, 2010).
Google Scholar 

49.
Drewe, J. A., Tomlinson, A. J., Walker, N. J. & Delahay, R. J. Diagnostic accuracy and optimal use of three tests for tuberculosis in live badgers. PLoS ONE 5, e11196 (2010).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

50.
Walker, T. M. et al. Whole-genome sequencing to delineate Mycobacterium tuberculosis outbreaks: a retrospective observational study. Lancet Infect. Dis. 13, 137–146 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

51.
Álvarez, J. et al. Research in Veterinary Science Bovine tuberculosis : within-herd transmission models to support and direct the decision-making process. Res. Vet. Sci. 97, S61–S68 (2014).
PubMed  Article  Google Scholar 

52.
Rossi, G. et al. Epidemiological modelling for the assessment of bovine tuberculosis surveillance in the dairy farm network in Emilia-Romagna (Italy). Epidemics 11, 62–70 (2015).
PubMed  Article  Google Scholar 

53.
Kao, R. R., Roberts, M. G. & Ryan, T. J. A model of bovine tuberculosis control in domesticated cattle herds. Proc. R. Soc. B Biol. Sci. 264, 1069–1076 (1997).
ADS  CAS  Article  Google Scholar 

54.
Conlan, A. J. K. et al. Estimating the hidden burden of bovine tuberculosis in Great Britain. PLoS Comput. Biol. 8, e1002730 (2012).
MathSciNet  CAS  PubMed  PubMed Central  Article  Google Scholar 

55.
O’Hare, A., Orton, R. J., Bessell, P. R. & Kao, R. R. Estimating epidemiological parameters for bovine tuberculosis in British cattle using a Bayesian partial-likelihood approach. Proc. R. Soc. B Biol. Sci. 281, 20140248 (2014).
Article  Google Scholar 

56.
Rossi, G., Aubry, P., Dubé, C. & Smith, R. L. The spread of bovine tuberculosis in Canadian shared pastures: data, model, and simulations. Transbound. Emerg. Dis. 66, 562–577 (2019).
PubMed  Article  Google Scholar 

57.
R Core Team. R: A Language and Environment for Statistical Computing. (2018).

58.
Soetaert, K., Petzoldt, T. & Setzer, R. W. Solving differential equations in R: package deSolve. J. Stat. Softw. 33, 1–25 (2010).
Google Scholar 

59.
Lawes, J. R. et al. Bovine TB surveillance in Great Britain in 2014. Vet. Rec. 178, 310–315 (2016).
CAS  PubMed  Article  Google Scholar 

60.
Trewby, H. et al. Use of bacterial whole-genome sequencing to investigate local persistence and spread in bovine tuberculosis. Epidemics 14, 26–35 (2016).
PubMed  PubMed Central  Article  Google Scholar 

61.
Crispell, J. et al. Using whole genome sequencing to investigate transmission in a multi-host system: bovine tuberculosis in New Zealand. BMC Genom. 18, 180 (2017).
Article  CAS  Google Scholar 

62.
Salvador, L. C. M. et al. Disease management at the wildlife-livestock interface: using whole-genome sequencing to study the role of elk in Mycobacterium bovis transmission in Michigan. USA. Mol. Ecol. https://doi.org/10.1111/mec.15061 (2019).
Article  PubMed  Google Scholar 

63.
Brooks-Pollock, E. et al. Age-dependent patterns of bovine tuberculosis in cattle. Vet. Res. 44, 1 (2013).
Article  Google Scholar  More

1.
Dickman, A. J. & Hazzah, L. Money, myths and man-eaters: Complexities of human–wildlife conflict. In Problematic Wildlife (ed. Angelici, F. M.) 339–356 (Springer, Berlin, 2016). https://doi.org/10.1007/978-3-319-22246-2_16.
Google Scholar 
2.
Nyhus, P. J. Human–wildlife conflict and coexistence. Annu. Rev. Environ. Resour. 41, 143–171 (2016).
Article  Google Scholar 

3.
Carter, N. H. et al. Coupled human and natural systems approach to wildlife research and conservation. Ecol. Soc. 19, 43 (2014).
Article  Google Scholar 

4.
Dirzo, R. et al. Defaunation in the anthropocene. Science (80–) 345, 401–406 (2014).
ADS  CAS  Article  Google Scholar 

5.
Balasubramaniam, K. N. et al. Addressing the challenges of human–wildlife conflict research using coupled natural and human systems. Biol. Conserv. Rev.

6.
McDonald, A. M. H., Rea, R. V. & Hesse, G. Perceptions of moose-human conflicts in an urban environment. ALCES 48, 123–130 (2012).
Google Scholar 

7.
Berger-Tal, O. et al. A systematic survey of the integration of animal behavior into conservation. Conserv. Biol. 30, 744–753 (2016).
PubMed  Article  PubMed Central  Google Scholar 

8.
Snijders, L., Blumstein, D. T., Stanley, C. R. & Franks, D. W. Animal social network theory can help wildlife conservation. Trends Ecol. Evol. 32, 567–577 (2017).
PubMed  Article  PubMed Central  Google Scholar 

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

10.
Lischka, S. A. et al. A conceptual model for the integration of social and ecological information to understand human–wildlife interactions. Biol. Conserv. 225, 80–87 (2018).
Article  Google Scholar 

11.
Morrow, K. S., Glanz, H., Ngakan, P. O. & Riley, E. P. Interactions with humans are jointly influenced by life history stage and social network factors and reduce group cohesion in moor macaques (Macaca maura). Sci. Rep. 9, 1–12 (2019).
Article  CAS  Google Scholar 

12.
Dickman, A. J. From cheetahs to chimpanzees: A comparative review of the drivers of human-carnivore conflict and human-primate conflict. Folia Primatol. 83, 377–387 (2013).
Article  Google Scholar 

13.
Chiyo, P. I., Moss, C. J. & Alberts, S. C. The influence of life history milestones and association networks on crop-raiding behavior in male african elephants. PLoS One 7, 20 (2012).
Google Scholar 

14.
Maréchal, L. et al. Impacts of tourism on anxiety and physiological stress levels in wild male Barbary macaques. Biol. Conserv. 144, 2188–2193 (2011).
Article  Google Scholar 

15.
Barua, M., Bhagwat, S. A. & Jadhav, S. The hidden dimensions of human–wildlife conflict: Health impacts, opportunity and transaction costs. Biol. Cons. 157, 309–316 (2013).
Article  Google Scholar 

16.
Devaux, C. A., Mediannikov, O., Medkour, H. & Raoult, D. Infectious disease risk across the growing human–non human primate interface: A review of the evidence. Front. Public Health 7, 1–22 (2019).
Article  Google Scholar 

17.
Wolfe, N. D., Dunavan, C. P. & Diamond, J. Origins of major human infectious diseases. Nature 447, 279–283 (2007).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

18.
Oro, D., Genovart, M., Tavecchia, G., Fowler, M. S. & Martínez-Abraín, A. Ecological and evolutionary implications of food subsidies from humans. Ecol. Lett. 16, 1501–1514 (2013).
PubMed  Article  Google Scholar 

19.
LaBarge, L. R., Hill, R. A., Berman, C. M., Margulis, S. W. & Allan, A. T. L. Anthropogenic influences on primate antipredator behavior and implications for research and conservation. Am. J. Primatol. 82, 1–13 (2020).
Article  Google Scholar 

20.
Wong, B. B. M. & Candolin, U. Behavioral responses to changing environments. Behav. Ecol. 26, 665–673 (2015).
Article  Google Scholar 

21.
Hockings, K. J., Anderson, J. R. & Matsuzawa, T. Socioecological adaptations by chimpanzees, Pan troglodytes verus, inhabiting an anthropogenically impacted habitat. Anim. Behav. 83, 801–810 (2012).
Article  Google Scholar 

22.
Carne, C., Semple, S., MacLarnon, A., Majolo, B. & Maréchal, L. Implications of tourist–macaque interactions for disease transmission. EcoHealth 14, 704–717 (2017).
PubMed  PubMed Central  Article  Google Scholar 

23.
Morzillo, A. T., de Beurs, K. M. & Martin-Mikle, C. J. A conceptual framework to evaluate human–wildlife interactions within coupled human and natural systems. Ecol. Soc. 19, 44 (2014).
Article  Google Scholar 

24.
Balasubramaniam, K. N. et al. Using a coupled natural and human systems approach to understand the behavioral dimensions of human–primate interfaces.

25.
Fuentes, A. Ethnoprimatology and the anthropology of the human–primate interface. Annu. Rev. 41, 101–117 (2012).
Google Scholar 

26.
Radhakrishna, S. & Sinha, A. Less than wild? Commensal primates and wildlife conservation. J. Biosci. 36, 749–753 (2011).
PubMed  Article  Google Scholar 

27.
Soulsbury, C. D. & White, P. C. L. Human–wildlife interactions in urban areas: A review of conflicts, benefits and opportunities. Wildl. Res. 42, 541 (2015).
Article  Google Scholar 

28.
Fuentes, A. & Gamerl, S. Disproportionate participation by age/sex classes in aggressive interactions between long-tailed macaques (Macaca fascicularis) and human tourists at Padangtegal monkey forest, Bali, Indonesia. Am. J. Primatol. 66, 197–204 (2005).
PubMed  Article  Google Scholar 

29.
Hsu, M. J., Kao, C. C. & Agoramoorthy, G. Interactions between visitors and formosan macaques (Macaca cyclopis) at Shou-Shan Nature Park, Taiwan. Am. J. Primatol. 71, 214–222 (2009).
PubMed  Article  Google Scholar 

30.
Wolf, M., van Doorn, G. S., Leimar, O. & Weissing, F. J. Life-history trade-offs favour the evolution of animal personalities. Nature 447, 581–584 (2007).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

31.
Promislow, D. E. L. & Harvey, P. H. Living fast and dying young: A comparative analysis of life-history variation among mammals. J. Zool. 220, 417–437 (1990).
Article  Google Scholar 

32.
Clutton-Brock, T. Reproductive competition and sexual selection. Philos. Trans. R. Soc. B Biol. Sci. 372, 20160310 (2017).
Article  Google Scholar 

33.
Mitani, J. C., Gros-Louis, J. & Richards, A. F. Sexual dimorphism, the operational sex ratio, and the intensity of male competition in polygynous primates. Am. Nat. 147, 966–980 (1996).
Article  Google Scholar 

34.
Marty, P. R. et al. Individuals in urban dwelling primate species face unequal benefits associated with living in an anthropogenic environment. Rev. (2019).

35.
Fehlmann, G. et al. Extreme behavioural shifts by baboons exploiting risky, resource-rich, human-modified environments. Sci. Rep. 7, 1–8 (2017).
CAS  Article  Google Scholar 

36.
Beisner, B. A. et al. Human–wildlife conflict: Proximate predictors of aggression between humans and rhesus macaques in India. Am. J. Phys. Anthropol. 156, 286–294 (2015).
PubMed  Article  PubMed Central  Google Scholar 

37.
Fa, J. E. Visitor-directed aggression among the Gibraltar macaques. Zoo Biol. 11, 43–52 (1992).
Article  Google Scholar 

38.
Bruintjes, R. & Radford, A. N. Context-dependent impacts of anthropogenic noise on individual and social behaviour in a cooperatively breeding fish. Anim. Behav. 85, 1343–1349 (2013).
Article  Google Scholar 

39.
Lacy, K. E. & Martins, E. P. The effect of anthropogenic habitat usage on the social behaviour of a vulnerable species, Cyclura nubila. Anim. Conserv. 6, 3–9 (2003).
Article  Google Scholar 

40.
Ram, S., Venkatachalam, S. & Sinha, A. Changing social strategies of wild female bonnet macaques during natural foraging and on provisioning. Curr. Sci. 84, 780–790 (2003).
Google Scholar 

41.
Maréchal, L., MacLarnon, A., Majolo, B. & Semple, S. Primates’ behavioural responses to tourists: Evidence for a trade-off between potential risks and benefits. Sci. Rep. 6, 1–11 (2016).
Article  CAS  Google Scholar 

42.
Kaburu, S. S. K. et al. Interactions with humans impose time constraints on urban-dwelling rhesus macaques (Macaca mulatta). Behaviour 156, 1255–1282 (2019).
Article  Google Scholar 

43.
Marty, P. R. et al. Time constraints imposed by anthropogenic environments alter social behaviour in long-tailed macaques. Anim. Behav. 150, 157–165 (2019).
Article  Google Scholar 

44.
Balasubramaniam, K. N. et al. Impact of anthropogenic factors on affiliative behaviors among bonnet macaques. Am. J. Phys. Anthropol. 171, 704–717 (2020).
PubMed  Article  Google Scholar 

45.
Belton, L. E., Cameron, E. Z. & Dalerum, F. Social networks of spotted hyenas in areas of contrasting human activity and infrastructure. Anim. Behav. 135, 13–23 (2018).
Article  Google Scholar 

46.
Chilvers, B. L. & Corkeron, P. J. Trawling and bottlenose dolphins’ social structure. Proc. R. Soc. B 268, 1901–1905 (2001).
CAS  PubMed  Article  Google Scholar 

47.
Jacobs, A., Sueur, C., Deneubourg, J. L. & Petit, O. Social network influences decision making during collective movements in Brown lemurs (Eulemur fulvus fulvus). Int. J. Primatol. 32, 721–736 (2011).
Article  Google Scholar 

48.
Bode, N. W. F., Wood, A. J. & Franks, D. W. The impact of social networks on animal collective motion. Anim. Behav. 82, 29–38 (2011).
Article  Google Scholar 

49.
Sueur, C., MacIntosh, A. J. J., Jacobs, A. T., Watanabe, K. & Petit, O. Predicting leadership using nutrient requirements and dominance rank of group members. Behav. Ecol. Sociobiol. 67, 457–470 (2013).
Article  Google Scholar 

50.
Schino, G. Grooming, competition and social rank among female primates: A meta-analysis. Anim. Behav. 62, 265–271 (2001).
Article  Google Scholar 

51.
Balasubramaniam, K. N. & Berman, C. M. Grooming interchange for resource tolerance: Biological markets principles within a group of free-ranging rhesus macaques. Behaviour 154, 1145–1176 (2017).
Article  Google Scholar 

52.
Overduin-de Vries, A. M., de Vries, H., Vermande, M. M., Reijntjes, A. H. A. & Sterck, E. H. M. Both aggressive and affiliative behaviour facilitate resource access in high-ranking female long-tailed macaques (Macaca fascicularis). Behaviour 157, 267–287 (2020).
Article  Google Scholar 

53.
Kurita, H., Sugiyama, Y., Ohsawa, H., Hamada, Y. & Watanabe, T. Changes in demographic parameters of Macaca fuscata at Takasakiyama in relation to decrease of provisioned foods. Int. J. Primatol. 29, 1189–1202 (2008).
Article  Google Scholar 

54.
Bonacich, P. Some unique properties of eigenvector centrality. Soc. Netw. 29, 555–564 (2007).
Article  Google Scholar 

55.
Wey, T., Blumstein, D. T., Shen, W. & Jordán, F. Social network analysis of animal behaviour: A promising tool for the study of sociality. Anim. Behav. 75, 333–344 (2008).
Article  Google Scholar 

56.
Dunbar, R. I. M. Time: A hidden constraint on the behavioural ecology of baboons. Behav. Ecol. Sociobiol. 31, 35–49 (1992).
Article  Google Scholar 

57.
Sinha, A., Mukhopadhyay, K., Datta-Roy, A. & Ram, S. Ecology proposes, behaviour disposes: Ecological variability in social organization and male behavioural strategies among wild bonnet macaques. Curr. Sci. 89, 1166–1179 (2005).
Google Scholar 

58.
Thierry, B. Unity in diversity: Lessons from macaque societies. Evol. Anthropol. 16, 224–238 (2007).
Article  Google Scholar 

59.
Balasubramaniam, K. N. et al. Hierarchical steepness and phylogenetic models: Phylogenetic signals in Macaca. Anim. Behav. 83, 20 (2012).
Article  Google Scholar 

60.
van Schaik, C. P. The ecology of social relationships amongst female primates. In Comparative Socioecology (eds Standen, V. & Foley, R. A.) 195–218 (Blackwell, New York, 1989).
Google Scholar 

61.
Morales, J. C. & Melnick, D. J. Phylogenetic relationships of the macaques (Cercopithecidae: Macaca), as revealed by high resolution restriction site mapping of mitochondrial ribosomal genes. J. Hum. Evol. 34, 1–23 (1998).
CAS  PubMed  Article  Google Scholar 

62.
Gumert, M. D. A common monkey of Southeast Asia: Longtailed macaque populations, ethnophoresy, and their occurrence in human environments. In Monkeys on the Edge: Ecology and Management of longTailed Macaques and their Interface with Humans (eds Gumert, M. D. et al.) 3–43 (Cambridge University Press, Cambridge, 2011).
Google Scholar 

63.
Southwick, C. H. & Siddiqi, M. F. Primate commensalism: The Rhesus Monkey in India. Rev. Ecol. (Terre Vie) 49, 223–231 (1994).
Google Scholar 

64.
Sinha, A. & Mukhopadhyay, K. The monkey in the Town’s commons, revisited: An anthropogenic history of the Indian Bonnet Macaque. In The Macaque Connection: Cooperation and Conflict between Humans and Macaques (eds Radhakrishna, S. et al.) 187–208 (Springer, Berlin, 2013). https://doi.org/10.1007/978-1-4614-3967-7.
Google Scholar 

65.
Kaburu, S. S. K. et al. Rates of human–macaque interactions affect grooming behavior among urban-dwelling rhesus macaques (Macaca mulatta). Am. J. Phys. Anthropol. 168, 92–103 (2018).
PubMed  Article  Google Scholar 

66.
Altmann, J. Observational study of behavior: Sampling methods. Behaviour 49, 227–267 (1974).
CAS  Article  Google Scholar 

67.
Martin, P. & Bateson, P. Measuring Behaviour (Cambridge University Press, Cambridge, 1993).
Google Scholar 

68.
R Development Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, 2009).
Google Scholar 

69.
Fujii, K. et al. Perc: Using Percolation and Conductance to Find Information Flow Certainty in a Direct Network. (R Package Version 0.1., 2015).

70.
Handcock, M. S., Hunter, D. R., Butts, C. T., Goodreau, S. M. & Morris, M. Package ‘statnet’. R Packag. (2008).

71.
Shannon, P. et al. Cytoscape: A software environment for integrated models. Genome Res. 13, 2498–2504 (2003).
CAS  PubMed  PubMed Central  Article  Google Scholar 

72.
Burnham, K. P. & Anderson, D. R. Model Selection and Multimodel Inference (Springer, Berlin, 2002).
Google Scholar 

73.
Burnham, K. P., Anderson, D. R. & Huyvaert, K. P. AIC model selection and multimodel inference in behavioral ecology: Some background, observations, and comparisons. Behav. Ecol. Sociobiol. 65, 23–35 (2011).
Article  Google Scholar 

74.
Richards, S. A. Testing ecological theory using the information-theoretic approach: Examples and cautionary results. Ecology 86, 2805–2814 (2005).
Article  Google Scholar 

75.
Quinn, G. P. & Keough, M. J. Experimental Designs and Data Analysis for Biologists (Cambridge University Press, Cambridge, 2002).
Google Scholar 

76.
Bates, D., Mächler, M., Bolker, B. M. & Walker, S. C. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 20 (2015).
Article  Google Scholar 

77.
Bolker, B., Skaug, H., Magnusson, A. & Nielsen, A. Getting started with the glmmADMB package. R Packag. ver. 2.0–8 12 (2012).

78.
Mckinney, T. A classification system for describing anthropogenic influence on nonhuman primate populations. Am. J. Primatol. 77, 715–726 (2015).
PubMed  Article  PubMed Central  Google Scholar 

79.
Silk, J. B. Kin selection in primate groups. Int. J. Primatol. 23, 849–875 (2002).
Article  Google Scholar 

80.
Marty, P. R., Hodges, K., Agil, M. & Engelhardt, A. Alpha male replacements and delayed dispersal in crested macaques (Macaca nigra). Am. J. Primatol. 79, 1–8 (2017).
Article  Google Scholar 

81.
Janson, C. H. Social correlates of individual spatial choice in foraging groups of brown capuchin monkeys, Cebus apella. Anim. Behav. 40, 910–921 (1990).
Article  Google Scholar 

82.
Hemelrijk, C. K. Towards the integration of social dominance and spatial structure. Anim. Behav. 59, 1035–1048 (2000).
CAS  PubMed  Article  Google Scholar 

83.
VanderWaal, K. L., Atwill, E. R., Isbell, L. A. & McCowan, B. Linking social and pathogen transmission networks using microbial genetics in giraffe (Giraffa camelopardalis). J. Anim. Ecol. 83, 406–414 (2014).
PubMed  Article  Google Scholar 

84.
Lusseau, D. et al. Quantifying the influence of sociality on population structure in bottlenose dolphins. J. Anim. Ecol. 75, 14–24 (2006).
PubMed  Article  Google Scholar 

85.
Pfeiffer, M. B., Iglay, R. B., Seamans, T. W., Blackwell, B. F. & DeVault, T. L. Deciphering interactions between white-tailed deer and approaching vehicles. Transp. Res. Part D 79, 102251 (2020).
Article  Google Scholar 

86.
Hamilton, W. D. Geometry for the selfish herd. J. Theor. Biol. 31, 295–311 (1971).
CAS  PubMed  Article  Google Scholar 

87.
Hall, C. L. & Fedigan, L. M. Spatial benefits afforded by high rank in white-faced capuchins. Anim. Behav. 53, 1069–1082 (1997).
Article  Google Scholar 

88.
Johnson, M. L. Exploratory behavior and dispersal: A graphical model. Can. J. Zool. 67, 2325–2328 (1989).
Article  Google Scholar 

89.
Marty, P. R. Dispersal and immigration strategies in primate males. Malayan Nat. J. 69, 353–356 (2017).
Google Scholar 

90.
Marty, P. R., Hodges, K., Heistermann, M., Agil, M. & Engelhardt, A. Is social dispersal stressful? A study in male crested macaques (Macaca nigra). Horm. Behav. 87, 62–68 (2017).
CAS  PubMed  Article  Google Scholar 

91.
Brent, L. J. N., Semple, S., Dubuc, C., Heistermann, M. & MacLarnon, A. Social capital and physiological stress levels in free-ranging adult female rhesus macaques. Physiol. Behav. 102, 76–83 (2011).
CAS  PubMed  Article  Google Scholar 

92.
Brent, L. J. N. Friends of friends: Are indirect connections in social networks important to animal behaviour?. Anim. Behav. 103, 211–222 (2015).
PubMed  PubMed Central  Article  Google Scholar 

93.
Blomberg, S. P., Garland, T. & Ives, A. R. Testing for phylogenetic signal in comparative data: Behavioral traits are more labile. Evolution (N. Y.) 57, 717–745 (2003).
Google Scholar 

94.
Garamszegi, L. Z. & Møller, A. P. Effects of sample size and intraspecific variation in phylogenetic comparative studies: A meta-analytic review. Biol. Rev. 85, 797–805 (2010).
PubMed  Google Scholar 

95.
Balasubramaniam, K. N. et al. The influence of phylogeny, social style, and sociodemographic factors on macaque social network structure. Am. J. Primatol. 80, 20 (2018).
Article  Google Scholar 

96.
Balasubramaniam, K. N. et al. Affiliation and disease risk: Social networks mediate gut microbial transmission among rhesus macaques. Anim. Behav. 151, 131–143 (2019).
PubMed  PubMed Central  Article  Google Scholar  More

Outline
In this section, first, we outline the method as a whole to give the reader an orientation and context for the following two parts. Second, we discuss the input data (for quantification and monetization). Third, we explain the merging of all input data, and thus the calculation of the output data. Finally, we address the influence of uncertainties on our method. The reference year for this analysis is 2016, and the reference country is Germany, which is listed as the third most affected country in the Global Climate Risk Index 2020 Ranking70.
Method in short
We differentiate between two steps within this method of calculating food-category-specific externalities and the resulting external costs. These are first the quantification and second the monetization of externalities from GHGs (visualized in Fig. 3). We use this bottom-up approach following the example of Grinsven et al.10, who conducted a cost-benefit analysis of reactive nitrogen emissions from the agricultural sector. This two-stepped method also allows the adequately differentiated assessment for GHG emissions of various food categories.
Fig. 3: Visualization of the method.

The method includes quantifying and monetizing product-specific externalities. In the case of Germany, emission data were obtained from the Global Emission Model for Integrated Systems (GEMIS)44. We used production data from the German Federal Statistical Office88 and AMI89,90, and calculated the emission difference between organic and conventional production based on a meta-analytical approach (see “Results” subsection on input data for quantification). The category-specific emission data were calculated on the basis of these input data. The emission cost rate was obtained from the German Federal Environmental Agency (UBA)32. The category-specific external costs were determined on the basis of the previously developed price-quantity-framework (see “Results” subsection on input data for monetization).

Full size image

The quantification includes the determination of food-specific GHG emissions—also known as carbon footprints39—occurring from cradle to farmgate by the usage of a material-flow analysis tool. Carbon footprints are understood within this paper in line with Pandey et al.71 where all climate-relevant gases, which (in addition to CO2) include methane (CH4) and nitrous oxide (N2O), are considered. Their 100-year CO2 equivalents conversion factors are henceforth defined as 28 and 265, respectively72. Here, the material-flow analysis tool GEMIS (Global Emission model for Integrated Systems)44 is used, which offers data for a variety of conventionally farmed foodstuff. As GEMIS data focus on emissions from conventional agricultural systems, we carried out the distinction to organic systems ourselves. We determined the difference in GHG emissions between the systems by applying meta-analytical methods to studies comparing the systems’ GHG emissions directly to one another. Meta-analysis is commonly used in the agricultural context, for example, when comparing the productivity of both systems57,58,59 or their performance1.
For better communicability, we first aggregate the 11 food-specific datasets given in GEMIS to the broader food categories plant-based, animal-based, and dairy by weighting them with their German production quantities (cf. “Results“ subsection on quantification). On top of that, LUC emissions are calculated for conventional foodstuff.
Through monetization, these emission data are translated into monetary values, which constitute the category-specific external costs. The ratio of external costs to the foodstuff’s producer price represents the percentage which would have to be added on top of the current food price to internalize externalities from GHGs and depict the true value of the examined foodstuff.
Input data for quantification
Starting with the data on food-specific emissions, GEMIS is used because of its large database of life-cycle data on agricultural products with a geographic focus on Germany. GEMIS is a World-Bank acknowledged tool for their platform on climate-smart planning and drew on 671 references, which are traced back to 13 different databases. The German Federal Environmental Agency uses GEMIS as a database for their projects and reports establishing it to be an adequate tool for the German context especially73,74. This tool is provided by the International Institute for Sustainability Analysis and Strategy (IINAS). GEMIS offers a complete view on the life cycle of a product, from primary energy and resource extraction to the construction and usage of facilities and transport systems. As GEMIS only offers data for the year 2010, we conducted a linear regression on the basis of the prevailing emission trend for the German agricultural context in order to align the data with the reference year 201675. For this, annual German emission data from 2000 to 2015 from the Federal Environmental Agency of Germany was used76. On every level of the process chain, data on energy- and material-input, as well as data on output of waste material and emissions, are provided by GEMIS. These data consist partly of self-compiled data from IINAS and partly of data from third-party academic research or other life-cycle assessment tools. Specific information on the data sources is available for every dataset of a product. In this study, the system boundaries for assessing food-specific GHG emissions span from cradle to farmgate. This means that we consider all resource inputs and outputs during production up to the point of selling by the primary producer (farmgate). This includes emissions from all production-relevant transports as well as emissions linked to the preliminary building of production-relevant infrastructure.
We specify that for animal-based products, emissions from feed production, as a necessary resource input, are assigned to these animal-based products. Such emissions naturally should include LUC emissions. LUC emissions are of negligible proportion for locally grown products, as agricultural land area is slightly decreasing in Germany47. Thus, we have to focus solely on imported feed for conventional animal-based and dairy products. Organic feed is not considered as article 14d of the EU-Eco regulation stipulates that organic farms have to primarily use feed which they produce themselves or which was produced from other organic farms in the same region77. The region is understood as the same or the directly neighboring federal state46. Although the EU-Eco regulation does not completely rule out fodder imports from foreign countries, it limits its application significantly. Also, one has to consider that over 60% of the organic agricultural area belongs to organic farming associations78. These associations stipulate even stricter rules than the standard EU-eco regulation. Examples are Bioland, where imports from other EU and third countries are only allowed as a time-limited exception79, Naturland, where additionally imports of soy are banned completely80, or Neuland, that ban any fodder imports from overseas81. We thus assume that the emissions that could possibly be caused by organic farming in Germany through the import of feed constitute a negligibly small fraction of the total emissions of a product. Thus, we follow common assumptions from the literature82,83,84 and calculate no LUC emissions for organic products. For conventional products, we calculate LUC emissions by application of the method of Ponsioen and Blonk45. This method allows the calculation of LUC emissions for a specific crop in a specific country for a specific year. With regards to the year, we apply our reference year 2016. With regards to crop and country one has to keep in mind that in the case of Germany, the net imports of feed are the highest for soymeal, followed by maize and rapeseed meal, making up over 90% of all net positive feed imports85. Maize and rapeseed meal are both imported mainly from Russia and Ukraine (93% and 87% of all imports86). Taken together, the crop area of Russia and Ukraine is decreasing by 150,000 ha/year (data from 1990 to 2015 were used87). Following Ponsioen and Blonk45, we thus assume that there are no LUC emissions of agricultural products from these countries. This leaves us with soymeal, of which 97% are imported from Argentina and Brazil. We thus calculate LUC emissions of soymeal for Argentina and Brazil, respectively. Data are used from Ponsioen and Blonk45, except for the data of the crop area, where updated data from FAOSTAT are used in order to match the reference year. We then weigh those country-specific emission values according to their import quantity. This results in 2.54 kg CO2eq/kg soymeal. To incorporate this value into the conventional emission data from GEMIS, we map the LUC emissions to all the soymeal inputs connected to the food-specific products.
For aggregation to narrow categories, we categorize every dataset from GEMIS into one of the eleven narrow food categories. The choice of separation into these specific categories is based on the categorization of the German Federal Office of Statistics88 from which production data were obtained. According to one category’s yearly production quantity, we incorporate every food product into the weighted mean of its corresponding food category. Thus, the higher a food’s production quantity, the greater the weight of this product’s emission data in the broad category’s emission mean. All data on the production quantities refer to food produced in Germany in the year 2016. For this weighting and aggregation step, only production quantities used for human nutrition were considered, thus feed and industry usage of food are ruled out (in contrast to emission calculation, where feed is indeed considered). Besides the German Federal Office of Statistics88, the source for this data is the German Society for Information on the Agricultural Market (AMI)89,90. Only production data for conventional production is used. Thereby, we imply ratios of production quantities across the food categories for organic production that are equal to those of conventional production. This does not fully reflect the current situation of organic production properties but allows for a fair comparison between the emission data of organic and conventional food categories. Doing otherwise would create ratios between emission values of organic and conventional broad categories that would not be representative of the ratios between organic and conventional narrow categories. In Table 4, all production data are listed, whereby total production quantities in 1000 t can be found in the right column. Translating these into percentage shares, the column right to the narrow category’s column represents the shares of the specific foods inside the narrow categories, whereas the column right to the broad category’s column represents the shares of the narrow categories inside the broad categories. These shares are expressed in formula 2a and 2b (see “Method and data” subsection on output data) by the terms (frac{{p_{b,n,conv}}}{{P_{b,conv}}}) (share in broad categories) and (frac{{q_{b,n,i,conv}}}{{p_{b,n,conv}}}) (share in narrow categories).
We aggregate GEMIS emission data (qb,n,i,conv) to narrow (eb,n,conv) and broad categories (Eb,conv) by multiplying the respective emission data with the shares from Table 3 (cf. formula 2a and b, “Method and data“ subsection on output data). From these conventional emission values, we derive emissions for organic production. For narrow as well as broad categories, the respective conventional emission values are multiplied with the applicable emission differences Db,org/conv (cf. Table 2).
With these data, we aggregate the above mentioned eleven food categories to three broad categories: plant-based, animal-based, and dairy. Besides the obvious differentiation between animal- and plant-based products, dairy is considered separately from other animal-based products because of its relatively high production volume and its, in contrast to that, relatively low externalities. Because the weighted mean of the three main categories is affected by the production quantities of its corresponding subcategories, mapping dairy into the animal-based category would otherwise distort the emission data of this very category.
As outlined before, only data regarding externalities of conventional agricultural production are included in GEMIS and could therefore be aggregated. Nevertheless, by applying meta-analytical methods regarding the percentage difference of GHG emissions between conventional and organic production, we derive the emission data for organic production for each of the broad categories (plant-based, animal-based, and dairy). It has to be noted that LUC emissions are consistently excluded at this level of calculation. To derive emission differences between organic and conventional farming, research was conducted by snowball sampling from already existing and thematically fitting meta-analysis, by keyword searching in research databases, as well as forward and backward search on the basis of already-known sources. Criteria for selected studies were climatic and regulative comparability to Germany. In the selected studies, relative externalities between conventional and organic farming are compared in relation to the cropland. To cover a reasonably relevant period, we decided to search for studies published within the past 50 years (from 1969 to 2018) and could therefore identify fifteen relevant studies, spanning from 1995 to 2015. Four of these studies have Germany as their reference country while the other eleven focus on other European countries (Denmark, France, Ireland, Netherlands, Spain, UK; please consult Table 2 for specifics). The weighted mean of the individual study results amounts to the difference in GHG emissions between the two farming production systems. As the selected studies are based on geophysical measurements and not on inferential statistics, a weighting based on the standard error of the primary study results like in standard meta-analysis91 was not possible. We aimed for a system that weights the underlying studies regarding their quality and therefore including their results weighted accordingly in our calculations. Within the scope of classic meta-analyses92, the studies’ individual quality is estimated according to their reported standard error (SE), which is understood as a measure of uncertainty: the smaller the SE, the higher the weight that is assigned to the regarding the source. Due to the varying estimation methods of considered studies, the majority of considered papers does not report measures of deviation for their results. These state definite values; therefore, there is no information about the precision of the results at hand. Against this background, we have decided to use a modified approach to estimate the considered papers’ qualities93. Following van Ewijk et al.94 and Haase et al.95, we apply three relevant context-sensitive variables to approximate the standard error of the dependent variable and thereby evaluate the quality of each publication: the newer the paper (compared to the timeframe between 1995 and 2018), the higher we assume the quality of reported results. The more often a paper was cited per year (measured on the basis of Google Scholar), the higher the paper’s reputation. The higher the publishing journal’s impact factor (measured with the SciMago journal ranking), the higher its reputation and therefore, the paper’s quality. For every paper, the three indicators publishing year (shortened with PY in Table 2), citations/year (CY), and journal rank (SJR) rank a paper’s impact on a scale from 1 to 10, where 1 describes the lowest qualitative rank and 10 the highest. The sum of these three factors (SUM) then determines the weight of a paper’s result in the mean value (WEIGHT). The papers’ reported emission differences between organic and conventional (diff. org/conv) are weighted with the papers’ specifically calculated WEIGHTS and finally aggregated to the emission difference between both systems.
With this approach, we weight results of qualitatively valuable papers higher and are therefore able to reduce the level of uncertainty in the estimated values because standard errors could—due to inconsistencies in the underlying studies—not be used. The results of this meta-analytical approach are listed in Table 2 (cf. “Results” subsection on quantification); further details can be found in Supplementary Note 1 and Supplementary Table 1. The studies considered compare GHG emissions of farming systems in relation to the crop/farm area. However, since our study aims to compare GHG emissions in relation to the weight of foodstuff, we include the difference in yield (yield gap) between the two farming systems for plant-based products and the difference in productivity (productivity gap) for animal-based and dairy products. For plant-based products, the yield gap is 117%, meaning that conventional farming produces 17% more plant-based products than organic farming in a given area. This gap was derived from three comprehensive meta studies57,58,59 and weighted as just described for the emission difference between organic and conventional farming. For animal-based as well as dairy products, the productivity gap could be determined with the same studies used for the meta-analytical estimation of the emission differences22,23,24,25,28,95. The productivity gap is 179% for animal-based and 152% for dairy products. In line with Sanders and Hess63, the yield (or productivity) difference (frac{{yield_{conv}}}{{yield_{org}}}) affects the calculation of the food-weight-specific emission difference (frac{{GHG_{org;food;weight}}}{{GHG_{conv;food;weight}}} = D_{org/conv}) between both farming systems: the yield difference is hereby multiplied with the cropland-specific emission difference (frac{{GHG_{org;cropland}}}{{GHG_{conv;crioland}}}). Resulting from this, the emission difference can be formulated as follows:

$$D_{org/conv} = frac{{GHG_{org;food;weight}}}{{GHG_{conv;food;weight}}} = frac{{GHG_{org;cropland}}}{{GHG_{conv;cropland}}} times frac{{yield_{conv}}}{{yield_{org}}}$$
(1)

If the yield difference were not included, emissions from organic farming would appear lower than they actually are as organic farming has lower emissions per kg of foodstuff but also lower yields per area. With formula 1, we adjust for that.
Input data for monetization
Monetization of these externalities requires data on GHG costs as well as data on the food categories’ producer prices.
The cost rate for CO2 equivalents used in this study stems from the guidelines of the German Federal Environment Agency (UBA) on estimating external ecological costs32. They recommend a cost rate of 180 € per ton of CO2 equivalents. This value is very close to the value of the 5th IPCC Assessment Report (173.5 €/tCO2eq), where the mean of all (up to this point) available studies with a time preference rate of 1% was determined33. The cost rate from the German Federal Environment Agency’s guideline is based on the cost damage model FUND96 and includes an equity weighting as well as a time preference rate of 1% for future damages. In this model, different impact categories are considered in order to estimate external costs from GHG emissions. Damage costs can be differentiated as benefit losses such as lowered life expectancy or agricultural yield losses and costs of damage reduction such as medical treatment costs or water purification costs97. Following UBA, these damage costs are analyzed in the following categories: agriculture, forestry, sea-level rise, cardiovascular and respiratory disorders related to cold and heat stress, malaria, dengue fever, schistosomiasis, diarrhea, energy consumption, water resources, and unmanaged ecosystems96. Using a cost-benefit-analysis (CBA), an adequate level of emissions is reached when marginal abatement costs are equal with damage costs. In a CBA external damage, costs can therefore be conceptualized as a price surcharge necessary to effect their optimal reduction98.
For the pricing of the food categories, we determine the total amount of proceeds that farmers accumulate for their sold foodstuff in €99 for each category (producer price) divided by its total production quantity. Thereby we calculate the relative price per ton for each foodstuff. We solely refer to producer prices as the system boundaries only reach until the farmgate.
Calculating output data
Output data include the aggregation and separation of food-specific categories to the broader categories of animal- and plant-based products, as well as conventional and organic products. As previously explained, such aggregation and separation are needed because the underlying material-flow analysis tool only lists food-specific emission data for conventionally produced foodstuff. Combining the input data, we are now able to quantify and monetize externalities of GHGs for different food categories.
For quantification, we separate between the following two steps: first, the aggregation of emissions data to broader categories and second the differentiation between conventional and organic farming systems. We iterate these steps two times, once for broad categories of animal-based products, plant-based products, and dairy and once for more narrow categories of vegetables, fruits, root crops, legumes, cereal, and oilseeds on the plant-based side as well as milk, eggs, poultry, ruminant, and pig on the animal-based side. Figure 4 displays the whole process of quantification schematically before we describe it in detail in the following text.
Fig. 4: Visualization of the quantification process.

Quantification as well as corresponding input and output data are displayed. Data from the Global Emissions Model for Integrated Systems (GEMIS)44 (gb,n,i,conv) and production data88,89,90 (qb,n,i,conv) are combined, and emission data for broad (Eb,conv) and narrow (eb,n,conv) categories are derived for conventional production. Organic emission values are calculated by multiplication of conventional emission values (Eb,org and eb,n,org) with the emission difference (Db,org/conv) (cf. “Input data for quantification”).

Full size image

Concerning the reasoning behind the method, the question that might come to mind is why the differentiation between farming systems happens after the aggregation and not before. This is due to the fact that the proportional production quantities of specific food as well as food categories to each other differ from conventional to organic production. Let us imagine aggregation would take place after the differentiation of farming systems: for example, beef actually makes up over 50% of all produced food in the organic animal-based product category, while it only accounts for 25% of the conventional animal-based product category (cf. production values in Table 3). As beef production produces the highest emissions of all foodstuffs, these high emissions would be weighted far stronger in the organic category than in the conventional category and thereby producing a higher mean for the organic animal-based product category than for the conventional one. As can be seen from this example, the organic animal-based product category could have a higher mean of emissions than the conventional animal-based product category while still having lower emissions for each individual organic animal-based product than conventional production. Deriving GHG emissions of foodstuff before aggregating to broader categories would thus be problematic and create means not representative for the elements that make up the broader category. To prevent this problem, the chosen method in this paper is thus to first aggregate to the chosen level of granularity (broad or narrow food categories) and then to derive emissions of organic production from conventional production data.
The first step of aggregation consists first of aggregating food-specific emission data from GEMIS gb,n,i,conv to the narrow categories eb,n,conv and second aggregating emission data from the narrow categories to the broad categories Eb,conv. As mentioned before and remarked in the respective indices, all these data only refer to conventional production up to this point. For both steps, the method is identical. The aggregation to narrow categories is represented in (2a) where eb,n,conv stands for the emissions of the narrow category n, which itself is part of the broad category b. Input data from GEMIS are remarked as gb,n,i,conv, whereby the index i refers to the ith element of category n. It’s production quantity is qb,n,i,conv. pb,n,conv represents the production quantity of the narrow category n. I (and N in formula 2b) represents the highest index of an element in a narrow (or a broad) category.

$$e_{b,n = x,conv} = mathop {sum }limits_{i in n = x}^I g_{b,n,i,conv} times frac{{q_{b,n,i,conv}}}{{p_{b,n,conv}}}$$
(2a)

The aggregation to broad categories is described by formula 2b whereby Eb,conv are the emissions and Pb,conv the production quantity of broad category b.

$$E_{b = x,conv} = mathop {sum }limits_{n in b = x}^N e_{b,n,conv} times frac{{p_{b,n,conv}}}{{P_{b,conv}}}.$$
(2b)

In the second step, we calculate emission values for organic production by multiplying the calculated emission difference Db,org/conv between both farming systems (cf. “Input data for quantification”) with the conventional emission values. These organic emission values are denoted as Eb,org for broad categories and eb,n,org for narrow categories.
To calculate the costs Cb of category-specific emissions, we multiply the cost rate P for CO2 equivalents with the category-specific emission data Eb or eb,n (depending on whether broad or narrow categories are observed). Further, we determine percentage surcharge costs (Delta _b) by setting these costs in relation to the producer price ppb of the respective food category: (Delta _b = frac{{C_b}}{{pp_b}}) (the calculation is analogue for narrow categories). These surcharge costs represent the price increase necessary to internalize all externalities from GHG emissions for a specific food category.
Dealing with uncertainties
Due to the interdisciplinarity and novelty of our study, we connect several methodological approaches and refer to various sources for data. Against this background, we had to accept some uncertainties while assembling and using the developed framework for our calculation. The studies included in our meta-analytical approach of calculating the difference between organic and conventional emission values, for one, are not fully consistent in the methodologies each of them uses (refer to Supplementary Table 1 for details). Furthermore, from the results of all included studies, it is apparent that there exists a wide range of emission differences between the farming practices, depending on the paper’s scope and examined produce21. We attempted to account for this by performing the studies according to their fit regarding the object of research (cf. “Input data for quantification”). Due to insufficient availability of the data for the emission differences between organic and conventional on the basis of each narrow category, an average for the emission difference was used. This possibly results in imprecisions during the internalization of the external costs on the level of all narrow categories. Therefore, we focus on the aggregated broad categories, as this uncertainty can be evaded here. Furthermore, the in literature reported price factor for CO2 equivalents is volatile over time, impacting the results of this paper. It is to be expected that the external costs of GHG emissions are likely to rise in the future (cf. subsection on research aim and literature review). Also, our study’s scope is confined to the assessment of the current production situation within the German agricultural sector. Therefore, we do not account for future developments regarding a changing agricultural production landscape after internalization of the accounted external costs. We do, however, discuss possible effects on demand patterns as well as the environmental and social performance of the agricultural sector in “Discussion”. Regarding the incorporated LUC emissions, there appears to be a lacking scientific consensus on a general method of calculation for such emissions45,100,101,102. We thus want to emphasize that these additional emissions should be treated with caution and are thereby displayed separately from the other data.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More