Introduction

The field of neuroscience has experienced remarkable advancements over the past decade, specifically propelled by the advent of innovative imaging techniques. The simultaneous application of positron emission tomography (PET) and magnetic resonance imaging (MRI) has revolutionized our ability to assess brain function across various physiological dimensions [1, 2]. This is especially true in the context of drug effect evaluations. Combining PET with functional MRI (fMRI) has opened many investigative avenues. Techniques such as resting-state functional connectivity (rs-FC) [3] or pharmacological MRI (phMRI) [4] are now complemented by the ability of PET to quantify molecular changes, including alterations in receptor and transporter availability [5].

While fMRI provides high spatial and temporal resolution, the interpretation of its readout necessitates caution. The blood-oxygen-level-dependent (BOLD) signal in fMRI indirectly captures neuronal changes through neurovascular coupling [6], revealing only the hemodynamic consequences of molecular-level drug effects. The integration of simultaneous PET acquisition can bridge this interpretative gap by offering essential molecular insights, particularly regarding transporter or receptor alterations. Typically, PET, used independently or simultaneously with fMRI, has been used mainly in pharmacological studies to illustrate quantitative shifts in neuroreceptor or transporter availability [7]. Several studies have also explored the interregional coherence of PET tracer signals [810], an approach akin to fMRI-derived rs-FC. While subject-level metabolic connectivity using [18F]FDG-PET has been established through the temporal correlation of regional PET signals [1, 8], studies employing transporter or receptor tracers have focused predominantly on interregional binding coherence across subjects using static scans [9, 10]. The concept and feasibility of molecular connectivity (MC) [11, 12] through the temporal correlation of the dynamic binding potentials of transporter or receptor tracers has yet to be explored.

In our study, we investigated the feasibility of deriving MCs from dynamic [11C]DASB-PET scans acquired simultaneously with fMRI in rats. We divided the rats into two cohorts, a baseline group and a pharmacological application group, which were exposed to 3,4-methylenedeoxymetamphetamine (MDMA). The baseline cohort assessed the feasibility of the novel methodology, contrasting it with traditional fMRI-derived rs-FC with a specific focus on temporal stability. We postulated that dynamic [11C]DASB PET temporal fluctuations could be harnessed for connectivity data like that used for hemodynamic rs-FC using seed-based and independent component analysis (ICA).

The second cohort, subjected to an MDMA challenge, allowed us to evaluate the utility of our novel approach. We aimed to outline the effects of MDMA by integrating the innovative MC concept with established analysis techniques. Our primary objective in this research was to elucidate the potential of PET-derived MCs in conjunction with simultaneous PET/fMRI, exploring the avenues this methodology could open for future diagnostic and drug development studies.

Materials and Methods

Our study reevaluates two datasets previously published by our group [13, 14] to explore FC and MC simultaneously at baseline [14] and following MDMA administration [13]. Please refer to these earlier publications for detailed descriptions of the animal handling methods, experimental setups, and data acquisition procedures.

Animals

A total of 41 male Lewis rats were obtained from Charles River Laboratories (Sulzfeld, Germany). Thirty rats (354 ± 37 g, corresponding to 11 weeks of age) underwent baseline [11C]DASB PET/fMRI scans without any pharmacological intervention, whereas 11 rats (365 ± 19 g, corresponding to 11 weeks of age) underwent [11C]DASB PET/fMRI scans, including an acute MDMA challenge. All experiments were in compliance with German federal regulations for experimental animals and received approval from the Regierungspräsidium Tübingen.

Simultaneous PET/fMRI Experiments

The rats were subjected to simultaneous PET/fMRI experiments involving 1.3% isoflurane anesthesia, tail vein catheterization, positioning on a temperature-controlled bed, and monitoring of vital signs. The scans were performed using a 7T small-animal ClinScan MRI scanner (Bruker Biospin, Ettlingen, Germany) with a custom-developed PET insert [15]. For [11C]DASB PET, we employed a bolus plus constant infusion protocol (kbol = 38.7 min) and reconstructed the scans in 1-minute time frames. The scanning protocol and sequence parameters are outlined in detail in our previous publication [13]. The MDMA cohort received a pharmacological MDMA challenge of 3.2 mg/kg intravenously 40 minutes after tracer injection.

Data Preprocessing and Analysis

Data preprocessing followed established protocols, including steps such as realignment, mask creation, coregistration, spatial normalization, signal cleaning, and spatial smoothing, as detailed in our previous work [13]. For the MDMA dataset, PET scans were analyzed for early and late effects post-challenge using the general linear model (GLM) available in SPM. For both datasets, the baseline was defined as 30 to 40 minutes after the start of the scan. For the fMRI data, a first-level analysis was applied to the individual scans without a high-pass filter (the filter was set to “Inf”). Statistical parametric maps were generated after GLM parameter estimation using contrast vectors. The group-level analysis involved a one-sample t-test on the subject-level statistical parametric maps (p < 0.05, one-sided, familywise error/FWE adjusted).

Static PET scans were generated by summing dynamic frames over defined periods for 10-minute periods after the MDMA challenge (50-60 minutes to investigate early effects and 70-80 minutes to investigate late MDMA effects). Two-sample t-maps were calculated between the normalized [11C]DASB uptakes of (1) the baseline scan period and the early effect period and (2) the early effect period and the late effect period (p < 0.05, FWE-adjusted).

All group-level t-maps were subjected to voxelwise signal quantification to determine the regional contributions of 48 regions selected according to the Schiffer atlas [16]. The average t-scores and standard deviations of all voxels were calculated.

Functional Connectivity Analysis

FC was determined using a seed-based analysis approach. The mean time series of the preprocessed BOLD-fMRI signals for each dataset across all regions (refer to Supplementary Table 1 for the list of regions) were extracted using the SPM toolbox Marseille Boîte À Région d’Intérêt (MarsBaR). Pairwise Pearson’s correlation coefficients were calculated between each pair of mean regional time series for every dataset. The Pearson’s r coefficients were converted into z values using Fischer’s transformation for group-level analysis. Fischer’s z-transformed correlation coefficients were then used to generate mean correlation matrices for both cohorts [17].

Molecular connectivity analysis

The mean [11C]DASB signal from the preprocessed PET datasets was extracted from the designed regions, including the 48 regions used for fMRI data analysis and the cerebellum using MarsBaR. Binding potentials were calculated framewise for all dynamic PET scans using the DVR-1 (Equation 1) to generate regional BPND values with the cerebellar gray matter as a reference region, which our earlier studies demonstrated to be the most appropriate for this tracer in rats [18, 19]:

where:

  • BPND is the binding potential

  • VT is the total volume of distribution

  • VND is the volume of distribution in a reference tissue

  • DVR is the relative volume of distribution

To calculate the MC, we discarded the first 20 minutes of every scan, which were dominated by perfusion effects, and applied a detrending approach to the remaining 60 minutes to obtain temporally stable values (for further details, please refer to the Supplementary Methods and Supplementary Figure 1). The BPND time courses were then used to calculate the MC as described above for fMRI: ROI-to-ROI subject-level correlation matrices between all regional time courses were generated, and z-transformed correlation coefficients were used to calculate mean correlation matrices.

Independent Component Analysis

Group ICA (GIFT toolbox, MIALAB, University of New Mexico, Albuquerque, NM, USA) was used for ICA in the baseline group. The fMRI and PET preprocessed datasets were investigated between 30 and 80 minutes after the start of data acquisition. For fMRI, we selected ten independent components, while we started with two components for PET and increased the number to ten components to thoroughly dissect the varying components within the signal. The components were thresholded at a z-value ≥ 1.96 (p-value ≤ 0.05) [20], and average z-scores and standard deviations were calculated for each component. The physiological significance of the components was further explored by contrasting them with the regional [11C]DASB changes induced by MDMA. Accordingly, the z scores of the independent components generated in the baseline cohort were correlated with the early and late regional [11C]DASB changes induced by MDMA, measured using t scores.

Results

Comparability of MC and FC in spatial contexts and over time

We first aimed to evaluate whether MCs align spatially with FC, possess similar graph theory properties and provide consistent temporal readouts throughout the scan duration in the baseline group (Figure 1).

Evaluation of the seed-based MC.

(A) Correlation matrix indicating whole-brain FC (beneath the diagonal) and MC (above the diagonal). Correlations not significant with multiple comparison corrections were set to zero (p < 0.05, FWE correction). (B) Scatter plot and correlation between MC and FC edges. (C) Small-world coefficients for all subjects and group-level one-sample t-test against the value of 1 (SW > 1 indicates small-world properties). Comparison of (D) FC and (G) MC early (20-40 minutes after the start of the scan, below the diagonal) and late (60-80 minutes after the start of the scan, above the diagonal). The similarities of early and late readouts were quantified for both (E) FC and (H) MC. The temporal stability of both (F) FC and (I) MC was evaluated using a sliding window approach, including 20-minute windows between 20 and 80 minutes after the start of the scan. Abbreviations: FC = fMRI-derived hemodynamic functional connectivity, MC = [11C]DASB PET-derived molecular connectivity.

We found a moderate but significant correlation between the edge-level MC and FC (r = 0.51, p < 0.001; Figure 1A and B). Furthermore, both connectomes revealed small-world properties at the group level, with coefficients higher than 1 (Figure 1C). At the subject level, three molecular and four functional connectomes fell below the threshold of 1 for the small-world coefficient. Significant consistency was observed between early and late scan-derived connectomes (Figure 1D for FC, G for MC), with FC having a slight edge (Figure 1E, r = 0.96) over MC (Figure 1H, r = 0.8). While both FC and MC maintained steady correlation intensities, there was a negligible decline over the scan duration (Figure 1F and 1I).

Deciphering the spatial characteristics of FC and MC using ICA

After establishing the feasibility of obtaining temporally stable readouts using the ROI-to-ROI approach, we employed a data-driven approach using ICA to compare the spatial characteristics of FC and MC (Figure 2).

Group independent component analysis for FC and MC.

(A) ICA was performed over 10 components for fMRI. (B) ICA performed over 2 components for [11C]DASB PET and regional quantification of the two derived components. The ICA was repeated over 10 components. Four and three components showed good overlap with the two components defined above. All the components were thresholded at z > 1.96 (p ≤ 0.05). Abbreviations: FC = fMRI-derived hemodynamic functional connectivity, MC = [11C]DASB PET-derived molecular connectivity.

We extracted ten group ICs from the fMRI data (Figure 2A), revealing known canonical resting-state networks, such as the posterior default-mode-like network (IC1-5, red), sensorimotor networks (IC1-5, green and purple), the anterior default-mode-like network (IC6-10, yellow) and the visual network (IC6-10, red). Given the unpredictability of the number of ICs suitable for MC IC extraction, we started with two components (Figure 2B). IC1 (orange) comprised both subcortical and cortical anterior brain regions, including the nucleus accumbens, amygdala, cingulate cortex, caudate putamen, orbitofrontal cortex and medial prefrontal cortex, whereas IC2 (blue) primarily received contributions from deeper posterior areas, such as the midbrain, thalamus, hypothalamus, periaqueductal gray cortex and medulla. Interestingly, when we extracted 10 independent components to mimic the number of components used for the FC data, the initial anterior component split into four different ICs, and the initial posterior IC split into three different ICs. Relatively clear spatial segregation can be observed for the newly formed ICs, for instance, the three posterior components extracted from specific regions (green, mainly from the medulla; red, from the hypothalamus; and part of the midbrain, blue, from the midbrain and thalamus).

MDMA-induced changes in ICA-derived molecular connectivity

Next, we aimed to explore the relationship between molecular changes in SERT availability and the molecular connectome derived from the ICA induced by acute MDMA administration (Figure 3).

Comparison of MDMA-induced [11C]DASB alterations.

(A) Left panel: Dynamic binding potentials of regions comprising the SERT subcortical network, defined by IC 1 in the validation cohort. Right panel: Dynamic binding potentials of regions comprising the SERT salience network, defined by IC 2 in the validation cohort (continuous lines indicate the means, and dotted lines indicate standard deviations). (B) Overlap between independent components extracted from the validation cohort (IC 1 = SERT subcortical network, IC 2 = SERT salience network;) and the early and late effects of MDMA. (C) Pairwise correlations between regional z scores of the ICs extracted from the validation cohort and regional t scores of early and late MDMA effects. (*** indicates p < 0.001, ns = nonsignificant). Abbreviations: SERT = serotonin transporter, ICA = independent component analysis; for abbreviations of regions, please refer to the Supplementary Information).

Two ICs extracted from the MC showed good overlap with regions associated with the SERT salience network (IC1) and the SERT subcortical network, comprising regions with intrinsically high SERT densities (IC2). Therefore, we defined the regions contributing strongly to IC2 as the SERT salience network (CPu, Cg, NAc, Amyg, Ins, mPFC) and those with strong signals in IC1 as the SERT subcortical network (VTA, Th, MB, PAG, Hypo). Interestingly, MDMA induced immediate strong decreases in all SERT subcortical network regions, and salience areas exhibited a delay of approximately 10 minutes (Figure 3A). Voxel-level analysis revealed clear spatial overlaps between early MDMA-responsive regions and those from the posterior IC, with delayed regions mirroring the anterior IC reminiscent of the SERT salience network (Figure 3B). To quantify the striking spatial similarity between the baseline independent components and the spatiotemporal characteristics of [11C]DASB changes after MDMA exposure, we found highly significant correlations between the z scores of the posterior and anterior ICs and the t scores of late and early MDMA effects on [11C]DASB alterations (p < 0.0001, Figure 3C).

Finally, we investigated the relationships between SERT availability changes and MC reductions following acute MDMA challenge (Figure 4).

Comparison of MC and BPND changes following MDMA.

A Reductions in BPND following MDMA (orange) and MC strength (blue) were compared. B The correlation coefficient of the regional T scores was low (r = 0.29).

We found that the decreases in MC and [11C]DASB BPND following MDMA application (Figure 4A) showed only a low correlation (r = 0.29, Figure 4B). While decreases in SERT availability exhibited a strong anterior-posterior gradient, which was most pronounced in areas with high SERT availability, such as the MB, VTA, Pons or PAG, the MC encompassed regions across the brain to similar extents. Specifically, while the significance of [11C]DASB reductions in Ins was very low across the investigated regions, Ins presented the greatest effects among regions for MC. In contrast, strong SERT occupancy effects in IC and SC did not translate into very prominent reductions in the respective global MCs of the two regions.

MDMA-induced changes in seed-based molecular connectivity

Next, we performed a seed-based analysis to compare changes in FC and MC after an MDMA pharmacological challenge on the salience network and regions with high SERT binding (Figure 5).

MDMA effects on seed-based FC and MC of the salience and subcortical networks.

(A) FC and (C) MC brain networks depicting the edge and node strengths of the salience network and subcortical network at baseline (20-40 min after the start of the scan) and after MDMA (60-80 min after the start of the scan). (B) FC and (D) MC time-resolved salience and subcortical network strengths computed by sliding windows. Individual values are provided in the bar graphs for the baseline and final post-MDMA periods. Asterisks indicate significant (p < 0.05, FDR-corrected) changes from baseline (time point zero, corresponding to 20–40 minutes after the start of the scan).

We observed a decrease in SN FC (16% at the end of the scan, p = 0.02, did not survive FDR correction) and almost constant FC of the subcortical network (<5% decrease, p = 0.79) following MDMA, as shown in Figure 5A at the edge and node levels and Figure 5C at the network level. For MCs, we observed profound reductions in the subcortical network (Figure 5C and D, p = 0.009), emphasizing the acute and spatially specific effect of MDMA on MCs.

Discussion

The mammalian brain operates on diverse physiological, spatial and temporal scales. FC via BOLD-fMRI offers insights into coherent functional brain networks, but its complexity and indirect link to neural activity highlight the need for more direct methodologies. In this context, the concept of MC using [11C]DASB PET, as introduced in our study, provides a more direct and complementary perspective on brain organization and its response to external stimuli, such as MDMA.

Physiological basis

Our findings suggest that [11C]DASB binding reflects the interplay of serotonin levels and SERT dynamics [21]. In support of the competition model, evidence indicates that endogenous serotonin competes with tracers for binding sites, affecting tracer binding [22]. However, contrasting results from various studies highlight the complexity of this interaction [23] [24] [25]. The internalization model suggests that serotonin levels influence SERT internalization, impacting [11C]DASB binding [21, 26, 27]. While supporting evidence, further exploration is needed to fully understand these dynamics, especially during the resting state [28]. Some models on this topic have proposed a regulatory function of the raphé nuclei in maintaining serotonin fluctuations over several temporal scales at rest [29]. Remarkably, fast microdialysis has resolved multiple spontaneous surges of up to 1500-fold of the basal serotonin occurring during 30-minute intervals [30]. Additionally, the same study indicated that SERT expression is essential for spontaneous surges and that reduced SERT drastically decreases serotonin spiking. Thus, it is feasible that the correlated temporal fluctuations captured by dynamic [11C]DASB at least partly reflect the role of the serotonergic system in the resting activity of the brain. In addition, SERT regulation occurs over multiple time scales, ranging from milliseconds to hours, depending on the mechanism involved [31]. Rapid changes in SERT surface expression can be mediated by protein-protein interactions or posttranslational modifications [32, 33], such as phosphorylation, which occur on a timescale of milliseconds to minutes. These processes dynamically modulate surface availability and function, allowing fine-tuned regulation of serotonin uptake even under resting-state conditions. Additionally, while slower processes involving endocytosis, recycling, and degradation typically occur over minutes to hours, subtle fluctuations in SERT trafficking and activity can still occur under basal conditions. These minor yet biologically relevant changes likely reflect ongoing homeostatic regulation essential for maintaining serotonergic balance. Therefore, tracer fluctuations observed during resting-state measurements should not be dismissed, as they may represent meaningful variations in SERT regulation that contribute to the fine control of serotonin clearance.

Implications for the whole-brain serotonergic system

Our graph theory analysis revealed comparable whole-brain organization between the hemodynamic (BOLD-derived) and serotonergic (PET-derived) connectomes, although we found only a moderate correlation between the edges of the two measurements. Notably, the [11C]DASB-based MC results showed a distinct binding pattern compared with BOLD-fMRI-derived FC. For FC, the ICA revealed classic RSNs, including the default-mode network, sensory network and motor network, consistent with previous studies in rats [34] and humans [35]. In contrast, the PET data revealed two distinct anatomical components in the serotonergic system. One component included key subcortical areas, such as the brainstem, parts of the midbrain and the thalamus, regions rich in SERT availability and central to serotonergic regulation. The other component comprises regions of the limbic system, such as the striatum, amygdala, insula, cingulate cortex and prefrontal cortex, critical targets of serotonergic projections from the raphe nuclei. These components suggest a broader spatial reach of the serotonergic system than the typically described RSN, reflecting its unique neurochemical architecture. In our analysis, PET ICs appeared less bilateral than fMRI ICs. This is likely due to the lower temporal resolution of PET (100 frames) than of fMRI (3000 frames), resulting in a reduced signal-to-noise ratio (SNR) and potentially affecting the stability and symmetry of the independent components.

Our results demonstrate that compared with FC, MDMA induces more pronounced changes in MCs, particularly in regions associated with the SERT subcortical network. The distinct temporal dynamics of BPnd variations between these components may reflect the hierarchical organization of the serotonergic system. Specifically, the raphe nuclei, as the primary source of serotonin, are likely to exert more immediate modulation on posterior subcortical structures (IC2), whereas downstream effects on limbic and cortical regions (IC1) may occur more gradually. While these findings align with current neuroanatomical and molecular knowledge, the precise biological mechanisms driving these temporal differences remain unclear. Future investigations are warranted to elucidate these mechanisms. Future studies combining direct measurements of serotonin levels with neuroimaging data will be critical to fully understanding these components’ distinct roles and temporal profiles in regulating serotonergic function.

Relevance to prior research

Our findings align with Salvan et al. [36], who integrated molecular maps into fMRI data and demonstrated how individual serotonergic receptors contribute to network-level activity. Despite our differing methodologies, the receptor activity patterns found, may also correspond to the independent components we found from the [11C]DASB PET data. The dual regression approach reported by Salvan et al. was first reported in the context of pharmacology to delineate the receptor-specific effects of MDMA in humans [37]. The authors found that MDMA specifically decreased FC in the 5HT1A maps in areas that could be ascribed to the human salience network, such as the insula and a collection of medial cortical regions. While not reaching significance, we observed a reduced salience FC and MC in our data.

Previous studies evaluating the effects of MDMA on FC have indicated relatively few effects on FC, as confirmed by the results of this study [38, 39]. Recently, a novel approach [37] revealed decreased connectivity induced by MDMA in areas associated with SERT and 5-HT1A availability. The increased activity in limbic and cortical structures when controlling for vascular effects is accompanied by decreased salience connectivity, indicating that, while neurons become more active through the drug, they do so in an incoherent manner, which would be in line with the hyperactive yet abnormal behavior reported for MDMA abuse. Importantly, the potent vascular effects that play a role in modulating the amplitude of the BOLD-fMRI signal may also influence FC readouts, although the extent of such effects is difficult to estimate [13].

Over the past decade, PET studies using [11C]DASB have focused on the associations of serotonin transporter (SERT) availability across different brain regions, revealing altered interregional SERT connections in patient cohorts, posttreatment changes, and a predictive capacity for treatment response [9, 10, 40] [41] [42]. Our study builds upon this foundation, revealing significant within-subject temporal associations in [11C]DASB binding and demonstrates more pronounced alterations in MC than in FC induced by drugs, thus highlighting the enhanced utility of our method.

Study limitations and future directions

One limitation of our study is that our experimental protocols predate the recently published consensus recommendations for rat fMRI [43], particularly concerning anesthesia and preprocessing pipelines. Using isoflurane anesthesia, although common at the time of data acquisition, introduces a potential confound due to its known neuronal activity effects. However, we previously demonstrated that isoflurane at 1.3% maintains stable physiological parameters and avoids burst suppression [44], a concern at higher doses. Furthermore, other studies have reported that low-dose isoflurane remains feasible for resting-state functional connectivity studies [45]. While isoflurane, a GABA-A agonist, could theoretically interact with the mechanisms of MDMA in the brain, we found no evidence in the literature suggesting significant cross-talk between these substances. Future studies employing medetomidine-based protocols may help minimize this potential confounding factor.

With respect to data preprocessing, we retained the same pipeline used in our prior publications [13, 14] to maintain methodological consistency. While we recognize the advantages of adopting standardized preprocessing, as outlined in the consensus guidelines, this approach ensures comparability with our previous analyses. To facilitate further investigation, we have made the full dataset publicly available (see Data Availability Statement), enabling reanalysis with updated pipelines or additional explorations of this dataset.

Nonetheless, we provide a strong rationale for considering such analyses in future studies. First, we show that at rest, the data are reliable, temporally stable and exhibit similar graph theory metrics to traditionally calculated functional connectomes. Second, despite comparable network-level organizational properties, at rest, MC is only moderately correlated with FC, indicating the complementary nature of both readouts. Third, resting-state MC ICs correlate well with SERT occupancy changes induced by the MDMA challenge. Finally, we show that the changes that MDMA elicits on the MC are complementary to standardly calculated [11C]DASB BPND alterations. Our data indicate that while changes in BPND are more pronounced in regions with higher baseline SERT availabilities, MC reveals a more globally distributed measure of tracking serotonergic changes since regions, such as the insula, with relatively low SERT expression, are strongly affected. Additionally, our multimodal imaging approach, although powerful, cannot fully decipher the mechanisms of interregional coherence in PET time courses. Furthermore, studies employing more direct sensitive methods to measure neurotransmitter release could provide deeper insights into the molecular processes underpinning our observations.

Conclusion

This study significantly contributes to integrating molecular data into connectomic frameworks, demonstrating that subject-level MCs are reliable and complementary to FC in both resting and pharmacologically challenged states. Our research provides a strong foundation for future investigations into the value and generalizability of PET-derived MCs, particularly for understanding drug-induced brain-wide molecular network changes.

Additional information

Funding

This research was supported by funds from the Eberhard Karls University Tübingen Faculty of Medicine (fortüne 2209-0-0 and 2409-0-0) to HFW, from Bundesministerium für Bildung und Forschung (BMBF, Grant No. 01GQ1415) to HFW and from the Carl Zeiss Foundation to KH.

CRediT authorship contribution statement

Tudor M. Ionescu: Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing – original draft, Visualization. Mario Amend: Investigation, Writing – review & editing. Rakibul Hafiz: Methodology, Software, Writing – review & editing. Bharat B. Biswal: Conceptualization, Writing – review & editing, Supervision, Project administration, Funding acquisition. Andreas Maurer: Writing – review & editing. Hans F. Wehrl: Conceptualization, Methodology, Writing – review & editing, Supervision, Project administration, Funding acquisition. Kristina Herfert: Writing – review & editing, Visualization, Supervision, Project administration, Funding.

Data availability

The data are openly available at the DRYAD repository (https://doi.org/10.5061/dryad.6djh9w1bf).

Additional files

Supplemental methods