Ecology

Data collectionPeer-reviewed articles published up to Dec 31, 2021 were searched through the Web of Science (http://webofknowledge.com), Google scholar (http://scholar.google.com), and the China National Knowledge Infrastructure (CNKI, http://www.cnki.net). Here we employed a combination of the following search terms: “(tree biomass OR aboveground biomass OR plant biomass OR plant productivity) and (allometric biomass equation OR allometric model OR productivity model OR biomass equation OR biomass model)”. To avoid potential selection bias and duplicates, we conducted a cross-check between the references of relevant articles, which resulted in the selection of 729 relevant articles from the thousands of the appearing articles initially. Subsequently, eligible articles were selected using the following criteria: (1) Allometric models built for specific species with confirmed locations without disturbances were selected, generalized species, large-scales (e.g., province or nation), as well as recently disturbed tree models were excluded. (2) The method employed to develop the model was destructive harvesting and weighing, with at least twenty sample trees, were selected; articles were excluded that did not include measurements and used less than twenty sample trees. (3) The model forms were W = a*Db and LnW = a + b*Ln(D) or W = a*(D2H)b and LnW = a + b*Ln(D2H), where W is the aboveground biomass, and D is the diameter at breast height, H is the tree height, were selected. Consequently, we excluded articles with other variables and other forms of models. Finally, 426 articles remained from the original 729 (Supplementary Fig. S1).We then distilled data from the articles for the following variables: (1) Allometric models, in the form of W = a*Db and LnW = a + b*Ln(D), W = a*(D2H)b and LnW = a + b*Ln(D2H) including the parameters a, b in the D range and H range. (2) Tree species corresponding to the models, including families, genera, and species. (3) Location data, including longitude, latitude, and study sites. (4) Climate data, including mean annual temperature (MAT, °C) and mean annual precipitation (MAP, mm) of the tree species location. (5) Terrain data, including slope and aspect. (6) Soil data, including soil organic carbon (SOC), clay, and soil type.Since not all articles provided the location, climate, soil, and terrain data of the studies, we estimated the missing data as follows, (1) we supplemented the longitude and latitude with the study location using Google Earth. (2) We extracted the missing climate data by using geographic coordinates from WorldClim version 2.0 (http://worldclim.org/current)16. (3) We obtained the shuttle radar topographic mission DEM data with 30 m resolution from NASA, and used SAGA-GIS software to derive various terrain data from the DEM such as altitude, slope, and aspect17, 18. (4) The missing soil data was derived from the Regridded Harmonized World Soil Database v1.219. In particular, we established the soil type according to Soil Taxonomy to increase the accuracy of the analysis and prediction. Furthermore, if the experiments were performed at multiple sites in one study, they were treated as independent observations. In light of above criteria, 817 allometric models in the form of W = a*Db or LnW = a + b*Ln(D) and 612 allometric models in the form of W = a*(D2H)b or LnW = a + b*Ln(D2H) were collected from the 426 articles.Allometric modelThe relationship between the diameter and aboveground biomass was in the form of the power function20:$$begin{array}{c}Wi=atimes D{i}^{b},end{array}$$
(1)
where Wi is the dry mass of the ith tree (kg), Di is diameter at breast height (cm), and a and b are the parameters of the model.$$Wi=atimes (D{i}^{2}Hi{)}^{b},$$
(2)
where Wi is the dry mass of the ith tree (kg), Di is diameter at breast height (cm), Hi is the tree height (cm), and a and b are the parameters of the model.However, a heteroscedasticity exists when directly fitting the tree biomass. The logarithmic transformation of Eq. (1) or Eq. (2), is convenient to facilitate model fitting and deal with heterocedasticity21. The logarithmic transformation allometric model:$$begin{array}{c}Lnleft(Wiright)=a+btimes Lnleft(Diright),end{array}$$
(3)

was used in this function, where a (Eq. 3) represents Ln(a) (Eq. 1), and b (Eq. 3) is the same as b (Eq. 1), respectively.$$begin{array}{c}Lnleft(Wiright)=a+btimes Lnleft(D{i}^{2}Hright),end{array}$$
(4)
was used in this function, where a (Eq. 4) represents Ln(a) (Eq. 2), and b (Eq. 4) is the same as b (Eq. 2), respectively. To unify the models, we transformed the collected Eqs. (1) to (3) and Eqs. (2) to (4).Data analysisTo establish the relationship between variables with parameters a and b for making a parameter prediction on a global scale, Random Forest (RF) (an example of a machine learning model) was employed, which consists of an ensemble of randomized classification and regression trees (CART)21. In short, the RF will generate a number of trees and aggregate these to provide a single prediction. In regression problems the prediction is the average of the individual tree outputs, whereas in classification the trees vote by majority on the correct classification22, 23. Generated trees called ntree are based on a bootstrapped 2/3 sample of the original data to decrease correlations by choosing different training sets in the RF modeling process15. In addition to this normal bagging function, the best split at each node of the tree was searched only among a randomly selected subset (mtry) of predictors24. The tree growing procedure is performed recursively until the size of the node reaches a minimum, k, which is parameterized by the user. For the rest of the original data, RF provides a believable error estimation using the data called Out-Of-Bag (OOB), which is employed to obtain a running unbiased estimate of the classification error as trees are added to the forest15.Predictive variable selectionThe variables included stand factors such as density, family, and diameters, as well as non-stand factors such as MAT, MAP, and SOC. Considering that the prediction was on a global scale, the first step was to exclude the factors that it was not possible to completely extract. Next, we selected variables through the following22: (1) the RF classifier was initially applied using all of the predictor variables, and variable importance was used to rank them based on the mean decrease in accuracy. (2) Removing the least important variables by the variable importance ranking, (3) the training data were then partitioned five-fold for cross-validation and the error rates for each of the five cross-validation partitions were aggregated into a mean error rate, and 20 replicates of the five-fold CV were performed25.By means of the above, eleven variables, including family, genus, species, MAT, MAP, altitude, aspect, SOC, slope, clay, and soil type, were remained to predict parameters. Since the combinations of variables were different, five combinations were performed to make predictions from the eleven variables above. Among the five combinations, each were used by RF to predict and select via the model evaluation index VaR explained and the mean of squared residual (Supplementary Table S1).Optimization of Random Forest parametersRF depends primarily on three parameters that are set by users. (1) ntree, the number of trees in the forest. (2) nodesize, the minimum number of data points in each terminal node. (3) mtry, the number of features tried at each node. To obtain the optimization of RF parameters, we set ntree = 1000, 2000, 3000 and the selection criterion was that ntree was small enough to maximize computational efficiency as well as produced stable OOB error25. As for nodesize, we used 3, 5, 7, and 5 as the default for regression RF, given that the mtry value always is always one third of the number of variables. Here we also set the mtry values (ranging from 2 to 4), which were tested, and we accessed the OOB error rates from 50 replicates for each mtry value25. The primary tuning parameter above were optimized, as well as each combination of the three RF parameters through a grided search, which were used to predict and set RF parameters according to the predictive effect of each combination (Supplementary Table S2).All above data analysis were conducted in R 4.0.326. And the output is the spatial pattern of allometric model parameters at 0.5° resolution.Predicted parameter validationFurther to assess the accuracy of the predicted parameters, we applied them to estimate the AGB at six sites. And the actual AGB of the sites had been obtained via destructive sampling from 209 plots, which were located in Hubei, Liaoning, Gansu, Hebei and Heilongjiang provinces, and Inner Mongolia autonomous region from 2009 to 201327 (Table 1). First, we selected the sample trees according the dominant, average and inferior tree outside the plot. Then the sample trees were felled as carefully as possible and tree height (H), tree diameter in the breast (DBH) and live crown length were recorded. To divide trees into several sub-samples, including branches, leaves, stem wood and stem bark, all of the branches were removed and leaves were picked. Besides, stem was divided into 1 m sections and bark of the stem was removed. Finally, all sub-samples of aboveground part of trees were oven-dried at 80 °C until a constant weight was reached and the sum of all the sun-samples weight was the actual AGB. Through the above process, 249 actual AGB data were obtained. Meanwhile, the predicted parameters of the models together with the DBH and H estimated the predicted AGB. The actual AGB data of 249 sample trees were compared with the predicted AGB by making fitting curves between them in R to show the availability of predicted parameters according root mean square error (RMSE) and R2.Table 1 The basic features of the sampling sites.Full size tableThe experimental research and field studies on plants in this study, including the collection of plant material, complied with the relevant institutional, national, and international guidelines and legislation. And we ensured that we have permission for the plant sampling, all of the steps were allowed in our study for the plant research. In addition, plant identification in this study was conducted by X.Z according to World Plants (https://www.worldplants.de) in the herbarium of School of Forestry & Landscape of Architecture, Anhui Agricultural University, and the voucher specimen of all plant material has been deposited in a publicly available herbarium. More

The research in this study was performed in agreement with permission M45-14 issued by Malmö/Lund Ethical Committee for Animal Research, Sweden, which granted capture and blood sampling of wild birdsSamplesNine willow warblers, determined to be males (based on a wing length > 69 mm), were caught opportunistically with mist nets during the time of autumn migration in September 2016 at Krankesjön, 15 km East of Lund, Southern Sweden. While most of the individuals were phenotypically similar to willow warblers breeding in Southern Scandinavia, some were slightly larger and had a greyer plumage, which is more commonly seen in Northern Scandinavia12. The set of samples thus potentially contained willow warblers of each of the two major migratory phenotypes. Blood from each bird was collected through a puncture of the brachial vein and was stored in two vails containing SET buffer and 70% ethanol, respectively. An aliquot of the blood was used for DNA extraction with a phenol-chloroform protocol. From the extracted DNA, we genotyped the samples for two loci located on chromosomes 1 and 5, respectively (NBEA and FADS2)45,46, and for a bi-allelic marker within the divergent region on chromosome 3 (AFLP-ww1)47. Based on the genotyping results we selected two samples that were homozygous northern or homozygous southern for all three loci, respectively. We also included a sample from a chiffchaff Phylloscopus collybita (female) for de novo genome sequencing of a closely related outgroup species, as well as an additional willow warbler (DD81063, male) to confirm breakpoint differences with linked read sequencing. Both of these birds were opportunistically caught at the same site as above during autumn migration in 2019, and collection of blood followed the same approach as for the other birds.Optical mapsDNA from the northern and southern willow warbler was extracted from blood stored in ethanol using a Plug Lysis protocol (v.30026D; Bionano Genomics, CA, USA). The blood was first separated from the ethanol through gentle centrifugation and embedded in molten 2% agarose plugs (DNA plug kit; Bio-Rad, CA, USA). The solidified plugs were submerged in Lysis Buffer solution (Bionano Genomics) and 66.8 µl per ml Buffer Puregene Proteinase K (Qiagen,MD, USA) for 2 h at 50 °C. The plugs were subsequently washed in 1× Wash buffer (Bio-Rad DNA plug kit) followed by TE buffer. In the following step, the plugs were treated with RNase (Qiagen, 20 µl in 1 ml TE buffer) for 1 h at 37 °C, followed by another washing step using the same buffers as in the previous step. Next, the plugs were melted for 2 min at 70 °C and treated with GELase (Epicenter, WI, USA) for 45 min at 43 °C. The DNA was then purified from digested agarose using drop dialysis against TE buffer on a 0.1 µm dialysis membrane (MF-Millipore, Merck KGaA, Germany) for 2.5 h.Optical maps for each of the two samples were produced using Bionano Genomic’s commercial Irys system48. BspQ1 was determined to be the most suitable nicking enzyme after using the software LabelDensityCalculator v.1.3.0 and Knickers v.1.5.5 to analyze a previous short-read assembly13. Bionano Genomic’s IrysPrep Labeling-NLRS protocol (v.30024) was used for the NLRS reaction. For this step, DNA was treated with Nt.BspQ1 (NEB, MA, USA) to create single-stranded nicks in a molecule-specific pattern. These were then labeled with Bionano Genomic’s (CA, USA) labeling mix (NLRS kit), aided by Taq Polymerase (NEB), and repaired using Bionano Genomics’s repair mix (NLRS kit), in the presence of Thermopol Rxn buffer, NAD+, and Taq DNA Ligase (NEB). Finally, the DNA backbone was stained using DNA stain from Bionano Genomics’s NLRS kit. Each sample was then loaded on two IrysChips (Bionano Genomics) each, and the DNA with stained BspQ1 nicks was visualized using an Irys instrument, following Bionano Genomics’s Irys user guide (v.30047). This resulted in 200 and 182 Gb of data for the northern and southern sample, respectively.Genome maps were assembled de novo using Bionano Genomic’s in house software IrysView v.2.5.1, with noise parameter set to “autonoise” and using a human arguments xml file. The genome map was then further refined by re-assembling all data, but using the first assembly version as a reference. The final assemblies were both 1.3 Gb in total size, with an average coverage of 92.3 and 96.4×, and N50 of 0.93 and 0.95 Mb, for the northern and southern sample, respectively.Linked read sequencingFor the southern sample and sample DD81063, DNA for chromium sequencing (10× Genomics, CA, USA) was extracted from blood stored in SET buffer using a MagAttract HMW DNAkit (Qiagen) at Scilifelab, Stockholm, Sweden. For the northern sample the extraction for bionano optical maps was used. The libraries of the northern and southern sample were each sequenced on a separate lane of a HiSeqX (Illumina, CA, USA) and the DD81063 sample was sequenced on a NovaSeq6000 (Illumina). For all samples sequencing was performed using a 2 × 150 bp setup.Northern willow warbler de novo assemblyLibrary preparation for long read sequencing was done on DNA previously extracted for the optical map and followed Pacific Bioscience’s (CA, USA) standard protocol for 10–20 kb libraries. No shearing was performed prior to the library construction, but the library was size selected using the BluePippin pulse field size selection system (Sage Science, MA, USA), with a size cut-off >25 kb. The library was sequenced on eight SMRT cells on a Sequel platform (Pacific Biosciences). The sequencing yielded 63.66 Gbp of data comprised of 4,690,365 subreads with a mean length of 13,573 bp (range: 50–170,531 bp).The Pacbio reads were assembled de novo in HGAP449 in the SMRT Link package with default settings except for specifying an expected genome size of 1.2 Gbp and setting the polishing algorithm to “Arrow”. We ran Falcon unzip50 on the assembly to obtain partially phased primary contigs and fully phased haplotigs. Within the software, Arrow was used to polish the assembly using reads assigned to each haplotype. We evaluated two unzipped assemblies based on 30× or 40× coverage of seed reads in the preassembly step in HGAP4. A lower coverage threshold will lead to longer reads in the initial assembly step, which may increase the contiguity of the assembly, but will on the other hand, limit the number of reads that can be used in the phasing and polishing step. Although the unzipped assemblies were very similar, the 40× version was chosen for downstream analyses as it was slightly more contiguous and contained a higher number of single-copy bird orthologues as determined by BUSCO version 3.0.251.The assembly was further polished with Pilon 1.2252 with Illumina chromium reads from the same sample. The Illumina reads were mapped to the assembly using bwa version 0.7.17-r118853 and duplicated reads were marked using picardtools 2.10.3 (http://broadinstitute.github.io/picard). Pilon was run by only correcting indels and in total the software made 1,043,827 insertions and 275,457 deletions, respectively, of which the vast majority (94%) were single basepair changes. The Illumina polishing had a pronounced effect on the number of single-copy bird orthologues that could be detected in the primary contigs (Supplementary Table 1).For further assembly steps, we extracted the Illumina-polished primary Pacbio contigs (N = 2737, N50 of 2.1 Mb and a length of 1.29 Gb). These contigs showed an unexpectedly high level of duplicated single-copy orthologues (7.4%), which suggested partial or complete overlap between some contigs. As a first step to reduce the redundancy and increase the contiguity of the assembly, we hybridized the primary contigs to the optical map of the same sample using bionano solve version 3.2.2 (BioNano Genomics) with default settings except for specifying aggressive scaffolding parameters. The hybrid scaffolding resulted in 19 cuts to the bionano maps and 259 cuts to the Pacbio contigs and created 363 super-scaffolds. Most of the gaps between the contigs in the super-scaffolds were estimated to be negative (i.e., some overlap between sequences). However, in the hybrid assembly, sequences on either side of these gaps were not collapsed and thus formed false segmental duplications. To remedy this problem we extracted 304 sets of overlapping contigs (“supercontigs”) and used GAP5 in the staden package 2.0.0.b1154 to find potential joins between the contig ends. Using this approach, we merged contigs at 558 (87%) of the putative overlaps. The mean alignment length in the overlaps was 111 kb (range: 0.259–661 kb) with a mean sequence divergence of 3.28% (range: 0.31–15.55%). The highest divergence was caused by the presence of large indels. By trimming off one or both ends of the contigs at the gaps (mean 23 kb, range: 0.6–60 kb), we were able to close 23 further gaps. For the remainder of gaps, GAP5 failed to find potential joins between contigs or the ends supposed to be joined were considered to have too high divergence. The new assembly, including supercontigs consisted of 2401 contigs with an N50 of 6.5 Mb and had a considerably lower amount of duplicated single-copy genes (4.6% vs 7.4%).To further reduce the redundancy, we used the purge haplotig pipeline55 (downloaded 2019-02-15) to remove contigs that could be mapped over most of their length to larger contigs and that showed limited diploid coverage. We first estimated coverage by mapping the Pacbio subreads used for the de novo assembly with minimap2 version2.13-r86056 using default settings for Pacbio reads (-x map-pb). To minimize the loss of repetitive sequences that could be separated and scaffolded by the bionano optical map, we used the first bionano hybrid assembly (363 superscaffolds and 1500 cut and unscaffolded contigs) as a reference for mapping. From the mapped data we detected a clear haploid and diploid peak and set a threshold of diploid coverage above 34× and below 85×. Any scaffold where less than 80% of its positions had diploid coverage was considered a putative haplotig and was mapped to other scaffolds using minimap2 within the software. We removed 1209 scaffolds (mean size: 107,655 bp, range: 598–495,788 bp) with a coverage to the best hit of at least 70% (mean: 97.4%). Using this approach, we specifically excluded contigs that could not be incorporated in superscaffolds. However, we also removed three contigs that each entirely made up short superscaffolds that could be uniquely assigned to larger superscaffolds and that had a high degree of haploid coverage. At this stage, we also removed five additional contigs shorter than 1000 bp that were the result of cutting the assembly with the bionano optical map. This led to an assembly with 1187 contigs, a length of 1.1 Gbp and a N50 of 7.9 Mb. The filtered assembly showed a large reduction in single-copy orthologue bird genes (1.3 vs 4.6%).To provide an intermediate level of scaffolding to the optical map, we mapped the 10× chromium reads of the same sample to the assembly using bwa and used arcs version 1.0.557 and LINKS version 1.8.658 for scaffolding. Arcs was run with default settings except for enabling gap size estimation (–dist_est) and LINKS was run by setting the number of supporting links to at least 5 (-l = 5) and the maximum link ratio between the two best contig pairs to 0.3 (-a = 0.3). The scaffolding resulted in 739 scaffolds with a N50 of 16.4 Mb and a length 1.12 Gb.As a final scaffolding step, we hybridized the 10× chromium-Pacbio scaffolds to the bionano optical map using the same settings as before. The hybrid scaffolding made 23 cuts to the optical map, 122 cuts to the scaffolds and resulted in 497 scaffolds with an N50 of 16.8 Mb. Two contigs representing the divergent region on chromosome 1 had been scaffolded together by arcs but were separated and not re-scaffolded with other sequences in the bionano hybrid assembly. Since the mismatched end of the optical map was short, located at a large gap, and the gene order is the same as seen in other bird genomes, we decided to keep the scaffold generated by arcs.For this round of hybrid scaffolding, there were 52 gaps that were estimated to be negative. Using the same approach as when creating supercontigs, we were able to close 10 of these gaps. We additionally closed gaps using PBJelly59 from PBSuite 15.8.24 with default settings except for specifying –spanOnly –capturedOnly”. The software filled 97 gaps, extended one end of 12 gaps, extended both ends of 18 gaps and overfilled 28 gaps (extended both ends but detected no overlap despite the extension is larger than the predicted gap).We further checked for potential misjoins between scaffolds that originate from different chromosomes. To this end, we used SatsumaSynteny 2.060 to produce whole-genome alignments between the assembly and the genomes of chicken (version GRCg6a) and zebra finch (version taeGut3.2.4), both downloaded from Ensembl (www.ensembl.org). Using this approach, we detected a scaffold that showed good alignments to both chromosomes 10 and 23 in both of the other species. We considered this join unlikely and decided to split the scaffold.Next, we performed a second round of polishing with the 10× chromium Illumina data from the same sample. For this round, since we had fewer than 500 scaffolds, we used the longranger 2.1.14 align pipeline61 to map reads in a barcode-aware way. Pilon was then run with the same settings as before and resulted in the correction of 417,032 indels, of which 78.7% were single-basepair changes. The second round of polishing considerably increased the number of single-copy bird orthologues that could be identified in the assembly (Supplementary Table 1).The mitochondrial genome was not found in the original Pacbio genome assembly. We obtained this genome by adding the complete mitochondrial sequence from a previous short-read assembly13. We then used bwa to map the 10× chromium reads from the northern sample to the assembly and extracted alignments on the mitochondrial sequence. Next, freebayes was used with a haploid setting to detect differences present in the aligned reads. The raw variant file was filtered with vcftools for sites with a quality less than 30 and for two intervals with excessive read coverage (possibly reads from unassembled NUMTs). The filtered variant file contained 11 substitutions and three indels, and was used with bcftools version 1.1462 to create a new mitochondrial reference.For the extraction and removal of sequences in the different assembly steps we used kentUtils 370 (https://github.com/ucscGenomeBrowser/kent). Summary statistics for each assembly (e.g., N50) were calculated using the assemblathon_stats.pl script63.Southern willow warbler and chiffchaff de novo assembliesThe southern willow warbler and the chiffchaff were each sequenced on two lanes on a Sequel II (Pacific Biosciences) using a high-fidelity (HiFi) setup. Sequencing libraries for the southern willow warbler was prepared from a previous extraction used for optical maps (see above), whereas for the chiffchaff, DNA was extracted from blood using a Nanobind extraction kit (Circulomics, MD, USA). The southern willow sample yielded 2,576,876 HiFi reads with a mean length 19,303 bp and representing 49.7 Gbp. The chiffchaff sample yielded 2,612,165 HiFi reads with a mean length of 19,829 bp and representing 51.8 Gbp.The HiFi reads were assembled de novo using hifiasm version 0.15.5-r35064 with default settings and primary contigs were selected for downstream analyses. For the chiffchaff hifiasm assembly, we removed the first 6 Mb part of a contig overlapping with another contig and removed a short interval at the end of a contig containing adaptor sequences. For the southern willow warbler, the primary contigs (N = 540, Supplementary Table 1) were hybridized to the optical map of the same sample using the same pipeline as for the northern sample. Although we had access to chromium data from the same sample, we did not include it to perform an intermediate scaffolding step (as we did for the northern willow warbler assembly) because the long-read assembly was already highly contiguous. The hybridization step made 39 cuts to the contigs and 20 cuts to the optical maps, resulting in an assembly with 111 superscaffolds and 439 non-scaffolded contigs. We decided to ignore an optical map-supported fusion of contigs that mapped to separate chromosomes in other bird species, as this fusion was made in a large repetitive region. We further excluded a 45 bp sequence resulting from the hybrid assembly cutting and masked four short intervals containing adaptor sequences. The assembly of the mitochondrion in the southern assembly followed the same pipeline as used for the northern assembly (see above). In this case, 10 substitutions and two indels were added to the mitochondrial sequence from the previous short-read assembly based on alignments of linked reads from the southern sample.Repeat annotationWe used Repeatmodeler version 1.0.865 for de novo identification of repeats in the southern assembly. The repeats detected by repeatmodeler were combined with 1,023 bird-specific repeats into a custom library. Next, we used repeatmasker version 4.0.766 with the custom library and by using a more sensitive search (-s flag) to annotate repeats in the genome. Bedtools v2.29.267, together with the annotated repeats, was used to create a softmasked version of the southern assembly, which was used in the gene annotation step. The same repeat library was also used to annotate repeats in the de novo assembly of the northern sample. For the chiffchaff assembly we used the same annotation approach as for the southern willow warbler, but included a species-specific library generated with repeatmodeler, and also included a tandem-repeat associated sequence associated with the divergent regions on chromosomes 1 and 3 from the willow warbler library. Intervals with tandem repeats in divergent regions were also analyzed with tandem repeats finder version 4.0.968 using default settings except for specifying a maximum period size of 2000 bp.Duplicated intervals within divergent scaffolds were identified with Minimap2 and subsequently aligned with EMBOSS Stretcher 6.6.0 (https://www.ebi.ac.uk/Tools/psa/emboss_stretcher/).RNA sequencingWe used total RNA extracted from whole brain from six samples used in an earlier study quantifying differential expression in migratory and breeding willow warblers69 (Supplementary Table 3). The quality of the RNA was checked with a Bioanalyzer version 2100 (Agilent, CA, USA). All of the extractions had a RNA Integrity Number (RIN) of at least > 7.10. RNA libraries for sequencing were prepared using a TruSeq Stranded mRNA Sample prep kit with 96 dual indexes (Illumina) according to the instructions of the manufacturer with the exception of automating the protocols using an NGS workstation (Agilent) and using purification steps as described in Lundin et al70. and Borgström et al71. The raw RNA data was trimmed using cutadapt version 1.872 within Trim Galore version 0.4.0 (https://github.com/FelixKrueger/TrimGalore) with default settings.We used Stringtie version 1.3.373 to create transcripts from the RNAseq data. These transcripts were not used directly in the generation of gene models, but used in the manual curation step as potential alternative transcripts. For the software, we first mapped the reads with Hisat2 version 2.1.074 using default settings for stranded sequence libraries and downstream transcript analyses.Gene annotationWe used Augustus version 3.2.375 to create gene models using hints provided from RNAseq data and protein data from other bird species. For the RNAseq data, we mapped the trimmed reads to the assembly using STAR version 2.7.9a76. Accessory scripts in the Augustus package were used to filter the alignments for paired and uniquely mapped reads and for extracting intron hints. We additionally generated coverage wig files for each strand from the filtered alignment file using the software stranded-coverage (https://github.com/pmenzel/stranded-coverage) and used these as input for the august wig2hints.pl to generate exonpart hints.For homology evidence, we downloaded a set of bird proteins from NCBI (https://www.ncbi.nlm.nih.gov/). This data set included 49,673 proteins from chicken, 41,214 proteins from zebra finch and 38,619 proteins from great tit. We also downloaded an additional dataset from Uniprot (www.uniprot.org) that consisted of 3175 manually reviewed bird proteins and 204 and 12,263 bird proteins that were not manually reviewed but supported by protein or transcript data, respectively. The protein data was mapped to the genome using exonerate version 2.4.077. We used the script align2hints.pl from braker 2.1.678 to generate CDSpart, intron, start and stop hints from the data.Augustus was run with species-specific parameters (see training Augustus below) and with default settings except for specifying “softmasking=true”, “–alternatives-from-evidence=true”, “–UTR = on”, “–gff3=on” and “–allow_hinted_splicesites=atac”. In the extrinsic configuration file, we changed the malus for introns from 0.34 to 0.001, which increases the penalty for predicted introns that are not supported by the extrinsic data (RNAseq and protein hints). The prediction resulted in 28,491 genes and 35,389 transcripts.The Augustus-derived gene models were assigned names based on overlap with synteny-transferred zebra finch genes. For this purpose, we used SatsumaSynteny with default settings to obtain whole-genome alignments between our assembly and the zebra finch genome version bTaeGut1.4.pri79. Based on the alignment, we used kraken80 (downloaded 2020-04-14) to transfer the zebra finch genome annotations (NCBI Release 106) to the willow warbler assembly. We then extracted the CDS from the Augustus gene models and the kraken genes and used bedtools intersect to quantify the overlap. The gene models were also searched against the longest translation of each of the chicken, zebra finch and great tit Parus major genes used as evidence for the gene prediction step and against 86,131 swissprot vertebrate proteins using blastp 2.5.0+81 with an E value threshold of 1e−5. Gene models that were not annotated through synteny were assigned a gene name based on the blast results. Protein domains in the gene models were annotated with interproscan v 5.30–69.082. To reduce the number of false positive predictions we removed 5697 genes that were not supported by synteny to zebra finch genes, showed no significant similarity to vertebrate proteins or did not contain any annotated protein domains.We used Webapollo 2.6.583 to manually curate gene models in the previously identified divergent chromosome regions and in other regions where differences were present. In the curation step, we specifically validated the support for the coding sequence and the UTR and also removed genes that were likely to be pseudogenes based on a truncated coding sequence compared to homologous genes in other vertebrates, had no support from synteny in other bird species and/or that were located in repeat-rich regions.Training AugustusWe used a previous repeat-masked short-read assembly13 and the trimmed RNAseq data used in this study to obtain species-specific parameters for Augustus. The RNAseq data was assembled into transcripts using Trinity version 2.0.284 to create a de novo and a genome-guided assembly that together were comprised of 1,929,396 transcripts. The genome-guided transcript assembly was based on RNAseq mapped to the genome using GSNAP version 2016-07-1185 with default settings. We used PASA version 2.0.286 to create high-quality transcripts, which were imported into Webapollo. To assess the completeness of the transcripts, we compared them to synteny-transferred models from the chicken genome using Kraken. We selected 1249 transcripts that appeared complete, were not overlapping with other genes and showed less than 80% amino acid similarity to another gene in the training set. From this set, we excluded 21 genes that were giving initial training errors, which gave us a training set of 1228 genes. This gene set was randomly split into 1028 training genes and 200 genes used for testing. For training, we used the optimize_augustus.pl script with default settings except for the flag –UTR = on.Whole-genome resequencing and variant callingWe used the whole-genome resequencing data from nine samples of each migratory phenotype provided in Lundberg et al13. and sequenced an additional two high-coverage samples from each migratory phenotype (Supplementary Table 4). Sequencing libraries for the new samples were prepared with a TruSeq DNA PCR-Free kit (Illumina) with a targeted insert size of 670 bp or with a Truseq DNA nano (Illumina) with a targeted insert size of 350 bp. All of the new samples were sequenced on a HiSeqX (Illumina). The raw reads were trimmed with trimmomatic 0.3687 with the parameters “ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:30”.Quality-trimmed reads were mapped to the southern assembly using bwa mem with default settings except for specifying -M flag to ensure compatibility with the downstream duplicate removal steps and converted into binary alignment map (bam) files using samtools. For samples sequenced across multiple lanes, reads from each lane were mapped independently and the resulting bam files were merged with samtools. Read duplicates were removed with the markduplicates tool provided in picardtools.From the aligned whole-genome resequencing data set, we called variants with freebayes v1.1.0 using default settings and parallelizing the analyses of separate scaffolds using GNU parallel88. Vcflib version 2017-04-0489 was used to filter the raw set of variants for sites with quality score >30 and for alternate alleles that were supported by at least one read on each strand (SAF  > 0 & SAR  > 0) and had at least one read balanced to the right and the left (RPL  > 0 & RPR  > 0). Next, we used vcftools 0.1.1690 to filter genotypes with a coverage of at least 5x and removed sites a maximum of four genotypes missing in each of the populations. The variants were also filtered for collapsed repeats by removing sites with a mean coverage of more than twice the median mean coverage (30×). We next used vcflib to decompose haplotype calls and complex alleles into indels and SNPs and removed any variants that were overlapping with annotated repeats. This gave us a final of 51 million variants of which 45 million were bi-allelic SNPs. We used vcftools to calculate FST91 for each variant and for bi-allelic SNPs in non-overlapping windows of 10 kb. As many rare variants segregate in the willow warbler populations, which may downwardly bias differentiation estimates92, we focused on variants with a minor allele frequency of at least 0.1.Coverage for each resequenced sample was calculated in non-overlapping 1 kb windows using bedtools and only included properly paired reads with a mapping quality of at least 1. The raw coverage values for each sample were normalized by its median coverage across all windows.Structural variant callingWe used a combination of delly 0.9.193 and GraphTyper 2.7.494 to call structural variants in the resequenced samples. To identify a set of high confidence variants, we first mapped the long reads from the northern willow warbler to the southern assembly using minimap 2.22-r110156 with default settings for Pacbio reads and from the alignments called variants using delly. Next, GraphTyper was used to genotype the resequenced samples for the delly variants in the scaffolds containing the divergent chromosome regions. The raw set of variants were filtered to contain only sites with a “PASS” flag and, for each variant, the aggregated genotype, which is the genotype model out of breakpoint alignments and coverage that has the highest genotyping quality, was chosen for downstream analyses. Genetic differentiation (FST) was calculated in vcftools and variants with FST ≥ 0.7 between homozygotes in each divergent chromosome region were extracted and checked for overlap with genes and gene features using bedtools. To get more reliable differentiation estimates, we only included sites where at least 80% of the southern and northern homozygotes had genotypes.Inversion genotypes for resequenced samplesThe resequenced samples were assigned a genotype of southern and northern haplotypes for each of the divergent regions based on a multidimensional scaling (MDS)-based clustering in invclust95 of SNP array genotypes in Lundberg et al.13. To obtain genotypes of the SNPs included on the array in the resequenced samples, we mapped the SNP array probe sequences to the northern assembly using gmap and from the alignments extracted the positions of the focal SNPs. Next, we used freebayes to genotype the resequenced samples for these positions and plink version 1.996 to combine the genotypes with the genotypes from the SNP array. In the genotyping step, we also included mapped 10× chromium libraries for the northern and southern reference samples and the additional willow warbler sample. From the combined dataset, we extracted genotypes for SNPs located in each of the divergent regions and used invclust to assign each sample a genotype of inverted and non-inverted haplotypes. The inverted and non-inverted haplotypes were recoded as southern or northern haplotypes based on their frequency in each subspecies.Breakpoint analysesWe used MUMmer 4.0.0rc197 to align the genomes of the southern and northern willow warblers, and the southern willow warbler genome to the genomes of the chiffchaff, zebra finch (3.2.4) and collared flycatcher FicAlb (1.5)98.To provide further evidence of breakpoints, we mapped the 10× chromium reads of each sample to both the northern and the southern assembly and called structural variants using the longranger wgs pipeline. For the southern genome, we selected the 499 largest scaffolds and concatenated the rest into a single scaffold to make it compatible with the software. We also checked for differences in linked read molecule coverage between the samples. For this purpose, the raw reads of each sample were first processed with longranger basic for quality trimming and barcode processing. The trimmed reads were mapped to the assemblies using bwa mem using a -C flag to extract the barcode information of each read and alignments converted into bam files using samtools. To estimate coverage of barcodes, we first used the tigmint-molecule script from tigmint 1.1.299 to obtain positional information of barcodes (molecules) in each divergent region. The software was run with default settings except for only using reads with a mapping quality of at least 1 and only to report molecules that were estimated to be at least 10 kb. We next used bedtools to count the number of overlapping molecules in 1 kb windows.We explored differences between optical maps by using the runSV.py script in bionano solve with the southern optical map as a query and the northern assembly as target and the reciprocal analysis with the northern optical map as a query and the southern assembly as a target. We also used the bionano solve hybrid assembly pipeline to visualize differences between the optical maps and the genome assemblies at breakpoint regions.Functional annotation of differencesWe used bedtools to quantify the distance between breakpoint intervals and annotated genes. To provide a functional annotation of the SNPs and short indels, we selected variants that showed a FST ≥ 0.7 between southern and northern homozygotes for each of the region and used these as input to Snpeff 5.0.0e100 together with the annotation and reference genome. We used Snpsift 5.0.0e101 to select variants that were predicted to have a moderate to high effect on genes. Gene ontology terms for the genes were extracted from orthologous genes in other bird genomes in ensembl (www.ensembl.org) or through domain searches of the proteins with interproscan.Age estimation and demographic analyses of divergent regionsIn order to estimate the timing of the inversion events, we used high-coverage resequencing data from two southern samples, two northern samples and, as an outgroup, one dusky warbler Phylloscopus fuscatus (Supplementary Table 4). The willow warbler samples were chosen so that they were either homozygous southern or northern for all of three divergent regions. The dusky warbler library was prepared using a TruSeq Nano DNA library prep kit for Neoprep (Illumina) according to the instructions of the manufacturer and sequenced on a HiSeq X (Illumina). Quality-trimming of the raw reads and mapping of the trimmed reads to the northern reference genome followed the same approach as used for the willow warbler resequencing samples (see above).Variants were called using freebayes and the raw set of variants were filtered using gIMble’s preprocess module (v0.6.0). Sample-specific callable sites were identified using gIMble preprocess and were defined as those with a minimum coverage of 8× and a maximum of 0.75 standard deviations above the mean coverage. Genic and repetitive regions of the genome were removed from the callable sites in order to limit downstream analyses to intergenic regions.Summary statistics of genetic variation (π and dxy) within the divergent regions were calculated using gIMble. Following this, net divergence (da) between northern and southern samples was calculated as dnorth–south − (πnorth + πsouth)/2. To convert the net divergence into years we used the germline mutation rate (4.6 × 10−9) estimated in the collared flycatcher21. Relative node depth (RND) using the dusky warbler (DW) as an outgroup was calculated as dnorth–south/(dDW-north + dDW-south)/2. For each divergent region, a blockwise site frequency spectrum (bSFS) was generated with gIMble using blocks of 64 bp in length. This length refers to the number of callable sites within a block, while the physical length of blocks was allowed to vary due to missing data but was limited to 128 bp. Downstream analyses that relied on a bSFS used a kmax of 2, meaning that only marginal probabilities were calculated for mutation counts >2. The composite likelihood (CL) of a model, given the bSFS of one of the divergent regions, was optimized using the Nelder-Mead algorithm with the maximum number of iterations set to 1000. Within the software we evaluated three different population models. The first model was a strict isolation model (SI), with parameters ancestral effective population size, effective population sizes for southern and northern willow warblers and divergence time. The second model was an isolation with migration model (IM1) that also included a migration rate from northern to southern samples, and the third model (IM2) instead had a migration rate from southern to northern willow warblers.Simulations were carried out by msprime 0.7.4102 through gIMble. The recombination rates used for these simulations were chromosome-specific estimates from a high-density recombination map of the collared flycatcher98 and were 2.04, 1.95, and 2.63 cM/Mb for chromosomes 1, 3, and 5, respectively. A total of 100 replicates were simulated for the optimized SI parameters of each region. These simulated bSFSs were then optimized under both an SI model as well as the best fitting IM model for that region. The improvement in CL between these models was used as a null distribution for testing whether improvements in CL observed for the real data were greater than expected given a history of no migration. For each parameter, we calculated 95% CI as Maximum Composite Likelihood (MCL) estimate ± 1.96 * standard deviation of simulations (Supplementary Table 7). As a result, our estimates of uncertainty are affected by the recombination rates that we assumed for simulations. We also used the results of simulations to quantify the potential bias in MCL estimates due to intra-block recombination (Supplementary Table 7). However, we did not attempt to correct for this bias as it is relatively small (e.g., the MCL divergence times are estimated to be biased upwards by 7, 24, and 10%) and our estimation of the bias itself is largely dependent on the recombination rates we assumed.MSMC224 was used to explore genome-wide changes in Ne through time. As input to the software, we used the callable intergenic bed file and filtered vcf file mentioned above, with the addition of further filtering the bed file to only include autosomal scaffolds ≥500 kb and excluding the divergent regions. The input files for MSMC2, i.e., an unphased set of heterozygous sites for each sample, were generated using the generate_multihetsep.py script from msmc-tools. MSMC2 was run with a starting ρ/μ of 1 for 30 expectation-maximum iterations. For both the demographic modeling and MSMC2, we used the collared flycatcher germline mutation rate21 and a generation time of 1.7 years11 to convert divergence times into years.To infer the effects of demographic events and selection, we also calculated several genetic summary statistics. To this end, we first imputed missing genotypes and inferred haplotypes for the filtered set of variants using beagle version 5.4103. From the full set of samples, we selected 10 and seven samples that were homozygous southern or northern for the three divergent regions, respectively, as determined from the MDS analysis (see above), and extracted bi-allelic SNPs. To identify ancestral and derived alleles, we extracted genotypes for the focal SNP positions from the aligned chiffchaff and dusky warblers reads using bcftools 1.1462 with the mpileup command. As a conservative approach, we considered any site with the presence of both the reference and alternate allele as heterozygous (regardless of their frequencies) and only included sites where the coverage was at least one-third of the mean coverage among all sites for each outgroup species. We next used a customized script to extract the sites from the original vcf files, and, if necessary, switch the reference and alternate allele and swap the genotypes accordingly. With the polarized genotype data, we used PopGenome 2.7.5104 to calculate Fay and Wu’s H and vcftools to get counts for the derived allele. We further used selscan 1.3.0105 to calculate XP-nsl106 between the southern and northern samples, Sweepfinder2107 to calculate a composite likelihood ratio (CLR) between a model where a selective sweep has had an effect on the allele frequency and a model based on the genome-wide allele frequency spectrum and used vcftools to calculate nucleotide diversity, Tajima’s D and linkage disequilibrium (D’).The use of the southern assembly as a reference could potentially lead to a mapping bias for reads from southern samples, particularly in regions of higher divergence between the subspecies. This, in turn, could have an effect on genetic summary statistics and demographic modeling estimates. To explore the effect of reference bias, we therefore also mapped the resequencing data to the northern assembly, performed variant calling and calculated nucleotide diversity and Tajima’s D in 10 kb windows. For the northern assembly, we also used the same demographic modeling as used for the southern assembly. Contrasting average genetic summary statistics and demographic parameter estimates, we found negligible differences between the two genome assemblies (Supplementary Table 10).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

Avdelas, L. et al. The decline of mussel aquaculture in the European Union: Causes, economic impacts and opportunities. Rev. Aquac. 13, 91–118. https://doi.org/10.1111/raq.12465 (2021).Article 

Google Scholar 
Tamburini, E., Turolla, E., Fano, E. A. & Castaldelli, G. Sustainability of Mussel (Mytilus galloprovincialis) farming in the Po River delta, northern Italy, based on a life cycle assessment approach. Sustainability 12, 3814. https://doi.org/10.3390/su12093814 (2020).Article 
CAS 

Google Scholar 
Shumway, S. E. et al. Shellfish aquaculture-In praise of sustainable economies and environments. J. World Aquacult. Soc. 34, 8–10 (2003).
Google Scholar 
Musella, M. et al. Tissue-scale microbiota of the Mediterranean mussel (Mytilus galloprovincialis) and its relationship with the environment. Sci. Total Environ. 717, 137209. https://doi.org/10.1016/J.SCITOTENV.2020.137209 (2020).Article 
ADS 
CAS 

Google Scholar 
Peharda, M., Župan, I., Bavčević, L., Frankić, A. & Klanjšček, T. Growth and condition index of mussel Mytilus galloprovincialis in experimental integrated aquaculture. Aquac. Res. 38, 1714–1720. https://doi.org/10.1111/J.1365-2109.2007.01840.X (2007).Article 

Google Scholar 
Sarà, G., Zenone, A. & Tomasello, A. Growth of Mytilus galloprovincialis (Mollusca, bivalvia) close to fish farms: A case of integrated multi-trophic aquaculture within the Tyrrhenian sea. Hydrobiologia 636, 129–136. https://doi.org/10.1007/S10750-009-9942-2/TABLES/4 (2009).Article 

Google Scholar 
Danovaro, R., Gambi, C., Luna, G. M. & Mirto, S. Sustainable impact of mussel farming in the Adriatic Sea (Mediterranean Sea): Evidence from biochemical, microbial and meiofaunal indicators. Mar. Pollut. Bull. 49, 325–333. https://doi.org/10.1016/j.marpolbul.2004.02.038 (2004).Article 
CAS 

Google Scholar 
Tancioni, L. et al. Anthropogenic threats to fish of interest in aquaculture: Gonad intersex in a wild population of thinlip grey mullet Liza ramada (Risso, 1827) from a polluted estuary in central Italy. Aquac. Res. 47(5), 1670–1674 (2016).Article 

Google Scholar 
Chary, K. et al. Integrated multi-trophic aquaculture of red drum (Sciaenops ocellatus) and sea cucumber (Holothuria scabra): Assessing bioremediation and life-cycle impacts. Aquaculture 516, 734621. https://doi.org/10.1016/j.aquaculture.2019.734621 (2020).Article 
CAS 

Google Scholar 
Purcell, S. W., Williamson, D. H. & Ngaluafe, P. Chinese market prices of beche-de-mer: Implications for fisheries and aquaculture. Mar. Policy 91, 58–65. https://doi.org/10.1016/j.marpol.2018.02.005 (2018).Article 

Google Scholar 
Morroni, L. et al. Sea cucumber Holothuria polii (Delle Chiaje, 1823) as new model for embryo bioassays in ecotoxicological studies. Chemosphere 240, 124819. https://doi.org/10.1016/j.chemosphere.2019.124819 (2020).Article 
ADS 
CAS 

Google Scholar 
Uthicke, S. & Karez, R. Sediment patch selectivity in tropical sea cucumbers (Holothuroidea: Aspidochirotida) analysed with multiple choice experiments. J. Exp. Mar. Biol. Ecol. 236, 69–87. https://doi.org/10.1016/S0022-0981(98)00190-7 (1999).Article 

Google Scholar 
MacTavish, T., Stenton-Dozey, J., Vopel, K. & Savage, C. Deposit-feeding sea cucumbers enhance mineralization and nutrient cycling in organically-enriched coastal sediments. PLoS ONE 7, 1–11. https://doi.org/10.1371/journal.pone.0050031 (2012).Article 
CAS 

Google Scholar 
Rakaj, A. et al. Towards sea cucumbers as a new model in embryo-larval bioassays: Holothuria tubulosa as test species for the assessment of marine pollution. Sci. Total Environ. 787, 147593. https://doi.org/10.1016/j.scitotenv.2021.147593 (2021).Article 
ADS 
CAS 

Google Scholar 
Purcell, S., Conand, C., Uthicke, S. & Byrne, M. Ecological roles of exploited sea cucumbers. Oceanogr. Mar. Biol. 54, 367–386. https://doi.org/10.1201/9781315368597-8 (2016).Article 

Google Scholar 
Zamora, L. N., Yuan, X., Carton, A. G., Slater, M. J. & Marine, L. Role of deposit-feeding sea cucumbers in integrated multitrophic aquaculture: Progress, problems, potential and future challenges. Rev. Aquac. 10, 57–74. https://doi.org/10.1111/raq.12147 (2016).Article 

Google Scholar 
Slater, M. J. & Carton, A. G. Survivorship and growth of the sea cucumber Australostichopus (Stichopus) mollis (Hutton 1872) in polyculture trials with green-lipped mussel farms. Aquaculture 272, 389–398. https://doi.org/10.1016/j.aquaculture.2007.07.230 (2007).Article 

Google Scholar 
Slater, M. J. & Carton, A. G. Effect of sea cucumber (Australostichopus mollis) grazing on coastal sediments impacted by mussel farm deposition. Mar. Pollut. Bull. 58, 1123–1129. https://doi.org/10.1016/j.marpolbul.2009.04.008 (2009).Article 
CAS 

Google Scholar 
Slater, M. J. & Carton, A. G. Sea cucumber habitat differentiation and site retention as determined by intraspecific stable isotope variation. Aquac. Res. 41, 695–702. https://doi.org/10.1111/j.1365-2109.2010.02607.x (2010).Article 
CAS 

Google Scholar 
Stenton-Dozey, J. Finding hidden treasure in aquaculture waste. Water Atmos. 15, 9–11 (2007).
Google Scholar 
Slater, M. J., Jeffs, A. G. & Carton, A. G. The use of the waste from green-lipped mussels as a food source for juvenile sea cucumber, Australostichopus mollis. Aquaculture 292, 219–224. https://doi.org/10.1016/j.aquaculture.2009.04.027 (2009).Article 

Google Scholar 
Stenton-Dozey, J. & Heath, P. A first for New Zealand: Culturing our endemic sea cucumber for overseas markets. Water Atmos. 17, 20–21 (2009).
Google Scholar 
Zamora, L. N. & Jeffs, A. G. Feeding, selection, digestion and absorption of the organic matter from mussel waste by juveniles of the deposit-feeding sea cucumber, Australostichopus mollis. Aquaculture 317, 223–228. https://doi.org/10.1016/j.aquaculture.2011.04.011 (2011).Article 

Google Scholar 
Zamora, L. N. & Jeffs, A. G. The ability of the deposit-feeding sea cucumber Australostichopus mollis to use natural variation in the biodeposits beneath mussel farms. Aquaculture 326, 116–122. https://doi.org/10.1016/J.AQUACULTURE.2011.11.015 (2012).Article 

Google Scholar 
Zamora, L. N. & Jeffs, A. G. A Review of the research on the Australasian Sea Cucumber, Australostichopus mollis (Echinodermata: Holothuroidea) (Hutton 1872), with emphasis on aquaculture. J. Shellfish Res. 32, 613–627. https://doi.org/10.2983/035.032.0301 (2013).Article 

Google Scholar 
Zamora, L. N. & Jeffs, A. G. Macronutrient selection, absorption and energy budget of juveniles of the Australasian sea cucumber, Australostichopus mollis, feeding on mussel biodeposits at different temperatures. Aquac. Nutr. 21, 162–172. https://doi.org/10.1111/ANU.12144 (2015).Article 
CAS 

Google Scholar 
Chatzivasileiou, D. et al. An IMTA in Greece: Co-culture of fish, bivalves, and holothurians. J. Mar. Sci. Eng. 10, 776. https://doi.org/10.3390/jmse10060776 (2022).Article 

Google Scholar 
Rakaj, A. et al. Spawning and rearing of Holothuria tubulosa: A new candidate for aquaculture in the Mediterranean region. Aquac. Res. 49, 557–568. https://doi.org/10.1111/are.13487 (2018).Article 
CAS 

Google Scholar 
Rakaj, A., Fianchini, A., Boncagni, P., Scardi, M. & Cataudella, S. Artificial reproduction of Holothuria polii: A new candidate for aquaculture. Aquaculture 498, 444–453. https://doi.org/10.1016/j.aquaculture.2018.08.060 (2019).Article 

Google Scholar 
González-Wangüemert, M., Aydin, M. & Conand, C. Assessment of sea cucumber populations from the Aegean Sea (Turkey): First insights to sustainable management of new fisheries. Ocean Coast. Manag. 92, 87–94. https://doi.org/10.1016/J.OCECOAMAN.2014.02.014 (2014).Article 

Google Scholar 
González-Wangüemert, M., Valente, S. & Aydin, M. Effects of fishery protection on biometry and genetic structure of two target sea cucumber species from the Mediterranean Sea. Hydrobiologia 743, 65–74. https://doi.org/10.1007/s10750-014-2006-2 (2015).Article 

Google Scholar 
González-Wangüemert, M., Domínguez-Godino, J. A. & Cánovas, F. The fast development of sea cucumber fisheries in the Mediterranean and NE Atlantic waters: From a new marine resource to its over-exploitation. Ocean Coast. Manag. 151, 165–177. https://doi.org/10.1016/j.ocecoaman.2017.10.002 (2018).Article 

Google Scholar 
González-Wangüemert, M. & Godino, J. Sea cucumbers as new marine resource in Europe. Front. Mar. Sci. 3, 112 (2016).
Google Scholar 
Domínguez-Godino, J. A., Slater, M. J., Hannon, C. & González-Wangüermert, M. A new species for sea cucumber ranching and aquaculture: Breeding and rearing of Holothuria arguinensis. Aquaculture 438, 122–128. https://doi.org/10.1016/J.AQUACULTURE.2015.01.004 (2015).Article 

Google Scholar 
Günay, D., Emiroğlu, D., Tolon, T., Özden, O. & Saygi, H. Growth and survival rate of Juvenile Sea Cucumbers (Holothuria tubulosa, Gmelin, 1788) at Various Temperatures. Turk. J. Fish. Aquat. Sci. 15, 533–541. https://doi.org/10.4194/1303-2712-v15_2_41 (2015).Article 

Google Scholar 
Tolon, T. Effect of salinity on growth and survival of the juvenile sea cucumbers Holothuria tubulosa (Gmelin, 1788) and Holothuria poli (Delle Chiaje, 1923). Fresenius Environ. Bull. 26, 3930–3935 (2017).CAS 

Google Scholar 
Tolon, T., Emiroğlu, D., Günay, D. & Hancı, B. Effect of stocking density on growth performance of juvenile sea cucumber Holothuria tubulosa (Gmelin, 1788). Aquac. Res. 48, 4124–4131. https://doi.org/10.1111/are.13232 (2017).Article 

Google Scholar 
Tolon, M. T., Emiroglu, D., Gunay, D. & Ozgul, A. Sea cucumber (Holothuria tubulosa Gmelin, 1790) culture under marine fish net cages for potential use in integrated multi-trophic aquaculture (IMTA). Indian J. Geol. Mar. Sci. 46, 749–756 (2017).
Google Scholar 
Neofitou, N. et al. Contribution of sea cucumber Holothuria tubulosa on organic load reduction from fish farming operation. Aquaculture 501, 97–103. https://doi.org/10.1016/j.aquaculture.2018.10.071 (2019).Article 

Google Scholar 
Sadoul, B. et al. Aquaculture Is Holothuria tubulosa the golden goose of ecological aquaculture in the Mediterranean Sea? Aquaculture 554, 738149. https://doi.org/10.1016/j.aquaculture.2022.738149 (2022).Article 
CAS 

Google Scholar 
Cutajar, K. et al. Culturing the sea cucumber Holothuria poli in open-water integrated multi-trophic aquaculture at a coastal Mediterranean fish farm. Aquaculture 550, 737881. https://doi.org/10.1016/j.aquaculture.2021.737881 (2022).Article 
CAS 

Google Scholar 
Grosso, L. et al. Integrated Multi-Trophic Aquaculture (IMTA) system combining the sea urchin Paracentrotus lividus, as primary species, and the sea cucumber Holothuria tubulosa as extractive species. Aquaculture 534, 736268. https://doi.org/10.1016/J.AQUACULTURE.2020.736268 (2021).Article 
CAS 

Google Scholar 
González-Wangüemert, M., Valente, S., Henriques, F., Domínguez-Godino, J. A. & Serrão, E. A. Setting preliminary biometric baselines for new target sea cucumbers species of the NE Atlantic and Mediterranean fisheries. Fish. Res. 179, 57–66. https://doi.org/10.1016/J.FISHRES.2016.02.008 (2016).Article 

Google Scholar 
Aydin, M. Biometry, density and the biomass of the commercial sea cucumber population of the Aegean Sea. Turk. J. Fish. Aquat. Sci 19, 463–474. https://doi.org/10.4194/1303-2712-v19_6_02 (2018).Article 

Google Scholar 
Whitlock, M. C. & Schluter, D. Analisi Statistica dei Dati Biologici, Zanichelli (2010)Hammer, O. & Harper, D. A. T. PAST PAleontological STatistics Version 3 Reference Manual (University of Oslo, 2013).Zhou, Y. et al. Feeding and growth on bivalve biodeposits by the deposit feeder Stichopus japonicus Selenka (Echinodermata: Holothuroidea) co-cultured in lantern nets. Aquaculture 256, 510–520. https://doi.org/10.1016/j.aquaculture.2006.02.005 (2006).Article 

Google Scholar 
Pensa, D. et al. Population status, distribution and trophic implications of Pinna nobilis along the South-eastern Italian coast. Npj Biodivers. https://doi.org/10.21203/rs.3.rs-1425249/v1 (2022).Article 

Google Scholar 
Francour, P. Predation on holothurians: A literature review. Invert. Bio. 116, 52–60. https://doi.org/10.2307/3226924 (1997).Article 

Google Scholar 
Mecheta, A. & Mezali, K. A biometric study to determine the economic and nutritional value of sea cucumbers (Holothuroidea: Echinodermata) collected from Algeria’s shallow water areas. Beche-de-mer Inf. Bull. 39, 65–70 (2019).
Google Scholar 
Sun, J., Hamel, J. F., Gianasi, B. L., Graham, M. & Mercier, A. Growth, health and biochemical composition of the sea cucumber Cucumaria frondosa after multi-year holding in effluent waters of land-based salmon culture. Aquac. Environ. Interact. 12, 139–151. https://doi.org/10.3354/aei00356 (2020).Article 

Google Scholar 
Boncagni, P., Rakaj, A., Fianchini, A. & Vizzini, S. Preferential assimilation of seagrass detritus by two coexisting Mediterranean sea cucumbers: Holothuria polii and Holothuria tubulosa. Estuar. Coast. Shelf Sci. 231, 106464. https://doi.org/10.1016/j.ecss.2019.106464 (2019).Article 
CAS 

Google Scholar 
Rakaj, A., and Fianchini, A. Mediterranean sea cucumbers—Biology, ecology, and exploitation, Chapter. In The World of Sea Cucumbers Challenges, Advances, and Innovations (Mercier, A., Hamel, J.-F, Suhrbier, A. & Pearce, C.) (2023)Massin, C. & Jangoux, M. Observations écologiques sur Holothuria tubulosa, Holothuria poli et Holothuria forskali (Echinodermata, Holothuroidea) et comportement alimentaire de H. tubulosa. Référ. Cah. Biol. Mar. 17, 45–59 (1976).
Google Scholar 
Coulon, P. & Jangoux, M. Feeding rate and sediment reworking by the holothuroid Holothuria tubulosa (Echinodermata) in a Mediterranean seagrass bed off Ischia Island, Italy. Mar. Ecol. Progr. Ser. 92, 201–204 (1993).Article 
ADS 

Google Scholar 
Belbachir, N., Mezali, K. & Soualili, D. L. Selective feeding behaviour in some aspidochirotid holothurians (Echinodermata: Holothuroidea) at Stidia, Mostaganem Province, Algeria (2014).Grosso, L. et al. Trophic requirements of the sea urchin Paracentrotus lividus varies at different life stages: comprehension of species ecology and implications for effective feeding formulations. Front. Mar. Sci. 9, 865450. https://doi.org/10.3389/fmars.2022.865450 (2022).Article 

Google Scholar 
Sun, Z. L., Gao, Q. F., Dong, S. L., Shin, P. K. & Wang, F. Estimates of carbon turnover rates in the sea cucumber Apostichopus japonicus (Selenka) using stable isotope analysis: The role of metabolism and growth. Mar. Ecol. Prog. Ser. 457, 101–112. https://doi.org/10.3354/meps09760 (2012).Article 
ADS 

Google Scholar 
Yuan, X. T. et al. Effects of aestivation on the energy budget of sea cucumber Apostichopus japonicus (Selenka) (Echinaodermata: Holothuroidea). Acta. Ecol. Sin. 27, 3155−3161. https://doi.org/10.1016/S1872-2032(07)60070-5 (2007).Article 

Google Scholar 
Liu, Y., Dong, S. L., Tian, X. L., Wang, F. & Gao, Q. F. Effects ofdietary sea mud and yellow soil on growth and energybudget of the sea cucumber Apostichopus japonicas (Selenka). Aquaculture 286, 266–270. https://doi.org/10.1016/j.aquaculture.2008.09.029 (2009).Article 

Google Scholar 
Brown, N. P. & Eddy, S. D. Echinoderm Aquaculture (Wiley, 2015).Book 

Google Scholar 
Qiu, T., Zhang, L., Zhang, T., Bai, Y. & Yang, H. Effect of culture methods on individual variation in the growth of sea cucumber Apostichopus japonicus within a cohort and family. Chin. J. Oceanol. Limnol. 32, 737–742. https://doi.org/10.1007/S00343-014-3131-5 (2014).Article 
ADS 

Google Scholar 
Zappes, I. A. et al. New data on Weddell seal (Leptonychotes weddellii) colonies: A genetic analysis of a top predator from the Ross Sea, Antarctica. PLoS ONE 12, 0182922. https://doi.org/10.1371/journal.pone.0182922 (2017).Article 
CAS 

Google Scholar 
Paltzat, D. L., Pearce, C. M., Barnes, P. A. & McKinley, R. S. Growth and production of California sea cucumbers (Parastichopus californicus Stimpson) co-cultured with suspended Pacific oysters (Crassostrea gigas Thunberg). Aquaculture 275, 124–137. https://doi.org/10.1016/j.aquaculture.2007.12.014 (2008).Article 

Google Scholar 
Dong, S. et al. Intra-specific effects of sea cucumber (Apostichopus japonicus) with reference to stocking density and body size. Aquac. Res. 41, 1170–1178. https://doi.org/10.1111/J.1365-2109.2009.02404.X (2010).Article 

Google Scholar 
Pei, S., Dong, S., Wang, F., Gao, Q. & Tian, X. Effects of stocking density and body physical contact on growth of sea cucumber, Apostichopus japonicus. Aquac. Res. 45, 629–636. https://doi.org/10.1111/ARE.12004 (2014).Article 

Google Scholar 
Xia, B., Ren, Y., Wang, J., Sun, Y. & Zhang, Z. Effects of feeding frequency and density on growth, energy budget and physiological performance of sea cucumber Apostichopus japonicus (Selenka). Aquaculture 466, 26–32. https://doi.org/10.1016/J.AQUACULTURE.2016.09.039 (2017).Article 

Google Scholar 
Domínguez-Godino, J. A. & González-Wangüemert, M. Does space matter? Optimizing stocking density of Holothuria arguinensis and Holothuria mammata. Aquac. Res. 49, 3107–3115. https://doi.org/10.1111/are.13773 (2018).Article 

Google Scholar 
Rugnini, L., Rossi, C., Antonaroli, S., Rakaj, A. & Bruno, L. The influence of light and nutrient starvation on morphology, biomass and lipid content in seven strains of green microalgae as a source of biodiesel. Microorganisms 8, 1254. https://doi.org/10.3390/microorganisms8081254 (2020).Article 
CAS 

Google Scholar  More

Our results revealed statistically significant relationships between CTs and BA worldwide (P  More

Seto, K. C., Güneralp, B. & Hutyra, L. R. Global forecasts of urban expansion to 2030 and direct impacts on biodiversity and carbon pools. Proc. Natl. Acad. Sci. USA 109, 16083–16088. https://doi.org/10.1073/pnas.1211658109 (2012).Article 
ADS 

Google Scholar 
Faeth, S. H., Bang, C. & Saari, S. Urban biodiversity: Patterns and mechanisms. Ann. N. Y. Acad. Sci. 1223, 69–81. https://doi.org/10.1111/j.1749-6632.2010.05925.x (2011).Article 
ADS 

Google Scholar 
Elmqvist, T., Zipperer, W. & Güneralp, B. Urbanisation, habitat loss, biodiversity decline: Solution pathways to break the cycle. In Routledge Handbook of Urbanisation and Global Environmental Change (eds Seta, K. et al.) 139–151 (Routledge, 2016).
Google Scholar 
Dirzo, R. et al. Defaunation in the Anthropocene. Science 345, 401–406 (2014).Article 
ADS 
CAS 

Google Scholar 
Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS One 12, e0185809. https://doi.org/10.1371/journal.pone.0185809 (2017).Article 
CAS 

Google Scholar 
Wagner, D., Grames, E. M., Forister, M. L., Berenbaum, M. R. & Stopak, D. Insect decline in the Anthropocene: Death by a thousand cuts. Biological sciences 118, e2023989118. https://doi.org/10.1073/pnas.2023989118 (2021).Article 
CAS 

Google Scholar 
Goulson, D., Nicholls, E., Botias, C. & Rotheray, E. L. Bee declines driven by combined stress from parasites, pesticides, and lack of flowers. Science 347, 6229. https://doi.org/10.1126/science.1255957 (2015).Article 
CAS 

Google Scholar 
Ollerton, J. (2021) Pollinators & pollination: nature and society. Pelagic publishing.IPBES (2016). The assessment report of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services on pollinators, pollination and food production. potts, S.G., Imperatriz-Fonseca, V.L and Ngo, H.T. (eds). Secretariat of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services, Bonn, Germany. 552 pages.Mallinger, R. E. & Gratton, C. Species richness of wild bees, but not the use of managed honeybees, increases fruit set of a pollinator dependent crop. J. Appl. Ecol. 52, 323–330 (2015).Article 

Google Scholar 
Kremen, C., Williams, N. M. & Thorp, R. W. Crop pollination from native bees at risk from agricultural intensification. Proc. Natl. Acad. Sci. U.S.A. 99, 16812–16816 (2002).Article 
ADS 
CAS 

Google Scholar 
Winfree, R., Fox, J. W., Williams, N. M., Reilly, J. R. & Cariveau, D. P. Abundance of common species, not species richness, drives delivery of a real-world ecosystem service. Ecol. Lett. 18, 626–635 (2015).Article 

Google Scholar 
Soroye, P., Newbold, T. & Kerr, J. Climate change contributes to widespread declines among bumble bees across continents. Science 367, 685–688 (2020).Article 
ADS 
CAS 

Google Scholar 
Matteson, K. C., Ascher, J. S. & Langellotto, G. A. Bee richness and abundance in New York City urban gardens. Ann. Entomol. Soc. Am. 101(1), 140–150. https://doi.org/10.1603/0013-8746(2008)101[140:BRAAIN]2.0.CO;2 (2008).Article 

Google Scholar 
Carré, G. et al. Landscape context and habitat type as drivers of bee diversity in European annual crops. Agr. Ecosyst. Environ. 133(1–2), 40–47. https://doi.org/10.1016/j.agee.2009.05.001 (2009).Article 

Google Scholar 
Goulson, D., Lye, G. C. & Darvill, B. Decline and conservation of bumble bees. Ann. Rev. Entomol. 53, 191–208. https://doi.org/10.1146/annurev.ento.53.103106.093454 (2008).Article 
CAS 

Google Scholar 
Bates, A. J. et al. Changing bee and hoverfly pollinator assemblages along an urban-rural gradient. PLoS One 6(8), e23459. https://doi.org/10.1371/journal.pone.0023459 (2011).Article 
ADS 
CAS 

Google Scholar 
Deguines, N., Julliard, R., De Flores, M. & Fontaine, C. Functional homogenization of flower visitor communities with urbanisation. Ecol. Evol. 6(7), 1967–1976. https://doi.org/10.1002/ece3.2009 (2016).Article 

Google Scholar 
Larsson, M. Higher pollinator effectiveness by specialist than generalist flower-visitors of unspecialized Knautia arvensis (Dipsacaceae). Oecologia 146(3), 394–403. https://doi.org/10.1007/s00442-005-0217-y (2005).Article 
ADS 

Google Scholar 
Pataki, D. E. et al. Coupling biogeochemical cycles in urban environments: Ecosystem services, green solutions, and misconceptions. Front. Ecol. Environ. 9, 27–36. https://doi.org/10.1890/090220 (2011).Article 

Google Scholar 
Mentens, J., Raes, D. & Hermy, M. Green roofs as a tool for solving rainwater runoff problems in the urbanized 21st century?. Landscape Urban Plann. 77, 217–226. https://doi.org/10.1016/j.landurbplan.2005.02.010 (2006).Article 

Google Scholar 
Oberndorfer, E. et al. Green roofs as urban ecosystems: Ecological structures, functions and services. Bioscience 57, 823–834. https://doi.org/10.1641/B571005 (2007).Article 

Google Scholar 
Braaker, S., Ghazoul, J., Obrist, M. K. & Moretti, M. Habitat connectivity shapes urban arthropod communities: The key role of green roofs. Ecology 95, 1010–1021. https://doi.org/10.1890/13-0705.1 (2014).Article 
CAS 

Google Scholar 
Colla, S. R., Willis, E. & Packer, I. Can green roofs provide habitat for urban bees (Hymenoptera: Apidae)?. Cities and the Environment 2(1), 1–12 (2009).Article 

Google Scholar 
Tonietto, R., Fant, J., Ascher, J., Ellis, K. & Larkin, D. A comparison of bee communities of Chicago green roofs, parks and prairies. Landsc. Urban Plan. 103, 102–108 (2011).Article 

Google Scholar 
Ksiazek, K., Fant, J. & Skogen, K. An asssement of pollen limitation on Chicago green roofs. Landsc. Urban Plan. 107, 401–408 (2012).Article 

Google Scholar 
MacIvor, J. S. Building height matters: Nesting activity of bees and wasps on vegetated roofs. Israel J. Ecol. Evol. 62, 88–96. https://doi.org/10.1080/15659801.2015.1052635 (2015).Article 

Google Scholar 
Kratschmer, S., Kriechbaum, M. & Pachinger, B. Buzzing on top: Linking wild bee diversity, abundance and traits with green roof qualities. Urban Ecosyst. 21, 429–441 (2018).Article 

Google Scholar 
MacIvor, J. S., Ruttan, R. & Salehi, B. Exotics on exotics: Pollen analysis of urban bees visiting Sedum on a green roof. Urban Ecosyst. 18, 419–430 (2014).Article 

Google Scholar 
Matteson, K. C. & Langellotto, G. A. Determinates of inner city butterfly and bee species richness. Urban Ecosyst. 13, 333–347. https://doi.org/10.1007/s11252-010-0122-y (2010).Article 

Google Scholar 
Geslin, B., Gauzens, B., Thébault, E. & Dajoz, I. Plant pollinator networks along a gradient of urbanisation. PLoS One 8, e63421. https://doi.org/10.1371/journal.pone.0063421 (2013).Article 
ADS 

Google Scholar 
Baldock, K.C.R, et al. (2015) Where is the UK’s pollinator biodiversity? The importance of urban areas for flower-visiting insects. Proc. R. Soc. B. https://doi.org/10.1098/rspb.2014.2849Theodorou, P. et al. Urban fragmentation leads to lower floral diversity, with knock-on impacts on bee biodiversity. Sci. Rep. 10, 21756. https://doi.org/10.1038/s41598-020-78736-x (2020).Article 
ADS 
CAS 

Google Scholar 
Lowenstein, D.M., Matteson, K.C., Xiao, I., Silva, A.M. and Minor, E.S (2014) Humans, bees, and pollination services in the city: The case of Chicago, IL (USA). Biodiversity Conservation 1–18. https://doi.org/10.1007/s10531-014-0752-0Winfree, R., Bartomeus, I. & Cariveau, D. Native pollinators in anthropogenic habitats. Annu. Rev. Ecol. Evol. Syst. 42, 1–22 (2011).Article 

Google Scholar 
Cariveau, D. P. & Winfree, R. Causes of variation in wild bee responses to anthropogenic drivers. Curr. Opin. Insect. Sci. 10, 104–109. https://doi.org/10.1016/j.cois.2015.05.004 (2015).Article 

Google Scholar 
Baldock, K. C. R. et al. A systems approach reveals urban pollinator hotspots and conservation opportunities. Nat. Ecol. Evol. 3, 363–373. https://doi.org/10.1038/s41559-018-0769-y (2019).Article 

Google Scholar 
Li, W. C. & Yeung, K. K. A. A comprehensive study of green roof performance from environmental perspective. Int. J. Sustain. Built Environ. 3, 127–134 (2021).Article 

Google Scholar 
Turner, M., Baker, W. L., Peterson, C. J. & Peet, R. K. Factors influencing succession: Lessons from large, infrequent natural disturbances. Ecosystems 1, 511–523. https://doi.org/10.1007/s100219900047 (1998).Article 

Google Scholar 
Molineux, C. J., Connop, S. P. & Gange, A. C. Manipulating soil microbial communities in extensive green roof substrates. Sci. Total Environ. 493, 632–638. https://doi.org/10.1016/j.scitotenv.2014.06.045 (2014).Article 
ADS 
CAS 

Google Scholar 
Macivor, S. & Ksiazek, K. Invertebrates on green roofs. Ecol. Stud. Anal. Synthes. 223, 333–355. https://doi.org/10.1007/978-3-319-14983-7_14 (2015).Article 

Google Scholar 
Madre, F., Vergnes, A., Machon, N. & Clergeau, P. A comparison of 3 types of green roof as habitats for arthropods. Ecol. Eng. 57, 109–117. https://doi.org/10.1016/j.ecoleng.2013.04.029 (2013).Article 

Google Scholar 
Lee, L. H. & Lin, J. C. Green roof performance towards good habitat for butterflies in the compact city. Int. J. Biol. 7, 103. https://doi.org/10.5539/ijb.v7n2p103 (2015).Article 
CAS 

Google Scholar 
Preston, F. W. The canonical distribution of commonness and rarity: Part I. Ecology 43(2), 185–215. https://doi.org/10.2307/1931976 (1962).Article 

Google Scholar 
Orford, K. A., Murray, P. J., Vaughan, I. P. & Memmott, J. Modest enhancements to conventional grassland diversity improve the provision of pollination services. J. Appl. Ecol. 53, 906–915. https://doi.org/10.1111/1365-2664.12608 (2016).Article 

Google Scholar 
Brenneisen, S. The Natural Roof (NADA): Research Project Report on the Use of Extensive Green Roofs by Wild Bees (University of Wädenswil, 2005).
Google Scholar 
Jacobs, J., Berg, M., Beenaerts, N. & Artois, T. Biodiversity of Collembola on green roofs: A case study of three cities in Belgium. Ecol. Eng. 177, 106572. https://doi.org/10.1016/j.ecoleng.2022.106572 (2022).Article 

Google Scholar 
McKinney, M.L., Sisco, N.D. (2018) Systematic variation in roof spontaneous vegetation: residential “low rise” versus commercial “high rise” buildings. Urban Nature SI, 73–88.Rotheray, G.E., & Gilbert, S.F. (2011) The natural history of hoverflies. Tresaith, UK: Forrest TextBenvenuti, S. Wildflower green roofs for urban landscaping, ecological sustainability and biodiversity. Landsc. Urban Plan. 124, 151–161. https://doi.org/10.1016/j.landurbplan.2014.01.004 (2014).Article 

Google Scholar 
Schneider, F. Beitrag zur Kenntnis der Generationsverhaltnisse und Diapause rauberischer Schwebfliegen (Syrphldae, Dipt.). Mittl. Schweiz Ent Ges 21, 249–285 (1948).
Google Scholar 
Rader, R., Edwards, W., Westcott, D. A., Cunningham, S. A. & Howlett, B. G. Pollen transport differs among bees and flies in a human-modified landscape. Divers. Distrib. 17, 519–529. https://doi.org/10.1111/j.1472-4642.2011.00757.x (2011).Article 

Google Scholar 
Burgio, G. & Sommaggio, D. Syrphids as landscape bioindicators in Italian agroecosystems. Agr. Ecosyst. Environ. 120, 416–422 (2007).Article 

Google Scholar 
Doyle, T. et al. Pollination by hoverflies in the Anthropocene. Proc. R. Soc. B 287, 20200508. https://doi.org/10.1098/rspb.2020.0508 (2020).Article 

Google Scholar 
Persson, A. S., Ekroos, J., Olssona, P. & Smith, H. G. Wild bees and hoverflies respond differently to urbanisation, human population density and urban form. Landsc. Urban Plann. 204, 103901. https://doi.org/10.1016/j.landurbplan.2020.103901 (2020).Article 

Google Scholar 
Verboven, H., Uyttenbroeck, R., Brys, R. & Hermy, M. Different responses of bees and hoverflies to land use in an urban–rural gradient show the importance of the nature of the rural land use. Landsc. Urban Plan. 126, 31–41. https://doi.org/10.1016/j.landurbplan.2014.02.017 (2014).Article 

Google Scholar 
Schönrogge, K. et al. Host propagation permits extreme local adaptation in a social parasite of ants. Ecol. Lett. 9, 1032–1040 (2006).Article 

Google Scholar 
Schweiger, O. et al. Functional richness of local hoverfly communities (Diptera, Syrphidae) in response to land use across temperate Europe. Oikos 116, 461–472 (2007).Article 

Google Scholar 
KMI: Koninklijk Meteorologisch Instituut (2022) Analyse van het jaar 2020 en 2021. Available from https://www.meteobelgie.be/klimatologie/waarnemingen-en-analyses/jaar-2020/2274-jaa-2020 (2020) https://www.meteobelgie.be/klimatologie/waarnemingen-en-analyses/jaar-2021/2291-analyse-van-het-jaar-2021 (2021). Accessed on 12/05/2022.Shrestha, M. et al. Fluorescent pan traps affect the capture rate of insect orders in different ways. Insects 10(2), 40. https://doi.org/10.3390/insects10020040 (2019).Article 

Google Scholar 
Cooper, R., & Whitmore, R.C. (1990) Arthropod sampling methods in ornithology, Avian Foraging: theory, methodology, and applications. Studies in Avian Biology 13, Cooper Ornithological Society, California.Oberprieler, S. K., Andersen, A. & Braby, M. F. Invertebrate by-catch from vertebrate pitfall rraps can be useful for documenting patterns of invertebrate diversity. J. Insect. Conserv. 23(3), 547–554. https://doi.org/10.1007/s10841-019-00143-z (2019).Article 

Google Scholar 
Skvarla, M. J., Larson, J. L. & Dowling, A. P. G. Pitfalls and preservatives: A review. J. Entomol. Soc. Ontario 145, 15–43 (2014).
Google Scholar 
Michez, D., Rasmont, P., Terzo, M. and Vereecken, N.J. (2019) Bees of Europe. Hymenoptera of Europe 1. N.A.P Editions.Williams, P.H., et al. (2012): Unveiling cryptic species of the bumblebee subgenus Bombus s. str. worldwide with COI barcodes (Hymenoptera: Apidae). Syste. Biodiversity. https://doi.org/10.1080/14772000.2012.66457Falck, S., & Lewington, R (2020) Bijen veldgids voor Nederland en Vlaanderen. Tirion.Koster, A. (2022) De Nederlandse wilde bijen en hun planten. http://www.denederlandsebijen.nl/. Accessed on 21/4/2022.Speight, M.C.D. & Sarthou, J.P. (2013) StN keys for the identification of adult European Syrphidae (Diptera) 2013/Clés StN pour la détermination des adultes des Syrphidae Européens (Diptères) 2013. Syrph the Net, the database of European Syrphidae, Vol. 74, 133pp, Syrph the Net publications, Dublin.Roback, P., Legler, J. (2021) Beyond Multiple Linear Regression: Applied Generalized Linear Models and Multilevel Models in R. Taylor & Francis Group, LLC.R Core Team (2020) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.Oksanen, J., et al. (2014) Vegan: community ecology package. R Package 280.Bengtsson, H. (2017). matrixStats: Functions that Apply to Rows and Columns of Matrices (and to Vectors). R Package Version 0.52.2.Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015) Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01.Wickham, H., François, R., Henry, L. and Müller, K. (2022). dplyr: A Grammar of Data Manipulation. https://dplyr.tidyverse.org, https://github.com/tidyverse/dplyr.Venables, W.N., & Ripley, B.D. (2002) Modern Applied Statistics with S, 4th ed. Springer, New York. ISBN 0–387–95457–0. https://www.stats.ox.ac.uk/pub/MASS4/.Ricotta, C. & Moretti, M. CWM and Rao’s quadratic diversity: A unified framework for functional ecology. Oecologia 167(1), 181–188 (2011).Article 
ADS 

Google Scholar 
Leclère, D. et al. Bending the curve of terrestrial biodiversity needs an integrated strategy. Nature 585(7826), 551–556. https://doi.org/10.1038/s41586-020-2705-y (2020).Article 
ADS 
CAS 

Google Scholar 
Drossart, M., et al. (2019) Belgian red list of Bees. Belgian Science Policy (BRAIN-be – (Belgian Research Action through Interdisciplinary Networks). Mons: Presse universitaire de l’Université de Mons.Fahrig, L. Why do several small patches hold more species than few large patches?. Glob. Ecol. Biogeogr. 29, 615–628. https://doi.org/10.1111/geb.13059 (2020).Article 

Google Scholar 
Ayers, A. C. & Rehan, S. M. Supporting bees in cities: how bees are influenced by local and landscape features. Insects 12, 128. https://doi.org/10.3390/insects12020128 (2021).Article 

Google Scholar 
Domínguez, M. V. S., González, E., Fabián, D., Salvo, A. & Fenoglio, M. S. Arthropod diversity and ecological processes on green roofs in a semi-rural area of Argentina: Similarity to neighbor ground habitats and landscape effects. Landscape and Urban Planning 199, (2020).Castagneyrol, B. & Jactel, H. Unravelling plant- animal diversity relationships: A meta-regression analysis. Ecology 93(9), 2115–2124 (2012).Article 

Google Scholar 
Harrison, T., Gibbs, J. & Winfree, R. Phylogenetic homogenization of bee communities across ecoregions. Glob. Ecol. Biogeogr. 27, 1457–1466. https://doi.org/10.1111/geb.12822 (2018).Article 

Google Scholar 
Wenzel, A., Grass, I., Belavadi, V. V. & Tscharntke, T. How urbanisation is driving pollinator diversity and pollination, a systematic review. Biol. Conserv. 241, 108321. https://doi.org/10.1016/j.biocon.2019.108321 (2020).Article 

Google Scholar 
Martins, K. T., Gonzalez, A. & Lechowicz, M. J. Patterns of pollinator turnover and increasing diversity associated with urban habitats. Urban Ecosyst. 20, 1359–1371 (2017).Article 

Google Scholar 
Bucholz, S. & Egerer, M. Functional ecology of wild bees in cities: Towards a better understanding of trait-urbanisation relationships. Biodiver. Conserv. 29, 2779–2801 (2020).Article 

Google Scholar 
Hernandez, J. L., Frankie, G. W. & Thorp, R. W. Ecology of urban bees : A review of current knowledge and directions for future study. Cities Environ. 2, 1–15 (2009).Article 

Google Scholar 
Cane, J. H. Bees, pollination, and the challenges of sprawl. In Nature in fragments: The legacy of sprawl (eds Johnson, E. A. & Klemens, M. W.) 109–124 (Columbia University Press, 2005).Chapter 

Google Scholar 
Koch, K. Wilde bijensoorten in een stedelijke omgeving: Stad Antwerpen. Antenna 4, 8–12 (2014).
Google Scholar 
Soper, J. & Beggs, J. Assessing the impact of an introduced bee, Anthidium manicatum, on pollinator communities in New Zealand. NZ J. Bot. 51(3), 213–228. https://doi.org/10.1080/0028825X.2013.793202 (2013).Article 

Google Scholar 
Bennet, D.G., Kelly, D., & Clemens, J. (2018). Food plants and foraging distances for the native bee Lasioglossum sordidum in Christchurch Botanic Gardens. New Zealand J. Ecol. 42(1), 40–47. https://doi.org/10.20417/nzjecol.42.1Vanormelingen, P., Remer, M., & D’Haeseleer, J. (2021) Wilde bijen en bebouwing: meer verliezers dan winnaars? Themanummer bijen in de stad en dorp, Hymenovaria, maart 2021.Rader, R. et al. Alternative pollinator taxa are equally efficient but not as effective as the honey-bee in a mass flowering crop. J. Appl. Ecol. 46(5), 1080–1087. https://doi.org/10.1111/j.1365-2664.2009.01700.x (2009).Article 

Google Scholar 
Garantonakis, N. et al. Comparing the pollination services of honey bees and wild bees in a watermelon field. Sci. Hortic. 204, 138–144. https://doi.org/10.1016/j.scienta.2016.04.006 (2016).Article 

Google Scholar 
Foldesi, R., Howlett, B. G., Grass, I. & Batary, P. Larger pollinators deposit more pollen on stigmas across multiple plant species – A meta-analysis. J. Appl. Ecol. 58(4), 699–707. https://doi.org/10.1111/1365-2664.13798 (2021).Article 

Google Scholar 
Howlett, et al. (2011). Can insect body pollen counts be used to estimate pollen deposition on pak choi stigmas? New Zealand Plant Protection 64, 25–31. https://doi.org/10.30843/nzpp.2011.64.5951Nelson, W., Barry Donovan, L. E. & Howlett, B. Lasioglossum bees – the forgotten pollinators. J. Apic. Res. https://doi.org/10.1080/00218839.2022.2028966 (2022).Article 

Google Scholar 
Passaseo, A., Pétremand, G., Rochefort, S. & Castella, E. Pollinators emerging from extensive green roofs: Wild bees (Hymenoptera: Antophila) and hoverflies (Diptera: Syrphidae) in Geneva (Switzerland). Urban Ecosyst. 23, 1079–1086. https://doi.org/10.1007/s11252-020-00973-9 (2020).Article 

Google Scholar 
Odanaka, K. A. & Rehan, S. M. Impact indicators: Effects of land use management on functional trait and phylogenetic diversity of wild bees. Agric. Ecosyst. Environ. 286, 106663 (2019).Article 

Google Scholar 
Wilson, C. J. & Jamieson, M. A. The effects of urbanisation on bee communities depends on floral resource availability and bee functional traits. PLoS ONE 14(12), e0225852. https://doi.org/10.1371/journal.pone.0225852 (2019).Article 
CAS 

Google Scholar 
Osborne, J. L. et al. Quantifying and comparing bumblebee nest densities in gardens and countryside habitats. J. Appl. Ecol. 45, 784–792. https://doi.org/10.1111/j.1365-2664.2007.01359.x (2007).Article 

Google Scholar 
Glaum, P., Simao, M. C., Vaidya, C., Fitch, G. & Lulinao, B. Big city Bombus: Using natural history and land-use history to find significant environmental drivers in bumble-bee declines in urban development. R Soc Open Sci. 4, 170156. https://doi.org/10.1098/rsos.170156 (2017) (PMID: 28573023).Article 
ADS 

Google Scholar 
Rasmont, P. et al. Climatic risk and distribution atlas of European bumblebees. Biorisk 10, 1–246 (2015).Article 

Google Scholar 
Roger, N. et al. Impact of pollen resources drift on common bumblebees in NW Europe. Glob. Change Biol. 23, 68–76 (2017).Article 
ADS 

Google Scholar 
Frankie, G. W. et al. Ecological patterns of bees and their host ornamental flowers in two northern California cities. J. Kansas Entomol. Soc. 78, 227–246 (2005).Article 

Google Scholar 
Lerman, S. B. & Milam, J. Bee fauna and floral abundance within lawn-dominated suburban yards in Springfield, MA. Ann. Entomol. Soc. Am. 109, 713–723 (2016).Article 
CAS 

Google Scholar 
Braaker, S., Obrist, M. K., Ghazoul, J. & Moretti, M. Habitat connectivity and local conditions shape taxonomic and functional diversity of arthropods on green roofs. J. Anim. Ecol. 86, 521–531. https://doi.org/10.1111/1365-2656.12648 (2017).Article 

Google Scholar 
Passaseo, A., Rochefort, S., Pétremand, G., & Castella, E. (2021) Pollinators on green roofs: Diversity and trait analysis of wild bees (Hymenoptera: Anthophila) and Hoverflies (Diptera: Syrphidae) in an urban area (Geneva, Switzerland). Cities and the Environment (CATE) https://doi.org/10.15365/cate.2021.140201Hennig, E. & Ghazoul, J. Pollinating animals in the urban environment. Urban Ecosyst. 15, 149–166. https://doi.org/10.1007/s11252-011-0202-7 (2012).Article 

Google Scholar 
Mecke R. (1996) Die fauna begrünter dächer: Ökologische untersuchung verschiedener dachflächer im Hamburger stadtgebiet. University of Hamburg, Diploma dissertation.Bevk, D. The diversity of pollinators on green roofs. Acta Entomol. Slovenica 29(1), 5–14 (2021).
Google Scholar 
Speight, M.C.D. (2011) Species accounts of European Syrphidae (Diptera), Glasgow 2011. Syrph the Net, the database of European Syrphidae, vol. 65, 285 pp., Syrph the Net publications, Dublin.Wotton, K. R. et al. Mass seasonal migrations of hoverflies provide extensive pollination and crop protection services. Curr. Biol. 29, 2167–2173 (2019).Article 
CAS 

Google Scholar 
Boyer, K. J., Fragoso, F. P., Mabin, M. E. D. & Brunet, J. Netting and pan traps fail to identify the pollinator guild of an agricultural crop. Nat. Res. Sci. Rep. 10, 13819. https://doi.org/10.1038/s41598-020-70518-9 (2020).Article 
CAS 

Google Scholar  More

Rahman, M. A. et al. Comparing the infiltration potentials of soils beneath the canopies of two contrasting urban tree species. Urban For. Urban Green. 38, 22–32. https://doi.org/10.1016/j.ufug.2018.11.002 (2019).Article 

Google Scholar 
Zölch, T., Henze, L., Keilholz, P. & Pauleit, S. Regulating urban surface runoff through nature-based solutions – An assessment at the micro-scale. Environ. Res. 157, 135–144. https://doi.org/10.1016/j.envres.2017.05.023 (2017).Article 
CAS 

Google Scholar 
Barron, O. V., Barr, A. D. & Donn, M. J. Effect of urbanisation on the water balance of a catchment with shallow groundwater. J. Hydrol. 485, 162–176. https://doi.org/10.1016/j.jhydrol.2012.04.027 (2013).Article 
ADS 

Google Scholar 
Rosenzweig, B. R. et al. The value of urban flood modeling. Earth’s Future 9, e2020EF001739. https://doi.org/10.1029/2020EF001739 (2021).Article 
ADS 

Google Scholar 
Pauleit, S., Fryd, O., Backhaus, A. & Jensen, M. B. In Encyclopedia of Sustainability Science and Technology (ed. Meyers, R. A.) 1–29 (Springer, 2020).
Google Scholar 
Rahman, M. A. et al. Traits of trees for cooling urban heat islands: A meta-analysis. Build. Environ. 170, 106606. https://doi.org/10.1016/j.buildenv.2019.106606 (2020).Article 

Google Scholar 
Ziter, C. D., Pedersen, E. J., Kucharik, C. J. & Turner, M. G. Scale-dependent interactions between tree canopy cover and impervious surfaces reduce daytime urban heat during summer. Proc. Natl. Acad. Sci. USA 116, 7575–7580. https://doi.org/10.1073/pnas.1817561116 (2019).Article 
ADS 
CAS 

Google Scholar 
Waldrop, M. M. News feature: The quest for the sustainable city. Proc. Natl. Acad. Sci. 116, 17134–17138. https://doi.org/10.1073/pnas.1912802116 (2019).Article 
ADS 
CAS 

Google Scholar 
Cleugh, H. A., Bui, E., Simon, D., Xu, J. & Mitchell, V. G. The Impact of Suburban Design on Water Use and Microclimate (2005).Chan, F. K. S. et al. “Sponge City” in China—A breakthrough of planning and flood risk management in the urban context. Land Use Policy 76, 772–778. https://doi.org/10.1016/j.landusepol.2018.03.005 (2018).Article 

Google Scholar 
Morgan, R. P. C. Soil Erosion and Conservation (Wiley, 2005).
Google Scholar 
Xu, C. et al. Surface runoff in urban areas: The role of residential cover and urban growth form. J. Clean. Prod. 262, 121421. https://doi.org/10.1016/j.jclepro.2020.121421 (2020).Article 

Google Scholar 
Ostoić, S. K. & van den Bosch, C. C. K. Exploring global scientific discourses on urban forestry. Urban For. Urban Green. 14, 129–138. https://doi.org/10.1016/j.ufug.2015.01.001 (2015).Article 

Google Scholar 
Rahman, M. A. et al. Tree cooling effects and human thermal comfort under contrasting species and sites. Agric. For. Meteorol. 287, 107947. https://doi.org/10.1016/j.agrformet.2020.107947 (2020).Article 
ADS 

Google Scholar 
Rötzer, T., Rahman, M. A., Moser-Reischl, A., Pauleit, S. & Pretzsch, H. Process based simulation of tree growth and ecosystem services of urban trees under present and future climate conditions. Sci. Total Environ. 676, 651–664. https://doi.org/10.1016/j.scitotenv.2019.04.235 (2019).Article 
ADS 
CAS 

Google Scholar 
Grote, R. et al. Functional traits of urban trees: Air pollution mitigation potential. Front. Ecol. Environ. 14, 543–550. https://doi.org/10.1002/fee.1426 (2016).Article 

Google Scholar 
Pace, R. et al. A single tree model to consistently simulate cooling, shading, and pollution uptake of urban trees. Int. J. Biometeorol. 65, 277–289. https://doi.org/10.1007/s00484-020-02030-8 (2021).Article 
ADS 

Google Scholar 
Kuehler, E., Hathaway, J. & Tirpak, A. Quantifying the benefits of urban forest systems as a component of the green infrastructure stormwater treatment network. Ecohydrology https://doi.org/10.1002/eco.1813 (2017).Article 

Google Scholar 
Rahman, M. A., Moser, A., Gold, A., Rötzer, T. & Pauleit, S. Vertical air temperature gradients under the shade of two contrasting urban tree species during different types of summer days. Sci. Total Environ. 633, 100–111. https://doi.org/10.1016/j.scitotenv.2018.03.168 (2018).Article 
ADS 
CAS 

Google Scholar 
Rahman, M. A., Smith, J. G., Stringer, P. & Ennos, A. R. Effect of rooting conditions on the growth and cooling ability of Pyrus calleryana. Urban For. Urban Green. 10, 185–192. https://doi.org/10.1016/j.ufug.2011.05.003 (2011).Article 

Google Scholar 
Schellekens, J., Scatena, F. N., Bruijnzeel, L. A. & Wickel, A. J. Modelling rainfall interception by a lowland tropical rain forest in northeastern Puerto Rico. J. Hydrol. 225, 168–184. https://doi.org/10.1016/S0022-1694(99)00157-2 (1999).Article 
ADS 

Google Scholar 
Guevara-Escobar, A., González-Sosa, E., Véliz-Chávez, C., Ventura-Ramos, E. & Ramos-Salinas, M. Rainfall interception and distribution patterns of gross precipitation around an isolated Ficus benjamina tree in an urban area. J. Hydrol. 333, 532–541. https://doi.org/10.1016/j.jhydrol.2006.09.017 (2007).Article 
ADS 

Google Scholar 
Xiao, Q. F. & McPherson, E. G. Surface water storage capacity of twenty tree species in Davis, California. J. Environ. Qual. 45, 188–198. https://doi.org/10.2134/jeq2015.02.0092 (2016).Article 
CAS 

Google Scholar 
Xiao, Q. F., McPherson, E. G., Ustin, S. L. & Grismer, M. E. A new approach to modeling tree rainfall interception. J. Geophys. Res. Atmos. 105, 29173–29188. https://doi.org/10.1029/2000jd900343 (2000).Article 
ADS 

Google Scholar 
Carlyle-Moses, D. E. & Gash, J. H. C. In Forest Hydrology and Biogeochemistry: Synthesis of Past Research and Future Directions (eds Levia, D. F. et al.) 407–423 (Springer, 2011).Chapter 

Google Scholar 
Hirano, T. et al. The difference in the short-term runoff characteristic between the coniferous catchment and the deciduous catchment: The effects of storm size on storm generation processes of small forested catchment. J. Jpn. Soc. Hydrol. Water Resour. 22, 24–39. https://doi.org/10.3178/jjshwr.22.24 (2009).Article 

Google Scholar 
Chandler, K. R. & Chappell, N. A. Influence of individual oak (Quercus robur) trees on saturated hydraulic conductivity. For. Ecol. Manage. 256, 1222–1229. https://doi.org/10.1016/j.foreco.2008.06.033 (2008).Article 

Google Scholar 
Stewart, I. D. A systematic review and scientific critique of methodology in modern urban heat island literature. Int. J. Climatol. 31, 200–217. https://doi.org/10.1002/joc.2141 (2011).Article 

Google Scholar 
Beck, H. E. et al. Present and future Köppen-Geiger climate classification maps at 1-km resolution. Sci. Data 5, 180214. https://doi.org/10.1038/sdata.2018.214 (2018).Article 

Google Scholar 
Moreno-de las Heras, M., Nicolau, J. M., Merino-Martín, L. & Wilcox, B. P. Plot-scale effects on runoff and erosion along a slope degradation gradient. Water Resour. Res. 46, W04503. https://doi.org/10.1029/2009WR007875 (2010).Article 
ADS 

Google Scholar 
Wu, L., Peng, M., Qiao, S. & Ma, X.-Y. Effects of rainfall intensity and slope gradient on runoff and sediment yield characteristics of bare loess soil. Environ. Sci. Pollut. Res. 25, 3480–3487. https://doi.org/10.1007/s11356-017-0713-8 (2018).Article 

Google Scholar 
Rutter, A. J., Kershaw, K. A., Robins, P. C. & Morton, A. J. A predictive model of rainfall interception in forests, 1. Derivation of the model from observations in a plantation of Corsican pine. Agric. Meteorol. 9, 367–384. https://doi.org/10.1016/0002-1571(71)90034-3 (1971).Article 

Google Scholar 
Gash, J. H. C. An analytical model of rainfall interception by forests. Q. J. R. Meteorol. Soc. 105, 43–55. https://doi.org/10.1002/qj.49710544304 (1979).Article 
ADS 

Google Scholar 
Véliz-Chávez, C., Mastachi-Loza, C. A., Gonzalez-Sosa, E., Becerril-Pia, R. & Ramos-Salinas, N. M. Canopy storage implications on interception loss modeling. Am. J. Plant Sci. 05, 3032–3048. https://doi.org/10.4236/ajps.2014.520320 (2014).Article 

Google Scholar 
Fan, J., Oestergaard, K. T., Guyot, A. & Lockington, D. A. Measuring and modeling rainfall interception losses by a native Banksia woodland and an exotic pine plantation in subtropical coastal Australia. J. Hydrol. 515, 156–165. https://doi.org/10.1016/j.jhydrol.2014.04.066 (2014).Article 
ADS 

Google Scholar 
Ghimire, C. P., Bruijnzeel, L. A., Lubczynski, M. W. & Bonell, M. Rainfall interception by natural and planted forests in the Middle Mountains of Central Nepal. J. Hydrol. 475, 270–280. https://doi.org/10.1016/j.jhydrol.2012.09.051 (2012).Article 
ADS 

Google Scholar 
Pereira, F. L. et al. Modelling interception loss from evergreen oak Mediterranean savannas: Application of a tree-based modelling approach. Agric. For. Meteorol. 149, 680–688. https://doi.org/10.1016/j.agrformet.2008.10.014 (2009).Article 
ADS 

Google Scholar 
Pypker, T. G., Bond, B. J., Link, T. E., Marks, D. & Unsworth, M. H. The importance of canopy structure in controlling the interception loss of rainfall: Examples from a young and an old-growth Douglas-fir forest. Agric. For. Meteorol. 130, 113–129. https://doi.org/10.1016/j.agrformet.2005.03.003 (2005).Article 
ADS 

Google Scholar 
Ringgaard, R., Herbst, M. & Friborg, T. Partitioning forest evapotranspiration: Interception evaporation and the impact of canopy structure, local and regional advection. J. Hydrol. 517, 677–690. https://doi.org/10.1016/j.jhydrol.2014.06.007 (2014).Article 
ADS 

Google Scholar 
Price, A. G. & Carlyle-Moses, D. E. Measurement and modelling of growing-season canopy water fluxes in a mature mixed deciduous forest stand, southern Ontario, Canada. Agric. For. Meteorol. 119, 69–85. https://doi.org/10.1016/S0168-1923(03)00117-5 (2003).Article 
ADS 

Google Scholar 
Fathizadeh, O., Hosseini, S. M., Zimmermann, A., Keim, R. F. & Darvishi Boloorani, A. Estimating linkages between forest structural variables and rainfall interception parameters in semi-arid deciduous oak forest stands. Sci. Total Environ. 601–602, 1824–1837. https://doi.org/10.1016/j.scitotenv.2017.05.233 (2017).Article 
ADS 
CAS 

Google Scholar 
Livesley, S. J., Baudinette, B. & Glover, D. Rainfall interception and stem flow by eucalypt street trees—the impacts of canopy density and bark type. Urban For. Urban Green. 13, 192–197. https://doi.org/10.1016/j.ufug.2013.09.001 (2014).Article 

Google Scholar 
Xiao, Q. & McPherson, E. G. Rainfall interception by Santa Monica’s municipal urban forest. Urban Ecosyst. 6, 291–302. https://doi.org/10.1023/B:UECO.0000004828.05143.67 (2002).Article 

Google Scholar 
Rohatgi, A. WebPlotDigitizer (4.4), 2020).Team, R. C. (R Foundation for Statistical Computing, 2020).García-Palacios, P., Gross, N., Gaitán, J. & Maestre, F. T. Climate mediates the biodiversity–ecosystem stability relationship globally. PNAS 115, 8400–8405. https://doi.org/10.1073/pnas.1800425115 (2018).Article 
ADS 
CAS 

Google Scholar 
Le Provost, G. et al. Land-use history impacts functional diversity across multiple trophic groups. PNAS 117, 1573–1579. https://doi.org/10.1073/pnas.1910023117 (2020).Article 
ADS 
CAS 

Google Scholar 
El Kateb, H., Zhang, H., Zhang, P. & Mosandl, R. Soil erosion and surface runoff on different vegetation covers and slope gradients: A field experiment in Southern Shaanxi Province, China. CATENA 105, 1–10. https://doi.org/10.1016/j.catena.2012.12.012 (2013).Article 

Google Scholar 
Oliveira, P. T. S. et al. The water balance components of undisturbed tropical woodlands in the Brazilian cerrado. Hydrol. Earth Syst. Sci. 19, 2899–2910. https://doi.org/10.5194/hess-19-2899-2015 (2014).Article 
ADS 

Google Scholar 
Hümann, M. et al. Identification of runoff processes – The impact of different forest types and soil properties on runoff formation and floods. J. Hydrol. 409, 637–649. https://doi.org/10.1016/j.jhydrol.2011.08.067 (2011).Article 
ADS 

Google Scholar 
Sun, D. et al. Soil erosion and water retention varies with plantation type and age. For. Ecol. Manage. 422, 1–10. https://doi.org/10.1016/j.foreco.2018.03.048 (2018).Article 

Google Scholar 
Jost, G., Schume, H., Hager, H., Markart, G. & Kohl, B. A hillslope scale comparison of tree species influence on soil moisture dynamics and runoff processes during intense rainfall. J. Hydrol. 420–421, 112–124. https://doi.org/10.1016/j.jhydrol.2011.11.057 (2012).Article 

Google Scholar 
Sadeghi, S. M. M., Attarod, P., Van Stan, J. T. & Pypker, T. G. The importance of considering rainfall partitioning in afforestation initiatives in semiarid climates: A comparison of common planted tree species in Tehran, Iran. Sci. Total Environ. 568, 845–855. https://doi.org/10.1016/j.scitotenv.2016.06.048 (2016).Article 
ADS 
CAS 

Google Scholar 
Pretzsch, H. et al. Climate change accelerates growth of urban trees in metropolises worldwide. Sci. Rep. https://doi.org/10.1038/s41598-017-14831-w (2017).Article 

Google Scholar 
Rahman, M. A., Moser, A., Rötzer, T. & Pauleit, S. Microclimatic differences and their influence on transpirational cooling of Tilia cordata in two contrasting street canyons in Munich, Germany. Agric. For. Meteorol. 232, 443–456. https://doi.org/10.1016/j.agrformet.2016.10.006 (2017).Article 
ADS 

Google Scholar 
Nytch, C. J., Meléndez-Ackerman, E. J., Pérez, M. E. & Ortiz-Zayas, J. R. Rainfall interception by six urban trees in San Juan, Puerto Rico. Urban Ecosyst. 22, 103–115. https://doi.org/10.1007/s11252-018-0768-4 (2018).Article 

Google Scholar 
Rahman, M. A. et al. Comparative analysis of shade and underlying surfaces on cooling effect. Urban For. Urban Green. 63, 127223. https://doi.org/10.1016/j.ufug.2021.127223 (2021).Article 

Google Scholar 
Chen, L., Zhang, Z. & Ewers, B. E. Urban tree species show the same hydraulic response to vapor pressure deficit across varying tree size and environmental conditions. PLoS One https://doi.org/10.1371/journal.pone.0047882 (2012).Article 

Google Scholar 
Moser-Reischl, A., Rahman, M. A., Pauleit, S., Pretzsch, H. & Rötzer, T. Growth patterns and effects of urban micro-climate on two physiologically contrasting urban tree species. Landsc. Urban Plan. 183, 88–99. https://doi.org/10.1016/j.landurbplan.2018.11.004 (2019).Article 

Google Scholar 
Hao, M. et al. Impacts of changes in vegetation on saturated hydraulic conductivity of soil in subtropical forests. Sci. Rep. 9, 8372. https://doi.org/10.1038/s41598-019-44921-w (2019).Article 
ADS 
CAS 

Google Scholar 
Peters, E. B., McFadden, J. P. & Montgomery, R. A. Biological and environmental controls on tree transpiration in a suburban landscape. J. Geophys. Res. Biogeosci. https://doi.org/10.1029/2009jg001266 (2010).Article 

Google Scholar 
Komatsu, H., Kume, T. & Otsuki, K. Increasing annual runoff—broadleaf or coniferous forests?. Hydrol. Process. 25, 302–318. https://doi.org/10.1002/hyp.7898 (2011).Article 
ADS 

Google Scholar 
Li, X. et al. Process-based rainfall interception by small trees in Northern China: The effect of rainfall traits and crown structure characteristics. Agric. For. Meteorol. 218–219, 65–73. https://doi.org/10.1016/j.agrformet.2015.11.017 (2016).Article 
ADS 

Google Scholar 
Lukaszkiewicz, J. & Kosmala, M. Determining the age of streetside trees with diameter at breast height-based multifactorial model. Arboricult. Urban For. 34, 137–143. https://doi.org/10.48044/jauf.2008.018 (2008).Article 

Google Scholar 
Buttle, J. M. & Farnsworth, A. G. Measurement and modeling of canopy water partitioning in a reforested landscape: The Ganaraska Forest, southern Ontario, Canada. J. Hydrol. 466–467, 103–114. https://doi.org/10.1016/j.jhydrol.2012.08.021 (2012).Article 

Google Scholar 
Yang, B., Lee, D. K., Heo, H. K. & Biging, G. The effects of tree characteristics on rainfall interception in urban areas. Landsc. Ecol. Eng. 15, 289–296. https://doi.org/10.1007/s11355-019-00383-w (2019).Article 
CAS 

Google Scholar 
Klamerus-Iwan, A. & Witek, W. Variability in the Wettability and Water Storage Capacity of Common Oak Leaves (Quercus robur L). Water 10, 695. https://doi.org/10.3390/w10060695 (2018).Article 
CAS 

Google Scholar 
Van Stan, J. T., Siegert, C. M., Levia, D. F. & Scheick, C. E. Effects of wind-driven rainfall on stemflow generation between codominant tree species with differing crown characteristics. Agric. For. Meteorol. 151, 1277–1286. https://doi.org/10.1016/j.agrformet.2011.05.008 (2011).Article 
ADS 

Google Scholar 
Selbig, W. R. et al. Quantifying the stormwater runoff volume reduction benefits of urban street tree canopy. Sci. Total Environ. 806, 151296. https://doi.org/10.1016/j.scitotenv.2021.151296 (2022).Article 
ADS 
CAS 

Google Scholar 
Centre for Watershed Protection. Review of the Available Literature and Data on the Runoff and Pollutant Removal Capabilities of Urban Trees (Center for Watershed Protection, 2017).
Google Scholar 
Berland, A. et al. The role of trees in urban stormwater management. Landsc. Urban Plan. 162, 167–177. https://doi.org/10.1016/j.landurbplan.2017.02.017 (2017).Article 

Google Scholar 
Pauleit, S. Urban street tree plantings: Indentifying the key requirements. Proc. Inst. Civ. Eng. Municipal Eng. 156, 43–50. https://doi.org/10.1680/muen.2003.156.1.43 (2003).Article 

Google Scholar 
Weller, M. Tree Inventory Data of Central European Cities—Studies on the Composition and Structure of Urban Tree Populations and Derivation of Ecosystem Services. MSC thesis, Technical University of Munich, Germany (2021). More

Land vulnerability of an area is directly related to the natural as well as anthropogenic activities involved in the geomorphological unit. Being one of the most vulnerable ecosystems, the estuaries and estuarine islands are delicately affected by both ecological processes of the sea and land and have pressures from multiple anthropogenic stressors and global climate change42,43,44. Ecological vulnerability and ecological sensitivity are similar and both originated from the concept of ecotone10,45. The geomorphologic concept of landscape sensitivity was first proposed by Brunsden and Thornes, who argued that the sensitivity indicated the propensity to change and the capacity to absorb the effects of disturbances10,46,47. Landscape sensitivity is studied by many researchers such as Allison and Thomas, Miles et al., Harvey, Knox, Usher, Haara et al., Thomas, Jennings and Yuan Chi8,47,48,49,50,51,52,53,54, through different case studies. Based on their findings Yuan Chi summarized the important characteristics of the landscape sensitivity are: a, the change of the landscape ecosystem; it involves the change likelihood, ratio, and component, as well as the resistance and susceptibility to the change, b, the temporal and spatial scales; which determine the occurrence, degree, and distribution of the change, c, the external disturbances that cause the change; the disturbances included natural and anthropogenic origins with different categories and intensities, and d, the threshold of the landscape sensitivity; it refers to the point of transition for the landscape ecosystem8. The environmental vulnerability of the Munroe Island has been studied based on the characterization of the geomorphological and sociocultural dynamics of the region based on the above characteristics.Bathymetric surveys in Ashtamudi lake and the Kallada riverThe present study shows that the geomorphic processes occurring on the Munroe Island are affected by anthropogenic disturbances in the morpho-dynamics of the Kallada river, Ashtamudi backwaters and associated fluvio-tidal interactions. A detailed bathymetric survey of both water bodies up to the tidal-influenced upper limit of the Kallada river27 was conducted with 200 m spaced grid references (Fig. 5). Bathymetry shows that the deepest point of the Ashtamudi backwater system is in Vellimon lake (13.45 m), the SE extension of Ashtamudi lake. The eastern side of Ashtamudi lake is deeper than the western side of this backwater system. The depth of the backwater decreases towards the estuary, and most parts of the lakebed are exposed here at the mouth of the inlet during the low tide. Compared to Ashtamudi lake, the Kallada river is deeper, and the riverbed area is recorded as the average depth is greater than 13 m. The deepest part of 14.9 m is recorded near Kunnathoor bridge, which is 12 km upstream from Munroe Island. Except for a few spots of hard (resistant) rocks, the river fairly and consistently follows a higher depth throughout its course.Figure 5Bathymetric profile of Ashtamudi lake and adjoining Kallada river (Figure was generated by Arc GIS 10.6).Full size imageOnce the Kallada river supplied very fertile alluvium during its flooding seasons (monsoon/rainy season), and most of this alluvium is deposited in the floodplains of the Munroe Island and the Ashtamudi lake. With a vast river catchment area from elevated lands of Western Ghats and a shorter course of 121 km33,55 and a higher elevation gradient of 12.6 m/km56, the Kallada river has a higher transporting capacity. The eroded surface and mined river/lakebeds at lower courses were replaced by the sediment load supplied by the Kallada river during each flood season until dam construction. During the focus group discussions with residents of the Island, they had described that they were crossing the Kallada river on foot in the 1990s or even earlier during the dry seasons. The construction of the Thenmala reservoir dam in 1980s across the river drastically choked the sediment supply of the Kallada river. In addition, excessive commercial sand mining without any regulation from the riverbeds of Kallada and Ashtamudi waterbodies accelerated the deepening of waterbodies. It increased the erosion of surface and subsurface soils through fluvial and hydraulic action. This, in turn, drastically reduced the deposition of fertile alluvium over the low-lying Munroe Island. The current bathymetry shows that the river channel has deepened its course to 14 m compared to 5–6 m of 1980s. When comparing the bathymetric data of 200127, it is interesting to note that no considerable changes occurred in the bathymetry of Ashtamudi lake over the last two decades.Dams indeed alter aquatic ecology and river hydrology, upstream and downstream, affecting water quality, quantity, breeding grounds and habitation22. The other significant impact of the damming of the Kallada river is the saline water intrusion towards upstream of Ashtamudi lake and the Kallada river. The freshwater discharge is regulated after the construction of the Thenmala reservoir, and the water is being diverted to the reservoir and associated canals. There is a decline in sedimentation over the floodplains and catchment area as a result of the increased tidal effects and associated running water dynamics, which may accelerate the erosion trend of the nearby places.Lithological characterization of the Munroe IslandThe Munroe Island is a riverine delta formation by the Kallada river at the conjunction of river and backwater systems. To understand the micro-geomorphological processes of the study area, the near-surface geology of the Munroe Island had been studied in detail with the help of resistivity meter surveys and borehole datalogs from different locations. As per the current resistivity survey, it is evident that the Munroe Island is formed by recent unconsolidated loose sediments more than 120 m thick succession below ground level (Figs. 6 and 7). The electrical resistivity tomography of identified locations within the deltaic region shows a meagre resistance value to its maximum penetration (Fig. 6), which proves that the sedimentary column with intercalations of sand and carbonaceous clays of varying thickness extends to a depth of 120 m, in turn indicating the process of enormous sedimentation happened during the recent geological period. Loose wet soils of saline nature records a lower resistance value for an electric circuit. The layers formed in the diagram (Fig. 6) represent the seasonal deposition of unconsolidated soils as thin sequence. The Mulachanthara station of the resistivity meter tomography, which is situated at a more stable location of the Island, has a higher resistivity value than the West Pattamthuruth location, which is located at the exact alluvial flood plain.Figure 6Electrical resistivity profiles of Munroe Island.Full size imageFigure 7Geomorphological map showing litho-log of north (Kannamkadu); middle (Konnayil Kadavu); and south (Perumon bridge) locations of Munroe Island (borehole data source: PWD, Govt of Kerala) (Software used: Arc GIS 10.6).Full size imageThe Public Works Department (PWD), Kerala State carried out soil profile studies through Soil Penetrating Test (SPT) borehole drilling method as part of constructing bridges at three different locations up to a depth of 62 m, i.e., one across the Kallada river (north side)57, one across Ashtamudi lake in southern Munroe Island58 and one at the central part of Munroe Island (across a canal)59 (Fig. 7). The hard rock is found only on the southern side of the lake at a depth of 45 m. The litho-log shows that unconsolidated loose sediments of significantly higher thickness occur in the entire Munroe Island (Fig. 7). Anidas Khan et al.60 studied the shear strength and compressibility characteristics of Munroe Island’s soil for two different locations with disturbed and undisturbed samples. They classified the soil of Mundrothuruth into medium compressibility clay (CI) and high compressibility clay (CH) with natural moisture contents of 44.5% and 74%, respectively. The unconfined compressive strengths of the undisturbed and remolded samples for the first location are 34.5 kN/m2 and 22.1 kN/m2, respectively, while they are 13 kN/m2 and 9 kN/m2 respectively for the second location60. Such compressive strength indicates that the soils of Munroe Island are soft or very soft in nature.Land degradation: a morphological analysisTo decrease the impact of the monsoon floods and to distribute the alluvium to the southern part of the island, Canol Munroe, the then Diwan of the Thiruvithamkoor Dynasty, made an artificial man-made canal during the 1820s connecting the Kallada river with the eastern extension of Ashtamudi lake, and this river is known as “Puthanar” (meaning a new river). During the last few decades, (after 1980s) the estuarine island ecosystem of Munroe Island has faced several structural deformities. The natural sedimentation and flooding happening in the Islands were very limited and hence, the normal events happened over the past several decades disturbed and significantly affected the land neutrality. These islands, once known as the region’s rice bowl, now devoid of any paddy cultivation mainly because of the increased soil salinity. According to the Cadastral map prepared by the revenue department (1960s) there were many paddy fields, locally named as Mathirampalli Vayal (Vayal is the local name for paddy field), Thekke Kothapppalam Vayal, Mattil Vayal, Kottuvayal, pallaykattu Vayal, Konnayil Vayal, Vadakke Kundara Vayal, Thachan Vayal, Thekke Kundara Vayal, Kizhakke Oveli Vayal, Thekke Oveli Vayal, Odiyil Vettukattu Vayal, Nedumala Vayal, Madathil Vayal, Karichal Vayal, Moonumukkil Vayal, Arupara Vayal, Kaniyampalli Vayal, Manakkadavu Vayal, Panampu Vayal, Pattamthuruth Vayal etc. The recent satellite images shows that no paddy cultivation exist now, which is further confirmed by the field observations conducted through our study. The annual report published by Gramapanchayat39 indicate that the paddy field of region was reduced from 227 to 8 acres (from 1950 to 1995) and now about in 2 acres only (2018). Most of the paddy fields of northern and northwestern regions are severely affected by land degradation due to erosion, saline water intrusion and flooding and are entirely or partially buried under the backwater system. Figure 8 depicts the morphological degradation of the severely affected areas of Munroe Island from 1989 to 2021 through different satellite images. Some paddy fields are converted into filtration ponds to take the benefit of frequent tidal flooding. The coconut plantations were later introduced in place of paddy fields, and they eventually replaced the paddy fields. However, during the last decades, it has been observed that these coconut plantations are also under threat mainly because of degradation of the soil fertility, which directly bears the quality and quantity of production (Fig. 9).Figure 8Morphological changes in the study area from the satellite images (a) 1989 (aerial photograph); (b) 2000 (Landsat); (c) 2011 (World View—II); (d) 2021 (Sentinel) (the modified maps of (a) is obtained from National remote Sensing Centre (NRSC), Hyderabad, (b) is downloaded from https://earthexplorer.usgs.gov/ (c) is obtained from Digital Globe through NRSC and (d) is downloaded from https://scihub.copernicus.eu/. Figures were generated using Arc GIS 10.6).Full size imageFigure 9Threatened coconut plantations indicating the low productive regime. Photographs taken by Rafeeque MK.Full size imageOver the study area the most affected alluvial plain of the Peringalam and Cheriyakadavu island are taken separately to study the morphological changes over the decades. This area is named Puthan Yekkalpuram (which means new alluvium land), and the north side of the Kallada river (the northward extension in the Mundrothuruth GP) is demarcated as old alluvium land (Pazhaya Yekkalpuram) as per the revenue department’s cadastral map. The study shows that total 38.73 acres of land has lost from the Peringalam and Cheriyakadavu Islands during the last 32 years, which is equivalent to 11.78% and 46.95% of the total geographical area of the Peringalam and Cheriyakadavu Islands, respectively. The land degradation details over the last three decades are given in the Table 2. Many other locations, such as Nenmeni and West Pattamthuruth, are also severely affected by land degradation. However, these areas are landlocked and less affected by running water or floods. Hence, the land degradation experienced is the settling of the topsoil and subsidence of structures such as houses and bridges. The sinking of basements of many houses and even the subsidence of railway platforms are well observed during field visits, indicating the alarming land degradation issues (Figs. 1 and 10) to be addressed its deserving importance. There are also clear indications of the gradual formation of new waterlogged areas in the islands, which may further deteriorate and forms the part of the backwater system which eventually affects total land area of the Munroe Island.Table 2 Land degradation of Peringalam and Cheriyakadavu region for the past 32 years.Full size tableFigure 10Various environmental degradations in Munroe Island. Photographs taken by Rafeeque MK.Full size imageThe island population also shows a negative growth over the years. According to the census report of 201138, the total population of Gramapanchayat has decreased to 9440 person/km2 in 2011 from 10,013 person/km2 of 2001 and 10,010 person/km2 of 1991 census reports. Frequent flooding (especially tidal flooding), the lack of drinking water, and migration in search of a better livelihood are the main reasons for the observed population reduction as revealed through the survey. The high intrusion of saline water into the cultivated land through tidal flooding and the lack of flushing of surface saline soils by monsoon floods (freshwater) decreased agricultural productivity of the area, and hence, now people are more dependent on fishing and backwater activities for their livelihood. Lack of proper transportation to the nearby markets limits their fishing activities to a daily subsistence level. Due to the flooding caused by subsidence/tidal surges and land degradation during the last few decades, more than 500 households have vacated their houses38,39.Tidal Flooding and Estuarine ProcessesIn Mundrothuruth, the major environmental degradation problems where occurring due to tidal flooding and saline water intrusion into the freshwater ecosystem. Mathew et al. studied the tidal and current mechanisms of the Ashtamudi backwater in 200161. They reported that the Kallada river plays a vital role in determining the eastern lake’s circulation pattern. In addition, the increased discharge from the north Chavara canal and the south Kollam canal also influences the local circulation of the Ashtamudi backwater. The current velocity reaches up to 100 cm/s at the estuary entrance, but it rapidly diminishes in the eastern parts, where the speed is generally less than 30 cm/s. One of the critical observations made during the field study, which corroborates with the acquaintance of local people as well, is that the flooding on Munroe Island is not related to the spring tide of the open ocean. The disappearance of the semidiurnal tide in the central lakes occurs due to frictional resistance and the time lags for the tide to travel across the estuary61. At the shorter semidiurnal period of approximately 12 h, the tide is more dissipated than the more extended constituents of 24-h duration. The survey conducted with the island inhabitants also reiterates these views.As per the experience of local inhabitants, tidal flooding in Munroe Island was not frequent in earlier times. The comparison of the bathymetry data collected during 200058 and 2017 (Fig. 5) in and around the regions of Munro Islands shows that there is not much change in bathymetry during the period. Hence, changes in basin geometry are not having a significant role in tidal dynamics in imparting the variations as observed. In addition to the bathymetric survey, the data on tide measurements at four locations corresponding to three seasons were also collected. The tide data measured during the pre-monsoon period is shown in Fig. 11a. The figure shows that the tidal range in the inland area is almost the same even during the spring and neap tides. As discussed earlier, the tidal flooding in Munro Island is not related to spring tide in the ocean, and there may be the influence of specific complicated dynamics in the basin for this flooding that needs to be studied more profoundly. Further the data pertaining to tidal dynamics were inadequate; we established three tide gauges in selected locations in and around Munro Island. From the analysis of tide gauge data, it is found that the signature of anomalous variability in water column height, which is not at all linked to the tidal dynamics.Figure 11(a) Salinity variation of bottom water at selected locations in Kallada river during monsoon and post monsoon. (b) Observed tide during pre-monsoon months.Full size imageThe water quality analysis for three time periods, during the year of the cyclonic storm, Okhi (2017), was conducted to understand river run-up impact on salinity in and around Munroe Island (Fig. 11). The riverbed is lowered below the baseline of erosion, and dense saline water is trapped in the deeps during high tide. This has been confirmed during the bathymetric survey of the Kallada river and Ashtamudi backwaters, which showed a significant increase in water depth, particularly within the river channel. The high-density saline water is trapped in the basins and trenches created in the river channel due to uncontrolled sand mining, which leads to the degradation of the quality of sediments and groundwater in the region. Nevertheless, the samples collected immediately after Okhi (when the dam’s shutter was opened due to heavy rainfall in the catchment area) show that the high runoff replaced the trapped saline water with fresh water. After ten days of the first sampling, the water became saline nature after the closure of the dam’s shutter. This proves that because of dam construction, the river runoff in the Kallada river was reduced significantly, and extensive human interactions especially sand mining activities increased the riverbed deepening and formation of pools beyond the base level of running water.Conservations and management strategiesConsidering the facts discussed above, the Munroe Island may continue to be badly affected unless suitable sustainable management strategies are not evolved. Construction and associated activities, such as the damming of reservoirs, sand mining and landfilling, are indispensable for any nation’s economic and social development. United Nations’s member states have formulated 17-point Sustainable Developmental Goals (SDGs) to better the world sustainably. Local and national governments pertaining to the Munroe Island need to develop a sustainable management plan to protect this Ramsar-listed wetland. The environmental issues of Mundrothuruth can be controlled, and land degradation may be monitored through a well-drafted working plan. All aspects of earth and social sciences may be integrated to draft such a management plan of reverse landscaping. The reverse landscaping (i.e., recalling the degrading landscape to its geomorphic isostatic state) method is a must-considered sustainable solution for land degradation and other environmental issues.The deep courses of Kallada river must be upwarped through a well-planned artificial sedimentation to eradicate the saline banks of deep basins. The sediments deposited in the Thenmala reservoir and the sediments removed through the digging of boat channels may be utilized in a periodic monitoring method. Sand mining from Ashtamudi lake and the Kallada river may be strictly controlled, and the minimum freshwater flow should be ensured. The construction methods practiced in Mundrothuruth are outdated and technically nonexistent. Well-studied engineering methods suitable for an environmentally fragile area must be implemented with a proper understanding of the soil characteristics, such as shear strength and compressibility rate, and hydrodynamics, such as tidal and fluvial actions. Soil fertility must be increased by supplying additional fertile soil and freshwater, at least for a minimum period. The inhabitants’ socioeconomic well-being is strengthened by advancing technology and providing easy access to the market and other social amenities. More

Establishing size-weight relationships for heavily targeted coral species is an important first step towards informing sustainable harvest limits19. Placing coral harvests into an ecological context is a core requirement for implementing a defensible stock assessment strategy, and this need is particularly critical given escalating disturbances and widespread reports of coral loss7,17,25. Using these relationships, managers can now easily sample and calculate biomass per unit area. It is important to point out that all sites sampled in our study represent fished locations, and there is no information available to test whether standing biomass has declined due to sustained coral harvesting at these locations. While these data may now provide a critical baseline for assessing the future effects of ongoing fishing, it is also important to sample at comparable locations where fishing is not permitted or has not occurred (where possible), to test for potential effects of recent and historical harvesting.Biomass per unit area data presented herein highlights the highly patchy abundance and biomass of targeted coral species14, which is evident based on the often vastly different mean and median values (Table 2). Examining biomass per unit area estimates for C. jardinei for example, which returned some of the highest biomass estimates, the 33.75 kg·m−2 maximum estimate from a transect stands as an extreme outlier, with 12 of the 16 other transects being below 0.2 kg·m−2. This indicates the challenges of managing species that occur in patchily distributed concentrations, particularly in a management area the size of the QCF. It is also important to note, these estimates are generated only on transects where the target species occurred, and therefore, should technically not be considered as an overall estimate of standing biomass. While the estimation of size-weight relationships is a step towards a standing biomass estimate, many challenges remain in terms of sampling or reliably predicting the occurrence of these patchily distributed species. Bruckner et al.14 attempted to overcome this management challenge in a major coral fishery region of Indonesia by categorising and sampling corals (in terms of coral numbers) in defined habitat types, and then extrapolating to estimated habitat area based on visual surveys and available data. This approach, utilising size-weight relationship derived biomass per unit area estimates (instead of coral numbers), may be a viable method for the QCF, however much more information is needed to understand the habitat associations (e.g., nearshore to offshore), and environmental gradients that influence the size and abundance of individual corals. Fundamentally, it is also clear that much more data is required to effectively assess the standing biomass of aquarium corals in the very large area of operation available to Australian coral fisheries.These corals are found in a range of environments, and it is important to consider available information on life history if attempting to use coral size-weight relationships to inform management strategies via standing biomass estimation. All corals in this study can be found as free living corals (at least post-settlement) in soft-sediment, inter-reefal habitats, from which they are typically harvested by commercial collectors19. However, only four of the 6 species are colonial (C. jardinei, D. axifuga, E. glabrescens, M. lordhowensis) while the remaining two species (H. cf. australis and T. geoffroyi) are more typically monostomatous or solitary. As indicated in previous work24, if larger colonial corals were to be fragmented during harvesting instead of removed entirely, fishery impacts would likely be lessened24. Given the power relationship between coral maximum diameter and weight, larger corals contribute disproportionately to the total available biomass of each species in a given area. The potential environmental benefit of leaving larger colonies (at least partially) intact is not limited to impacts on standing biomass, as this practice would likely be demographically beneficial given the greater reproductive potential (i.e., fecundity) of larger colonies, which also do not need to overcome barriers to replenishment of populations associated with new recruits (i.e., high mortality during and post-settlement26). This conclusion was drawn largely from data on branching taxa (e.g., Acropora), which are relatively resilient to fragmentation and commonly undergo fragmentation as a result of natural processes27,28,29. D. axifuga can be considered to exhibit a relatively similar branching growth form, however, the growth form of E. glabrescens and C. jardinei changes with size, moving from small discrete polyps to large phaceloid and flabello-meandroid colonies, respectively19. While larger colonies of E. glabrescens and C. jardinei may be relatively resilient to harvesting via fragmentation, the same may not be true for smaller colonies, or species with massive growth forms such as M. lordhowensis. Typically, for each species, the average reported weight was quite low, coinciding with the lower end of the sampled maximum diameter range. For colonial species, the harvested smaller maximum diameters (if fragments) are ideal from an ecological perspective as this will have the least impact possible on standing biomass, and may also leave a potentially mature breeding colony intact. Ultimately, in light of these considerations, the development of uniform and standardised industry-wide harvest guidelines to balance economic and ecological outcomes may be necessary. The development of these guidelines would require consultation with commercial harvesters, as well as considerable additional work in measuring ecological impacts and better understanding the cost of these impacts from an economic perspective. Conversely, if whole colonies are collected, which is necessarily the case for solitary species such as H. cf. australis and T. geoffroyi (and potentially smaller colonies of other species such as E. glabrescens and C. jardinei); smaller colonies may be collected before they reach sexual maturity, hindering their ability to contribute to population replenishment. Therefore, collection of small fragments should be encouraged for colonial species; while for monostomatous species where this is not possible, introduction of a minimum harvest size based on sexual maturity should be considered.Additionally, the need for further consideration of the selectivity of ornamental coral harvest fisheries3,4,30 when assessing standing biomass is evident. Due to various desirable traits, the majority of available biomass may not be targeted by collectors. As emphasised in this study, the focus on smaller corals is indicative of the trend towards collection of most of these species at the lower portion of their size range, at least compared to some of the maximum sizes recorded on transects (e.g., see Tables 1 and 2, section b). However, it is also important to consider that transects were conducted in areas subject to commercial collection and are likely to skew results and prevent clear conclusions relating to size selectivity. Sampling of unfished populations (i.e., any residing outside of permitted fishing zones) and/or spatial and temporal matching of catch data and transect data across a larger sample of operators will be required to properly address industry size selectivity trends. For instance, only 17.5% of C. jardinei corals measured on transects fell within the diameter range represented by data obtained from collectors, with 81.9% of corals measured on transects exceeding this range. If it is viable to collect fragments from larger colonies (which does appear to be the case for some corals such as C. jardinei), then a larger proportion of standing biomass outside of this size range could be targeted by fishers. As an additional consideration, only desirable colour morphs of these corals will be harvested, and due to lack of appropriate data, the prevalence of these morphs remains unclear. H. cf. australis and M. lordhowensis for example often occur in brown colour morphs, which are far less popular in markets where certain aesthetic qualities (e.g., specific, eye-catching colours or combinations of colours) are desired, such as the ornamental aquarium industry. Even without delving into further considerations such as heritability of phenotypic traits, management conclusions drawn from standing biomass estimates may be ineffective in the absence of efforts to account for selectivity in this fishery.The relationship between size and weight was found to differ between all corals, with the exception of C. jardinei and E. glabrescens. There can be some moderate similarity in skeletal structure between these two species, particularly between small colonies, reflecting the similar maximum diameter range of sampling in the current study. Subsequently, inherent physiological constraints may be imposed on corals that prevent the maintenance of growth rates between corals of smaller and larger sizes, for example, as the surface area to volume ratio declines with growth31. In the current study, all corals, with the exception of C. jardinei, showed evidence of allometric growth, as exhibited by an estimated exponent value different to 3. Sample size for C. jardinei was greatly limited, as this species typically forms extensive beds, and are rarely brought to facilities as whole colonies. Therefore, the lack of evidence for allometric growth may reflect higher error for the species coefficient parameter due to the comparatively small sample size for this species. This suggests that mass would not increase consistently with changes in colony size in 3 dimensions31, which seems likely considering the change in exhibited form described for E. glabrescens and C. jardinei previously. In the current context, this indicates that the estimated ‘a’ and ‘b’ constants are likely to vary as the sample range increases, reflecting the changes in the size-weight relationship between smaller and larger samples of these species. Therefore, ideally, these models should incorporate data that reflect the maximum diameter range of the species in the region of application to allow increased accuracy of biomass estimation. To achieve this will require additional fishery-independent sampling, as large colonies are rarely collected whole, though may be collected as fragments depending on the species. Sampling may be challenging for some species given the difficulty of physically collecting and replacing large whole colonies, particularly for inter-reefal species such as M. lordhowensis, which can occur in deep, soft sediment habitat, subject to strong currents. Importantly, obtaining ex situ or in situ growth rate data should be considered a priority for the management of heavily targeted species. This data is likely to be another necessary component (in conjunction with size-weight relationships) of any stock assessment model developed for LPS corals, and may also eliminate the need to collect large sample colonies to improve estimated size-weight relationships.The disproportionate focus on smaller corals (i.e., corals in the current study averaged between 4.28 and 11.48 cm in maximum diameter) is likely to lead to an underestimation of weight in corals at greater diameters when used as inputs for size-weight models. This may explain the apparent minor underestimation observed in some species (e.g., M. micromussa, T. geoffroyi). In the current context, this represents an added level of conservatism with estimates obtained from these equations. While the relationship between size and weight was particularly strong for some species, (mainly D. axifuga and T. geoffroyi), for other species, such as M. lordhowensis, growth curves tended towards underestimation at larger diameter values. As the mass of a coral is reflective of the amount of carbonate skeleton that has been deposited32, the coral skeleton may increase disproportionately to coral diameter if or when corals start growing vertically. For example, in massive corals such as M. lordhowensis, vertical growth (i.e., skeletal thickening) is often very negligible among smaller colonies, with thickening of the coral skeleton only becoming apparent once the coral has reached a threshold size in terms of horizontal planar area. Additional fisheries-independent sampling outside of the relatively narrow size range of harvested colonies will be required to address this source of error in future applications. Ecological context in the form of fishery independent data on stock size and structure is essential for effective management, especially in ensuring that exploitation levels are sustainable and appropriate limits are in place. Coral harvest fisheries offer managers an ecologically and biologically unique challenge, as the implementation of standard fisheries management techniques and frameworks is hampered by their coloniality and unique biology, as well as a general lack of relevant data for assessing standing biomass and population turnover, not to mention the evolving taxonomy of scleractinian corals33. Similarly, fishery-related management challenges such as the extreme selectivity in terms of targeted size-ranges and colour-morphs, plus the potentially vast difference in the impact of various collection strategies (i.e., whole colony collection vs fragmentation during collection) also complicates the application of typical fisheries stock assessment frameworks. The relationships and equations established in the current work offer an important first step for coral fisheries globally by laying the groundwork for a defensible, ecologically sound management strategy through estimation of standing biomass, thus bridging the gap between weight-based quotas and potential environmental impacts of ongoing harvesting. It is important to note that the species selected for the current work do not represent the extent of heavily targeted LPS corals. For example, Fimbriaphyllia ancora (Veron & Pichon, 1980), Fimbriaphyllia paraancora (Veron, 1990), Cycloseris cyclolites (Lamark, 1815), and Acanthophyllia deshayesiana (Michelin, 1850) are examples of other heavily targeted corals of potential environmental concern19, and management would also benefit from the estimation of size-weight relationships for these species. Moving forward, the next challenge for the coral harvest fisheries will be to comprehensively document and track the standing biomass of heavily targeted and highly vulnerable coral stocks, explicitly accounting for fisheries effects and also non-fisheries threats, especially global climate change. More