Study area<\/p>\n
In Peru, the study was carried out along an elevational gradient from 245\u2009masl to the tree line at 3,588\u2009masl in the Andes of south-east Peru (Kos\u00f1ipata valley), with continuous and mostly undisturbed wet rainforest\/cloud forest. The climatic gradient has three seasonal periods with a wet season (November to March), a dry season (May to July) and austral spring (September and October)39<\/a>. Mean annual temperatures range from 24.3\u2009\u00b0C in the lowlands to 6.7\u2009\u00b0C at 3,600\u2009masl (Fig. 1<\/a>). Mean annual precipitation levels are high with >1,500\u2009mm per year along the whole gradient, peaking at around 1,500\u2009masl with ~5,000\u2009mm (ref. 40<\/a>). Research was conducted on 26 study plots of approximately 100\u2009m\u2009\u00d7\u2009100\u2009m in ~250\u2009m elevation intervals and seven field stations distributed along the gradient. Four of the 26 plots were located inside Man\u00fa National Park. Nine of the plots matched the long-term research plots of the Andes Biodiversity and Ecosystem Research Group (ABERG) project41<\/a>.<\/p>\n In Kenya, the study was carried out along an elevational gradient from 11\u2009masl at Watamu to 3,450\u2009masl at Mount Kenya including forests, woodland, scrub and grassland in natural and semi-natural habitats. The climate is mostly semi-arid in the lowlands to humid at higher elevations42<\/a> and characterized by seasonality in precipitation. Two rainy seasons occur from March to May and from October to December, a distinct season during the boreal summer from June to September, and a pronounced dry season in January and February43<\/a>,44<\/a>. The forested parts of the Taita Hills region and Mount Kenya from ~1,300\u20132,500\u2009masl are characterized by a tropical montane forest climate with generally high humidity and constant high precipitation (>1,500\u2009mm). Mean annual temperatures range from 26.2\u2009\u00b0C at the lowest plot to 8.9\u2009\u00b0C at the highest plot (Fig. 1<\/a>). In total, 15 study plots were selected in similar elevation intervals along the elevation gradient.<\/p>\n Insect collection<\/p>\n Insects of all major orders (mainly Coleoptera, Diptera, Hymenoptera and Lepidoptera; additionally, Hemiptera and Orthoptera) were collected in the Neotropics (n\u2009=\u20094,690) in three seasons (September to December 2022, April to August 2023, and September to December 2023) with sweep nets in the understory of all study plots. In the Afrotropics, insects (n\u2009=\u20093,164) were collected in one season (March to June 2023), applying the same method. We stored live insects in Eppendorf tubes with moist tissue and sugar solution, protected from the sun, and transported them back to the field station, where thermal tolerance measurements of all collected individuals were carried out on the same day. Insect collection was focused on taxonomic breadth covering six major orders rather than on individual species. Note, that our results represent the more common insect communities. In the\u00a0Supplementary Methods<\/a>, we provide extensive robustness analyses to assess the effects of incomplete sampling or modifications to the phylogeny.<\/p>\n Measuring thermal limits<\/p>\n We measured critical thermal limits (CT) by exposing the insects in individual plastic tubes (2, 5 or 50\u2009ml, depending on body size) to decreasing (CTmin) or increasing (CTmax) temperatures using a programmable thermoblock (Eppendorf Thermostat C)7<\/a>,45<\/a>,46<\/a>. Each individual was only tested once. To avoid the effects of starvation, each tube was equipped with a piece of paper towel moistened with sugar water. For acclimatization to a common baseline, insects were first exposed to the starting temperature of 28\u2009\u00b0C for 10\u2009min. We then increased or decreased temperatures by 1\u2009\u00b0C every 2\u2009min\u2014that is, a ramping rate of 0.5\u2009\u00b0C\u2009min\u22121\u2014following standardized methods14<\/a>. After each 2-min interval, we checked all insects for mobility. The temperature at the point of lost mobility, even after tapping or gently shaking the tube, was noted as the upper or lower thermal limit (CTmax or CTmin)14<\/a>,15<\/a>. After the test, insects were stored in 96% ethanol for later genetic barcoding and specimen documentation. Since we did not directly measure the body temperature of the insects, CTmin and CTmax values determined with the above-mentioned protocol may not perfectly reflect body temperatures. Nevertheless, owing to the small size and their high surface-area-to-volume ratio, insects equilibrate to environmental temperatures quickly, typically within seconds6<\/a>.<\/p>\n All thermal limit data were originally stored as csv or xlsx files (Microsoft Excel v.2502) and imported into R v.4.347<\/a>. Since four different observers measured CTs in the field, we first tested for potential observer bias by evaluating the CTmax of lab-reared ants of one colony. Each observer measured CTmax of ten individual workers of the colony under the same conditions. We found no significant difference in mean CTmax across any of the observers (ANOVA, F3,36\u2009=\u20091.576, P\u2009=\u20090.212).<\/p>\n For investigating CT patterns along elevation gradients, we applied generalized additive models with elevation as a smooth term explanatory variable. We restricted the smooth term parameter k to 5 to avoid overfitting48<\/a>,49<\/a>. Thermal tolerance ranges were calculated by subtracting mean CTmin from mean CTmax for each study plot. We checked the trend of thermal ranges along elevation with an ordinary linear model. We additionally applied the same models for each major insect order. Thermal safety margins were calculated for both datasets by subtracting plot-specific variables of mean annual temperature and mean daily maximum air temperature of the warmest month from the CTmax values. Both temperature data variables were derived from the CHELSA data base (bio1 and bio5 of the BIOCLIM+ dataset)50<\/a>.<\/p>\n As a test for plastic responses, a subset of insects was first exposed to heat (40\u2009\u00b0C, n\u2009=\u2009777 insects) or cold (14\u2009\u00b0C, only in the neotropics, n\u2009=\u2009490 insects) shock for 10\u2009min prior testing, replacing the 28\u2009\u00b0C acclimatization period. The sublethal temperatures for the heat shock were identified in pre-experiments. In the Neotropics up to an elevation of 2,700\u2009masl, insects generally survived a 10\u2009min shock of 40\u2009\u00b0C; at higher elevations 35\u2009\u00b0C was chosen as shock temperature, because a large proportion of the pilot test insects did not survive a 40\u2009\u00b0C shock. After the shock, the standard protocol continued from 28\u2009\u00b0C with temperature changes in 2-min intervals.<\/p>\n Plastic capacities were analysed by comparing CT measurements between individuals that were exposed or not exposed to a heat or cold shock before CT measurements. We added heat shock or cold shock (yes\/no) as a variable for CTmax or CTmin, respectively, and calculated the means for both groups. Data were additionally grouped by elevation level (lowlands: <600\u2009masl, mid: \u2265600\u2009masl and <1,200\u2009masl, high \u22651,200\u2009masl), because we expected heat shock effects to vary depending on elevation. We hypothesized that a thermal shock would increase tolerance\u2014that is, have a positive effect. Therefore, for CTmax, effect calculations were done by subtracting the mean thermal tolerance of the control group from the mean thermal tolerance of the group exposed to a heat shock (\u2018shock yes\u2019\u2009\u2013\u2009\u2018shock no\u2019). For CTmin, the calculation was done in reverse (\u2018shock no\u2019\u2009\u2013\u2009\u2018shock yes\u2019), so that the effect direction remains the same, with a positive effect indicating an increased tolerance level.<\/p>\n Insect identification and phylogeny<\/p>\n All insects were morphologically sorted into orders and, where possible, families, and delimited into species-like units based on individual DNA barcoding. For this, we sampled tissue from all specimens that were tested for thermal tolerance and positioned them in a well filled with 30\u2009\u03bcl of absolute ethanol (99.9%) in a 96-well microplate. Sequencing preparation and conduction was done at the Canadian Centre for DNA Barcoding (CCDB). Libraries were created for DNA barcoding (standard 658\u2009bp COI (mitochondrial cytochrome c oxidase subunit I) barcoding) using SMRT sequencing technology on a PacBio Sequel IIe51<\/a>. All genetic sequences were uploaded together with specimen photographs and sampling information to BOLD (database code DS-A2TP; https:\/\/www.boldsystems.org\/<\/a>). The COI sequences had an average length of 557\u2009bp with a minimum overlap of 303\u2009bp (74% of all COI sequences were larger than 500\u2009bp and 50% larger than 600\u2009bp).<\/p>\n We applied the BOLD sequence cluster tool which uses the refined single linkage (RESL) algorithm with a pairwise distance model to sort all insects into unique operational taxonomic units (OTU), which we used as species-like units (called \u2018species\u2019 in the text)52<\/a>. Contaminants, records with stop codons and sequences with a length below 300 base pairs were excluded. In total, 4,300 barcoded individuals were included in the final phylogenetic analyses of which 2,330 were unique OTUs, resulting in an average of 1.8 individuals per OTU. The CTmax data included 2,246 measurements with 1.9 individuals per OTU, the CTmin data included 1,849 total observations with a mean of 2 individuals per OTU. We used the BOLD ID engine (minimum of 80% similarity to databased sequences) to allocate all sampled insects to family level53<\/a>.<\/p>\n We used R to create a phylogeny based on a family-level backbone tree and DNA sequences20<\/a>. For this, we downloaded an insect family-level backbone tree54<\/a>. For each family, a separate subtree was calculated first, using the AlignSeqs function from DECIPHER55<\/a> and a maximum likelihood model (TreeLine function). Statistical support for each branch is provided by aBayes values (see data repository). If a family consisted of only one OTU, a \u2018tree\u2019 with a single branch of fixed length was constructed. Next, all subtrees were added to the ultrametric backbone tree using the bind.tree function from the ape package56<\/a>. Branch lengths were calibrated with the bladj function from phylocom57<\/a>, that estimates node ages based on fossil calibration points while unknown ages are evenly distributed between known nodes54<\/a>. For clearer visualization (Fig. 1h<\/a>), terminal tip heights were equalized with the forceEqualTipHeights function from ips58<\/a>. Note that while this gives the appearance of an ultrametric tree, the branch lengths do not represent absolute time, but the result is a phylogeny with relative branch lengths. All scripts and details of phylogeny construction can be found in20<\/a>.<\/p>\n We reconstructed ancestral trait values of thermal tolerances using the fastAnc function from phytools59<\/a>, plotted them on the phylogenetic tree with ggtree60<\/a>, and calculated a phylogenetic correlogram using phylosignal to test for significance of the phylogenetic signal\u2014that is, if trait values of related OTU are more similar (or dissimilar) than expected by chance21<\/a>. The phylogenetic parameters Pagel\u2019s lambda and Blomberg\u2019s K were tested with the phylosig function and 10,000 randomizations. To disentangle the effects of adaptation or acclimatization to a local climate from those of phylogenetic relatedness, we applied a phylogenetic regression using the phylolm package, and function of the same name, with elevation as predictor and phylogenetic relationships entering a covariance matrix of the model either based on the assumption of a Brownian motion and an Ornstein\u2013Uhlenbeck model of trait evolution61<\/a>. We compare the Ornstein\u2013Uhlenbeck and Brownian motion models since they represent two mechanistically contrasting hypotheses about trait evolution and are most commonly used in literature9<\/a>. Under a Brownian motion model, traits are assumed to have evolved under random evolutionary drift, under an Ornstein\u2013Uhlenbeck model with stabilizing selection towards an optimum. We calculated the evolutionary optimum for CTmax under an Ornstein\u2013Uhlenbeck model\u2014that is, the parameter \u03b8, using the fitContinuous function of the geiger package62<\/a>. We compared model performance using the Akaike information criterion63<\/a>. Variance partitioning to extract partial r2 values was conducted on these models using the phylolm.hp package64<\/a>.<\/p>\n In order to verify the dependence of results concerning elevational trends in CT, on the reconstruction of ancestral traits, the estimation of a phylogenetic signal and of an upper boundary of CTmax (Ornstein\u2013Uhlenbeck model), we conducted an extensive set of robustness analyses considering various modifications of the phylogenetic tree and various data subsets, which is described in the\u00a0Supplementary Information<\/a>, with results presented in Extended Data Fig. 3<\/a> and Supplementary Table 1<\/a>.<\/p>\n Protein stability prediction<\/p>\n We predicted the thermal stability of proteins applying the deep learning model DeepSTABp26<\/a>. It uses a transformer-based protein language model to extract sequence embeddings, which are analysed with advanced deep learning techniques for large-scale protein stability predictions26<\/a>. To cover proteins from species across all insect orders, we downloaded all available genomes (protein format) from 677 insect species from InsectBase 2.065<\/a>. This data covered 20 orders, 158 families and 457 genera. Next, in R, we imported all genomic data translated into amino acid FASTA files and then randomly selected 1,000 proteins per species, a trade-off to receive as many overlapping proteins between species as possible while keeping the following analyses at a feasible computation time. We set up a local version of DeepSTABp, to predict melting points for this large number of proteins. Anaconda PowerShell was used to create a conda environment and Python programming language (v.3.13) to run the script66<\/a>. In DeepSTABp, we set growth conditions to \u2018cell\u2019 and a default temperature of 22\u2009\u00b0C (ref. 26<\/a>). For testing robustness of results, we conducted analyses with data from all genomes, only high-quality genomes (BUSCO\u2009>\u200989.9, N50\u2009>\u2009300\u2009kb), and only for proteins which were covered by all insect species. Furthermore, for comparison we used data of all proteins and of the 25% proteins per species that showed the highest thermal sensitivity (n\u2009=\u200988,643 proteins from 677 species).<\/p>\n To investigate the underlying, structural mechanism of heat tolerance in insects, we predicted CTmax by protein melting temperatures Tm. Differences in Tm across orders and families were tested with mixed effect models (lme), setting species-specific protein identity as random term. For the comparison of Tm with field-measured CTmax, we calculated mean CTmax values for each family of the Afrotropical and Neotropical data. We calculated the 25% quantile of Tm for each species (that is, a proxy for an average Tm of temperature sensitive proteins) of the genomic dataset and then averaged these values among all species per family. Using a linear model (lm), we tested the relationship between Tm of the temperature sensitive proteins (average of the 25% quantile of Tm) and the mean CTmax of families. Here, we first calculated an unweighted linear model and, second, a model weighted for the number of species with genomic data in each family.<\/p>\n Climate data<\/p>\n On all Neotropical plots, one TMS-4 soil and one air temperature logger (TOMST) recorded air temperature at 15\u2009cm and 2\u2009cm above the ground, as well as soil temperature and moisture at \u22126 cm depth in 15-min intervals for approximately one year (September 2022 to December 2023)67<\/a>. Air temperature was additionally measured at 1.5\u2009m height with iButton sensors (Analog Devices) in 240\u2009min intervals. The sensors were protected from rain and sunlight using white plastic dishes (diameter 18\u2009cm)68<\/a>. Air temperatures at 1.5\u2009m (Tair at 1.5\u2009m) were highly predictive of the temperature at 15\u2009cm height (Tair at 15\u2009cm)\u2014R\u00b2\u2009=\u20090.999, Tair at 1.5\u2009m\u2009=\u20090.26673\u2009+\u20090.99845\u2009\u00d7\u2009Tair at 15\u2009cm\u2014and due to the denser sampling intervals, data of the TMS-4 sensor were used for final statistical analyses. At one plot (ID 056), the TMS logger could not be recovered due to fallen trees; therefore, the mean annual temperature from iButton data was used to predict Tair at 15 cm for this plot instead. For the East African study plots, we did not have temperature loggers available but retrieved comparable data by modelling shaded air temperatures using NicheMapR. We simulated hourly temperatures at every plot across one year (micro_global function). This climatic model takes detailed environmental variables into account, such as elevation, slope, solar radiation and absorptivity69<\/a>.<\/p>\n Additionally, for all study plots climate data (BIOCLIM+: mean annual air temperature (bio1); mean daily maximum air temperature of the warmest month (bio5); mean daily minimum air temperature of the coldest month (bio6)) were extracted from CHELSA data base. This dataset is based on temperature measurements of climate station from the years 1981\u20132010. This was done to have climate data with higher comparability between the geographic regions and to model future temperatures at the study plots. CHELSA was chosen as it is particularly suitable for modelling climate along mountain slopes due to its terrain-based downscaling of global data (https:\/\/chelsa-climate.org<\/a>). Both bio1 and bio5 climatic variables were highly correlated to estimates of the same measures based on the field-measured climatic data (Extended Data Fig. 9<\/a> and Supplementary Information<\/a>).<\/p>\n Surface temperatures for all study plots were derived from the ECOSTRESS sensor, a radiometer mounted on the International Space Station that measures surface temperatures with a resolution of 70\u2009m (ref. 29<\/a>). As ECOSTRESS measures surface temperatures at a low temporal resolution and significant parts of the data were filtered due to low data quality (for example, due to cloud cover), surface temperature measurements of study plots were enriched by additional measurements of 50 randomly selected locations within a buffer of 5\u2009km around each plot. For all additional locations the elevational level was determined using the NASA SRTMGL1 digital elevation model (v.003) with a spatial resolution of 1 arc-second (approximately 30\u2009m). We applied a filter (binary \u20180b00\u2019\u2009=\u2009best quality retrieval) to only use high-quality surface temperature data. For better visualization, we only plotted temperatures above the 50th percentile of temperature values in Fig. 3<\/a>.<\/p>\n Heat coma models<\/p>\n For better estimating the effects of current and future temperatures on tropical insects, we converted the dynamic CT values to static CT values and calculated heat coma times tcoma for present and predicted future environmental temperatures19<\/a>. The static CT model includes a thermal sensitivity coefficient z, which describes how knockdown time changes with temperature. Across different insect species, z varies between 1 and 5, with an average value of 3 (ref. 19<\/a>). To provide additional estimations of z, we tested ramping rates from 1 to 0.06\u2009\u00b0C\u2009min\u22121 with three species of leaf cutter ants in the Neotropics (Acromyrmex coronatus, Acromyrmex octospinosus and Atta colombica) as model organisms, since they could be easily found along an elevation gradient from 270\u20131,349\u2009masl. In the ant data, we found z to range from 2.16\u20135.50 along the gradient, which overlapped with reports from the literature. For the final static value calculation, we used z\u2009=\u20093, following a conservative approach. Under the assumption of higher z values, tcoma would be shorter than the values reported in Fig. 3<\/a>.<\/p>\n For the calculation of the effect of environmental temperatures on tcoma of insect communities under current climate and climate change scenarios we extracted community-level CTmax values from the lowest plot and applied it across the whole elevational gradient, a rather conservative approach, assuming that the thermal sensitivity of insects across the elevation gradient can, by species turnover, adaptation, or acclimatization, increase to values of CTmax currently found in lowland species. We calculated tcoma assuming an average CTmax (median CTmax of all lowland insects), for a 25% quartile of CTmax (termed \u2018more heat-sensitive insects\u2019) and for a 10% quantile of CTmax (\u2018most heat-sensitive insects\u2019). In the Neotropics, the 10% quantile of all insects measured in the lowland CTmax was 40\u2009\u00b0C, the 25% quantile was 41\u2009\u00b0C, and the median was 42\u2009\u00b0C. In the Afrotropics, the 10% quantile of CTmax was 43\u2009\u00b0C, 45\u2009\u00b0C (25% quantile) and 47\u2009\u00b0C (median). For future climate projections (2071\u20132100) under the three shared socio-economic pathways SSP1-2.6 (sustainable scenario), SSP3-7.0 (medium-high scenario) and SSP5-8.5 (high scenario) we used the GFDL-ESM4 climate model predictions (highest priority; https:\/\/protocol.isimip.org\/#\/ISIMIP3b\/biodiversity<\/a>) of bio5 supplied by the CHELSA database50<\/a>. We extracted the mean daily maximum air temperature of the warmest month (bio5) and calculated the anomaly based on the difference between current and future bio5 temperatures (Extended Data Fig. 8<\/a>). We then estimated the future air temperatures by adding these anomalies to the current environmental temperatures measured in the field (Neotropics) and modelled microclimate (Afrotropics), as well as for the surface temperatures (ECOSTRESS) on each plot. Using this approach, we incorporated the fine-scaled temporal variation in environmental temperatures along the Afrotropical and Neotropical elevation gradients to future climate projections. Thereby, we assume that the differences of environmental temperatures at small temporal scale to bio5 remains the same under future climate projections (details in\u00a0Supplementary Information<\/a>). We acknowledge that surface temperatures may not warm at the same rate than air temperatures. The difference between surface and air temperatures is expected to increase with increasing temperatures, particularly in open areas70<\/a>. Thus, our method provides conservative estimates of future surface temperatures. For calculating the percentage of critical temperatures (leading to heat coma within an exposure time of 8\u2009h), we related for each study plot the number of critical temperatures to the total number of temperature values (including 100% of all temperature measurements\u2014that is, including night temperatures).<\/p>\n While reported results are mainly based on the GFDL-ESM4 climate model, all analyses were additionally calculated with a multi-(climate) model ensemble. The detailed methods are described in the\u00a0Supplementary Information<\/a> and results of the multi-model-ensemble are shown in Supplementary Fig. 1<\/a>.<\/p>\n Inclusion and ethics statement<\/p>\n Fieldwork in Peru and Kenya was conducted in collaboration with local research institutions, including the Universidad Peruana Cayetano Heredia and the Museo de Historia Natural in Lima, and the University of Embu in Kenya. All necessary permits were obtained from the respective national authorities (see Acknowledgements). Local scientists were involved in study design, data collection and manuscript preparation, and are included as co-authors. Field logistics were supported by local assistants, who are acknowledged accordingly. The research supported local capacity building through training initiatives and shared data access. The collection of biodiversity data addresses locally relevant priorities and contributes to future monitoring and conservation efforts.<\/p>\n Reporting summary<\/p>\n