s2fieldsofstudy:
Biology
Environmental Science

text:
Ecological and evolutionary trends of body size in Pristimantis frogs, the world's most diverse vertebrate genus

Body size is a key organismal trait. However, the environmental and evolutionary factors that drive body size patterns at the interspecific level remain unclear. Here, we explored these relationships between phenotype-environment using neotropical frogs of Pristimantis, the world’s most diverse vertebrate genus. We analyzed: (a) whether this group follows the Rensch’s rule, a trend of sexual size dimorphism (SSD) to increase with size when males are the larger sex; (b) whether environmental constraints have influenced body size variation; and (c) how the rates of body size evolution have varied over time. Analyses were based on two information sources, the first one including body sizes of ~ 85% (495 species) of known species in the genus, and a second one incorporating molecular phylogenetic information for 257 species. Our results showed that all Pristimantis species exhibited marked SSD but did not follow Rensch’s rule. We found that the models that best explained body size in males, females, and SSD contained environmental variations in temperature, precipitation, and elevation as predictors. In turn, body size has evolved toward an optimum, with a decelerating rate of evolution differentiated between the large Pristimantis clades.

Body size is perhaps the most visible trait for most organisms. In some cases, body size correlates with different ecological and evolutionary factors [1][2][3][4][5] , which has generated a commonly repeated question: what are the ecological and evolutionary processes that drive body size limits 6,7 . This question has been approached from a few different perspectives, including studies on allometry and trends correlated with body size 8 , relationships between physiological processes and body size 9,10 , and the processes responsible for evolutionary trends in body size 11,12 . Similarly, questions related to the interspecific differences in body size, particularly related to sexual size dimorphism (SSD; i.e., the differences in size between the sexes), continue to be an essential study subject for understanding how the sexes reach different optimal sizes and which ecological and evolutionary mechanisms drive these differences 13 . Sexual size dimorphism is widespread in different animal groups and is variable even between closely related groups 14,15 . SSD varies in terms of its direction due to differences in the body size of different taxa of vertebrate and invertebrates, with SSD patterns tending to lean towards increased body size in males (e.g., some birds and mammals), but also producing inverse cases where the SSD is female-biased 16,17 . This pattern, represented by the increase in SSD with increasing body size in males, is known as Rensch's rule 18,19 .
The drivers behind the variation in body size and its disparity between males and females have been associated with ecological, life history, and evolutionary factors 20,21 . Therefore, the integration of ecological information (e.g., elevation, temperature, and precipitation) and evolutionary information (e.g., ancestral states and rates of trait evolution) can play a decisive role in understanding the forces driving the variations and evolution of body size between sexes 22,23 . For example, Rensch's rule assumes that a larger body size in males is associated with a faster evolutionary rate (speed of genetic or morphological divergence in a lineage per million years), whereas a larger body size in females is associated with a slower evolutionary rate 14,15 , with a macroevolutionary trend where male-biased SSD drives the direction of SSD between related species 24 . However, other studies have rejected Rensch's rule, for example, in studies that have indicated inverse patterns where SSD is female-biased (e.g., [25][26][27] . Evolutionary contexts with an isometry between sexes indicate that changes in body size occurred at the same rate and thus equally influenced the evolution of SSD (e.g., 26,28 , whereas evolutionary contexts with mixed SSD patterns resulted in phylogenetic lineages with species associated with a female-biased SSD and others with a male-biased SSD 24,29 . Hypotheses associated with classical ecogeographic rules that relate geographical patterns in the variation of biological traits and environmental gradients are important for understanding interspecific patterns of body sizeassociated traits at different scales and ecological and evolutionary contexts 30 . Climatic variables have been found to be a major driver of broad-scale interspecific body size patterns for different taxonomic groups, including cold-blooded animals e.g., [31][32][33][34] . In fact, two hypotheses have been recurrent related to body size in amphibians, and both have to do with the importance of surface-to-volume ratios on the homeostatic capacity of organisms through its effects on thermoregulation and hydroregulation 31 . According to the "water availability hypothesis", larger organisms (i.e., those with a lower surface-to-volume ratio) would be favored in drier environments since they have lower rates of evaporative water loss and, hence, are more resistant to desiccation 35,36 . In thermoregulating anurans, "the heat balance hypothesis" suggests that the higher heat retention capacity associated with a lower surface-to-volume ratio would benefit organisms with large body sizes in cold environments 31 . This hypothesis, analogous to the classic heat conservation explanation of Bergmann 37 for size gradients in endotherms, can be extended to ectotherms below a certain body size threshold especially in the presence of compensatory physiological, morphological and/or behavioral mechanisms to reduce heating times 32,38 . For example, heating rates are faster in dark ectotherms than in light-colored ones, which is particularly advantageous for the survival of larger species under low ambient temperature conditions 39 . Finally, the habitat availability hypothesis associates the topography-macroclimate interaction, where body size variation may be influenced by habitat availability, where reduced areas (high mountains) may present decreases in body size in contrast to lowlands 40,41 .
In the context of body size evolution, it has been shown that there is no generalized pattern at the geographic and ecological level among anuran families 42 . This indicates that the determinants of changes in optimal body size are more related to clade-specific changes at a particular geographic scale for each anuran group. Thus, using an ecological and evolutionary approach, we evaluated the patterns of interspecific variation in body size using the amphibian genus Pristimantis as model study. Pristimantis is the most diverse vertebrate genus on the planet, with 584 known species 43 . Therefore, this genus is an excellent group for studying large-scale interspecific variations associated with body size due to its wide latitudinal and altitudinal distribution, from Central America to the Andean and Amazonian areas 44 , and from sea level up to 4500 m 45 . Pristimantis originated in South America during the early Miocene mainly associated with Andean orogenic events 46 . Studies on this genus have mainly been focused on taxonomic, systematic, and biogeographical contexts. However, the patterns related to interspecific variation in body size have only been partially evaluated 47 , with the existing information on body size being found across species descriptions, natural history notes, and selected morphological studies.
We used two approaches to investigate body size patterns. First, we collected body size data for all known Pristimantis species (584 species) without phylogenetic information. Second, we included phylogenetic information on 257 species for which body size data and molecular information were available for each species. We analyzed: (a) whether the SSD would conform to Rensch's rule, where SSD should increase with body size when SSD is male-biased; (b) whether male and female body size and SSD are driven by environmental factors, about which we predicted that bioclimatic patterns would be positively associated with variation in body size; and (c) how the evolutionary rates of body size changed throughout evolutionary time with different patterns between males and females, for which we predicted a pattern of deceleration in the body size evolution rate alongside a decrease in both male and female body size.

Results
General pattern of SSD. Maximum body size varied considerably among Pristimantis, from 12.6 mm (P. xeniolum) to 50.8 mm (P. labiosus) in males and 15.3 mm (P. coronatus) to 69.3 mm (P. lymani) in females (Supplementary Table S1). Our results showed that 468 Pristimantis species evaluated in this study exhibited marked sexual dimorphism related to body size (Fig. 1A), while 27 species are monomorphic in body size (SDI < 0.1) (Supplementary Table S1). In turn, none of the species analyzed have a male-biased SSD. The species with the lowest SSD was P. jabonensis (0.003) and the highest was P. latidiscus (1.06) (Supplementary Table S1). On average, SDD was female-biased with a near-normal distribution (Fig. 1B). The linear regression between the body sizes of males and females showed that the allometric slopes were less than 1 (β = 0.83922 ± 0.04967, 95% CI 0.8027025-0.8757443; R 2 = 0.80, F 1,493 = 2038, P < 0.001; Fig. 2); this variation in SSD did not follow Rensch's rule. Phylogenetically corrected analyses also showed a positive correlation between male and female sizes (R 2 = 0.72, F 1,255 = 686.4, P < 0.0001).
Environmental predictors of body size. Our analyses (with and without phylogeny) found significant patterns that showed that SSD and body size of males and females were associated with climatic variation related to the heat balance, water availability and habitat availability hypothesis. SSD analysis conducted by selecting models based on information theory identified temperature-associated variables (BIO1: annual mean temperature) as more relevant to explain SSD patterns in Pristimantis, with a greater relative importance of 0.99 supporting heat balance hypothesis ( Fig. 3A www.nature.com/scientificreports/ porated the variables BIO1, BIO2, BIO4 (heat balance hypothesis) and altitude (habitat availability hypothesis); for SSD, the best model included the variables BIO3 (heat balance hypothesis) and BIO15 (water availability hypothesis) ( Table 2).  (Table 3).

The evolution of body size in
Evolutionary rate of body size. Bayesian analysis of macroevolutionary mixtures (BAMM) supported changes in the evolution of body size, with the rates for males and females having a marginal probability of 0.988 (> 95% of all models sampled from the post-burn-in chain). A regimen change located at node 98 (20 Mya) between clades G1, G2, and G3 occurred in males and females. Regime change was associated with a reversal in the rate of body size evolution; thus, within these clades, the rates tended to decelerate over time (Fig. 5A,B; Supplementary  Table S7). In turn, males and females showed different patterns in the rates of body size evolution, where females

Discussion
Pristimantis body size patterns exhibited marked SSD. However, the genus conforms to a female-biased inverse Rensch's rule (b < 1). This is consistent with previously documented patterns for the three amphibian orders 47,48 . This deviation of SSD towards females could be related to fecundity selection, which describes the fitness advantages that result from the selection of traits that increase the number of offspring per reproductive period [48][49][50] . However, despite knowing that all Pristimantis species have direct egg development, the clutch size of most species is unknown, although estimates for a few species have been between 30 and 38 eggs per clutch 51-53 . Other  www.nature.com/scientificreports/ explanations have attributed the variations in body size between males and females to the energy requirements associated with reproductive behavior. In males, investment is usually related to territorial defense, agonistic, and advertising behavior, whereas females invest more in body growth 49 . Accordingly, variations in body size have several ecological and evolutionary implications. For example, larger amphibians tend to have lower metabolic rates for maintaining their temperature compared to smaller amphibians 54 . However, our results suggest that environmental conditions related with temperature (heat balance hypothesis) are associated with body size differences between males and females. SSD in Pristimantis exhibits a positive  www.nature.com/scientificreports/ relationship for temperature seasonality (BIO4), which is a measure of temperature change over the year 55 , and annual mean temperature (BIO1), which is the average maximum and minimum temperature for the warmest and coldest months, respectively 55 .
Previous studies have shown that temperature seasonality is related to body size in amphibians, and this relationship would be associated with the phenomenon that shapes latitudinal gradients 48,56 . However, over 60% of Pristimantis species are distributed in Andean areas, where local temperature variations tend to be relatively stable 57 . Therefore, the annual mean temperature may be more related to SSD patterns in Pristimantis, especially because this variable presented an importance of 0.99, which indicates that the variable has appeared with high weights repeatedly among the different candidate models. However, around 18% (89) of the species are restricted to the páramo zone (> 3000 m a.s.l.), which includes more fluctuating thermal landscapes that can influence the physiological patterns relevant to body size differentiation 58 that varies between 12 and 50 mm. Our data support a slight negative relationship between annual mean temperature (Bio1) and SSD (R 2 = 0.05, F 4,84 = 2.304, P = 0.009) for species above 3000 m a.s.l.
The climatic predictors that best explain body size for each sex are precipitation seasonality (water availability hypothesis) and elevation (habitat availability hypothesis). Although mean annual precipitation is used most often to describe different ecological relationships, the dynamic of temporal variation in precipitation under different geographical settings can influence ecosystem responses and the morphological trait patterns of species 59 . The elevation is also associated with body size in Pristimantis, especially in females, possibly due to the habitat availability hypothesis, which predicts larger body sizes in lowlands relative to mountains 41 . In turn, body size has been associated with habitat variability 60 and as a predictor of brood size 61 . For example, in highaltitude areas, it has been shown that breeding seasons are shorter for some amphibian species. Consequently, the number of clutches is lower 62 , which could result in smaller body sizes. However, these aspects remain to be explored in the case of Pristimantis, for which data associated with the number of clutches are still incipient. Therefore, possible relationships between clutch size, altitude, body size, and Andean bioclimatic landscapes remain as hypotheses to be evaluated.

Evolution of Pristimantis body size.
Body size of Pristimantis species ranged from 12.6 to 50.8 mm for males and between 15.3 and 73 mm for females. A recent study by 47 reported P. nanus (not included in this study) as the smallest species, with males being only 12.4 mm. This new species is part of the P. trachyblepharis species group, which is a group of micro-endemic taxa distributed in south-central Ecuador and northern Peru 47 . Our analyses did not identify any clades with a tendency to miniaturize, as size patterns varied within each clade ( Fig. 7; Supplementary Fig. S2). This implies that the species within each clade present multiple selective pressures that have determined the wide ranges of body size in Pristimantis. The tendency toward miniaturization or gigantism in anurans is more related to clade-specific patterns associated with responses to ecological and geographic factors 42 . However, among clades with different reproductive modes, miniaturization has been shown to be the main pattern in direct-developing amphibians (the only mode in Pristimantis) 42 .
The phylogenetic signal of body size for males and females was high (Pagel's λ > 0.85), which supports the findings previously reported by Zumel 47 and those reported for amphibians in general 63 . However, for other neotropical amphibians, such as the genus Boana which has a similar distribution to Pristimantis, a low phylogenetic signal for species has been reported, having large variations in body size with random patterns of trait variation 33 . Nevertheless, the strong phylogenetic signal found for Pristimantis sexes indicates that the data fit a Brownian motion model well, indicating that closely related species tend to have similar body size 64,65 .
Evolutionary rate of body size. Although our phylogeny is incomplete, body size evolution rates for male and female Pristimantis suggested a regimen change associated with progressive rate deceleration, with a break between clades G1-G2 and G3 (Fig. 5A,B; Table 3). It is important to clarify that for estimation of diversification rates (speciation-extinction analyses), incomplete sampling of taxa may have an impact on rate estimates 66 . However, for modeling trait evolution, it has no implications if taxon sampling is random with respect to phenotypes 67 . Changes in body size evolution rates of Pristimantis may be associated with historical biogeographic patterns that have shaped its diversification patterns. Clades G1-G2 showed the highest ranges of body size evolution rates, between 5.2750 and 5.6506 for males and 9.0240-11.1643 for females (Table 3). These early clades consist mainly of species from the ecoregions of the Northwestern Andes ecoregions, which have been suggested as a possible center of origin of Pristimantis 68 . Subsequent dispersal patterns with the colonization of ecoregions such as Chocó, Amazonia, the Colombian Eastern Andes, Central Andes, the Venezuelan Andes, the Atlantic Forest, Northern South American montane coastal areas, and the Guiana Shield have shaped clades G1-G2 (Acevedo et al., unpublished). Clade G3 (Fig. 5A,B; Table 3) exhibited evolution rates with a tendency to slow down between 3.2704 and 3.0947 Ma for males and 6.9802-6.4955 Ma for females (Table 3). Clade G3 comprises species from the Northwestern Andes, the Colombian Eastern Andes, Amazonia, Central Andes, Central America, and the Guiana Shield. The deceleration pattern is explained by the trend towards an evolutionary optimum in body size that Pristimantis species have experienced, where forces acted to reduce phenotypic variation. In contrast, the initial acceleration may be the response of the ancestral lineages to some selective pressure experienced in the early stages of genus diversification.
Females had higher evolutionary rates than males (opposite of Rensch's rule). For both sexes, there was a pattern of deceleration in the rates of body size evolution, and this was associated with reduced body size (Fig. 6A,B). This finding is in contrast with patterns found in mammals, for example, where the deceleration of body size evolution rates is associated with an increase in size 69 . However, at the macroevolutionary level for different amphibian groups, the patterns of the evolutionary rates of traits associated with body size are not yet clear. Therefore, future studies must test ecophysiological hypotheses, niches, and other ecological attributes to www.nature.com/scientificreports/ obtain a general understanding of the changes in the patterns of body size evolution rates. Our results provide an approach for understanding the relationship between body size and its evolutionary rates in Pristimantis.
In conclusion, our results support the evolution of body size and SSD variation in Pristimantis, suggesting that environmental variables can drive and/or maintain the degree of body size divergence between sexes. The rates of body size evolution show a deceleration over time with a slight increase in body size. These results contribute to knowing the factors involved in body size and its evolution, helping to generate future hypotheses related to sexual and natural selection or macroevolution tendency in SSD.

Methods
Datasets. We collected data on maximum body size, represented by snout-vent length (SVL), for adult males and females of all known Pristimantis species (584) described until July 2021. Maximum size was taken as a proxy for the potential size of organisms with indeterminate growth 54,70,71 , a common standard in studies that explore interspecific size comparisons in anurans 31,32,42 . Body sizes were log-transformed (log 10 ) to reduce the asymmetry of the data since they cover a wide range of body size values. The species included in this study represented the entire known distributional range for the genus, spanning from Central America to the Central Andes of Peru and Bolivia. Measurement data were compiled from the original descriptions of each species, field data, and museum records on the review and verification of specimens (Colección Herpetológica del Instituto de Ciencias Naturales, Universidad Nacional de Colombia; Colección Herpetológica de la Universidad de Pamplona; Museo de Historia Natural de la Universidad Distrital Francisco José de Caldas; and the Museo de Historia Natural de la Universidad de los Andes). Measurements were taken using a digital caliper to the nearest 0.01 mm following 72 . The limited information in the literature regarding body size for multiple males or females per Pristimantis species did not allow us to perform additional analyses to determine the influence of intraspecific sampling and its implications on SSD estimates. www.nature.com/scientificreports/ We built two databases for developing the two analytical approaches. The first approach focused on ecological analysis without including a molecular phylogeny, and the final database included 495 species (85% of the known species of the genus) the remaining 89 species only present data for one of the sexes, therefore they were excluded (Supplementary Table S1). For the second database, we used a subset of species on which we performed additional evolutionary analyses. Although there are molecular data available in different public repositories (e.g., Genbank) for approximately 300 Pristimantis species, we only included 257 species for which there were molecular data and body size measurements were available for both sexes (Supplementary Table S2).
To estimate the SSD for each species and to evaluate Rensch's rule, we used the size dimorphism indices (SDI's) proposed by Lovich and Gibbons 73 : Divergence time estimation. The divergence times for the 257 species (Supplementary Table S3) that were included in the body size matrix were obtained from a concatenated data set for 304 Pristimantis species for six partial mitochondrial genes (12S, 16S, CYTB, COI, ND1, and ND2) and two partial nuclear genes (RAG1 and TYR). Two hundred sixty-one species were retrieved from GenBank (as of August 2020). The remaining sequences comprised 43 previously unsequenced Pristimantis species (unpublished data), of these 24 were included in the final body size matrix.
We selected the best partition scheme and the corresponding substitution models using PartitionFinder 2 74 with the Bayesian information criterion and the greedy algorithm 74 . We assumed 19 possible partitions through the concatenated data matrix (by genes and codons).
Divergence times were estimated using a relaxed Bayesian clock implemented in BEAST 1.10.4 75 , using a Yule speciation process 76 with substitution rate variation. The analysis was performed using 100 million generations that were sampled every 1,000 generations. We based our analysis on four previously published divergence times: Cross-species analyses. We used a cross-species approach in which each species represented independent data 78 . We created an extensive compilation of the latitudinal and longitudinal coordinates of the 495 species. The geographical occurrence data were obtained through different methods, such as from the Global Biodiversity Information Facility (https:// www. gbif. org), the scientific literature, a review of biological collections (as done for body size), and our own field data. Records for each species were mapped in QGIS v3.14.16 79 and individually curated to correct georeferencing errors and eliminate erroneous locations. The final database contained 9237 geographic records (Supplementary Table S4 includes geographic information for all Pristimantis species).
We selected 24 climate variables at a resolution of 1 km 2 , taken from Worldclim (https:// www. world clim. org) and ENVIREM (https:// envir em. github. io). We excluded highly correlated variables (Pearson correlation coefficient ≥ 0.60) from the analyses to minimize multicollinearity using Spearman's rank correlation performed in the corrplot R package 80 (Supplementary Fig. S1; Table S5). The resulting variables were used to test the two main hypotheses: (a) heat balance 31,81 : annual mean temperature (BIO1), mean diurnal range (BIO2), isothermality (BIO3), and temperature seasonality (BIO4); and (b) water availability 36 : annual precipitation (BIO12), precipitation seasonality (BIO15), and moisture index. Additionally, we included altitudinal information to estimate the available habitat formed by the interaction between the topography and macroclimate 41,82 . For each species, the variables were averaged based on the geographic occurrences (Supplementary Table S1). The predictors were logarithmically transformed and then scaled to have a mean of zero and unit variance using the 'standardize' package in R 83 .
Phylogenetic comparative analyses. We estimated the phylogenetic signal of the residual errors for the body sizes of males and females using the statistic lambda λ (Pagel's λ). When λ = 1, it signifies that trait relationships are equal to the correlation of species imposed by their shared evolutionary history 84 . When λ = 0, it signifies patterns of trait similarity between species that are independent of phylogeny (non-phylogenetic regression) 85,86 . The phylogenetic signal tests were implemented in R using the 'phylosig' function of the 'phytools' package 87 .
We mapped the ancestral states of male and female body size that were transformed to the log10 scale for the internal nodes of the calibrated phylogeny obtained in BEAST. We used the 'contMap' function in the R 'phytools' package 87 . This function uses the maximum probability to estimate the ancestral states, with the trait data for each species at the tips of the tree 88 .

Bayesian analysis of macroevolutionary mixtures.
We performed an analysis in BAMM v. 2.5.0 89 to estimate and quantify the heterogeneity in the evolutionary rates of body size in males and females (rate of phenotypic evolution, ß). BAMM uses reversible-jump Markov chain Monte Carlo (rjMCMC) to choose between models that vary in the number of evolutionary regimes 90 www.nature.com/scientificreports/ rior = 1.0, betaInitPrior = 0.054371, betaShiftPrior = 0.051182) and females (poissonRatePrior = 1.0, betaInitPrior = 0.027243, betaShiftPrior = 0.051182) were estimated using the BAMMtools v. 2.1.7 package 89 in R before analysis in BAMM. We used a consensus tree obtained previously in BEAST, and the chains were run for 100 million generations with a sampling frequency of 10,000. The first 10% of samples were discarded as burn-in. An effective sample size (ESS) > 200 of the log-likelihood and the number of shift events present in each sample were evaluated using the R package 'coda' 92  Ecological and evolutionary analysis. To evaluate whether SSD patterns related to the increase, decrease, or isometry of the maximum log10-transformed Pristimantis body size, we performed simple linear regressions using the 'lm' function in R 93 for the database that included all species. We also performed a phylogenetic generalized least squares (PGLS) regression including the calibrated phylogenetic tree. For both analyses, we placed females on the x-axis (dependent variable) 15 . We examined the slope of the regressions of male and female body size according to three scenarios following Fairbairn 94,95 : (a) Rensch's rule, slope greater than one (b > 1); (b) inverse of Rensch's rule, slope less than one (b < 1); and (c) isometry, slope equal to one (b = 1). Additionally, we performed a simple linear regression comparing the SSD values against the SVL log10 of males and females to evaluate the increase or decrease in SSD according to sex.
To determine the predictors that best explained the interspecific variation in body size, we performed two analyses, one with the database for the 495 Pristimantis species (Supplementary Table S1) and another one including phylogenetic information and male and female body size of 257 species (Supplementary Table S2). Analyses were performed separately for males, females, and SSD. For the first analysis (495 species), we selected models based on information theory criteria 96 . We used the R 'glmulti' package 97 , which finds the set of confidence models among all possible models. The best models were found by a genetic algorithm (GA), from which the multi-model average was derived using the 'coef ' function. The GA incorporates an immigration operator, which allows the eliminated variables to be reconsidered. Immigration increases the level of randomization and, therefore, the probability of convergence of the model to a global optimum rather than a local one 96,97 . We obtained the importance of the evaluated predictors, which was equal to the sum of the weights / probabilities of the models in which the variable appeared. Then, a variable with large weights that appeared in a set of candidate models received a high importance value 98 . We set the parameter level to 1 to include only the main effects, meaning that 28 represented the 256 possible models that were generated in the set of considered candidates. We used the function 'crit = aicc' (AICc or AIC corrected) to select the model and the multi-model inference, which were based on different predictors, accounted for all possible models that corresponded to different biological hypotheses, and varied in the degree to which they fit the available data 99 .
The second analysis included a stepwise phylogenetic regression using the calibrated phylogeny and the databases for the 257 species that included the body size of males, females, and SSD in conjunction with the previously defined climatic variables (Supplementary Table S2). The stepwise search mode was backward elimination starting from the full model, including all independent variables (BIO1 + BIO2 + BIO3 + BIO4 + BIO12 + BIO15 + DEM). The models were compared based on their AIC values.
Ethics statement. All procedures related to sampling, processing, and obtaining genetic material of some amphibian species included in this study were evaluated and approved by the Research Ethics and Safety Unit of the Pontificia Universidad Católica de Chile (protocol ID 180611004). Collection permits were approved under the scientific research framework agreement No: 200 of the Corporacion Autonoma Regional de la Frontera. The methods developed in this manuscript do not involve live animal experimentation.

Data availability
The datasets generated and/or analysed during the current study are available in the Supporting information section. Supplementary Table S3 have the Pristimantis species included in the phylogenetic analyses with the respective Genbank accession numbers linked to the NCBI website.


s2fieldsofstudy: Biology Medicine

text:
eRNAs Identify Immune Microenvironment Patterns and Provide a Novel Prognostic Tool in Acute Myeloid Leukemia

Background: Enhancer RNAs (eRNAs) play an essential role in tumorigenesis as non-coding RNAs transcribed from enhancer regions. However, the landscape of eRNAs in acute myeloid leukemia (AML) and the potential roles of eRNAs in the tumor microenvironment (TME) remain unclear. Method: Gene expression data collected from The Cancer Genome Atlas (TCGA) project were combined with Histone ChIP-seq so as to reveal the comprehensive landscape of eRNAs. Single-sample gene set enrichment analysis algorithm (ssGSEA) and ESTIMATE were employed to enumerate immune cell infiltration and tumor purity. Results: Most prognostic eRNAs were enriched in immune-related pathways. Two distinct immune microenvironment patterns, the immune-active subtype and the immune-resistant subtype, were identified in AML. We further developed an eRNA-derived score (E-score) that could quantify immune microenvironment patterns and predict the response to immune checkpoint inhibitor (ICI) treatment. Finally, we established a prognostic nomogram combining E-score and other clinical features, which showed great discriminative power in both the training set [Harrell’s concordance index (C index): 0.714 (0.651–0.777), p < 0.0001] and validation set [C index: 0.684 (0.614–0.755), p < 0.0001]. Calibration of the nomogram was also validated independently. Conclusion: In this study, we systematically understood the roles of eRNAs in regulating TME diversity and complexity. Moreover, our E-score model provided the first predictive model for ICI treatment in AML.

INTRODUCTION
Acute myeloid leukemia (AML) is the most common leukemia of the hematopoietic system (De Kouchkovsky and Abdul-Hay, 2016;Estey, 2018), which is characterized by rapid proliferation and invasion of immature myeloid cells in the bone marrow (BM) and peripheral blood (PB), resulting in hematopoietic failure of normal blood cells. Chemotherapy and hematopoietic stem cell transplantation are the first-line treatment for AML patients (Döhner et al., 2010;De Kouchkovsky and Abdul-Hay, 2016;Tyner et al., 2018a;Estey, 2018), though only one-third of adult AML patients could get durable remission (Vago and Gojo, 2020). A high risk of relapse is considerable even after remission (Estey, 2018;Thol and Ganser, 2020). However, the cellular and molecular characteristics of AML have been investigated to promote current therapies (Yang and Wang, 2018;Vago and Gojo, 2020). For example, high counts of regulatory T cells were observed in both PB and BM of AML patients and inhibited effector T cells more effectively than in control (Szczepanski et al., 2009;Williams et al., 2019). Significantly, the immunosuppressive microenvironment of BM was closely associated with a poor prognosis (Vago and Gojo, 2020), which provided new insights into therapeutic targets. The combination of ICI treatment and chemotherapy showed great advantages. A more prolonged median overall survival (OS) of 6.3 months was yielded in Nivolumab combined with azacitidine group, compared with the azacitidine group in 70 older, relapsed or refractory (R/R) AML patients (Daver et al., 2019). Another PD1-blocking antibody, pembrolizumab, with decitabine or azacitidine, showed a similar response and survival advantage in R/R AML patients (Lindblad et al., 2018;Gojo et al., 2019). Nonetheless, a randomized phase II study of durvalumab with azacytidine verse azacytidine alone did not show a significant advantage in 214 newly diagnosed AML patients (Zeidan et al., 2019), suggesting new models to predict the response of immunotherapy are underexplored. eRNAs are transcribed in specific enhancer regions and considered as non-coding RNAs. It has been widely proven that eRNAs play essential roles in tumor proliferation, invasion, and migration. Several eRNAs, such as NET1e in breast cancer  and CCAT1 in squamous cell carcinoma cells (Jiang et al., 2018), showed strong interactions with target genes and contributed to cancer progression. Recent studies identified eRNAs as a better marker for active enhancers than H3K27ac, the histone modification mediated by P300 (Tyssowski et al., 2018). Due to their highly tissue-specific and disease state-specific characteristics, eRNAs have been demonstrated as an excellent marker to identify phenotypes such as ICI treatment response compared with mRNAs (Chen and Liang, 2020), which has been validated in a broad range of diseases including lung cancer (Qin et al., 2020;Ma et al., 2021), glioblastoma (Guo et al., 2021), and hepatocellular carcinoma . However, systematic studies for eRNAs in AML are still under exploration. Zhang et al. (2019) established a systematic strategy for investigating eRNA in the RNA-seq of large patient cohorts.
The FANTOM5 database has generated approximately 65,000 high-quality and in vivo-transcribed enhancer loci with CAGE-seq in hundreds of primary cells and cell lines (Andersson et al., 2014). The ENCODE and Blueprint databases provided a wealth of histone modification locus information, including H3K4me1, H3K4me3, and H3K27ac annotation (Inoue et al., 2017). The TCGA, Beat AML, and GEO databases provided detailed data on gene expression profiles and clinical annotations of AML patients (Weinstein et al., 2013;Tyner et al., 2018b). In conclusion, comprehensive high throughput sequencing provided detectable methods to identify eRNAs in AML. We screened the TCGA database to identify eRNAs closely correlated with clinical prognosis and TME abundance, thereby exploring the characteristics of different immune subtypes associated with eRNAs in AML patients. A scoring system was established via machine learning to investigate the association between immune microenvironment patterns and existing clinicopathological stratifications, demonstrating that the TME can supplement the current risk classification. As we expected, the scoring system named E-score was found to correlate with the response to anti-PD1/PDL1 agents computationally. Finally, we developed a clinical prognostic model for AML patients based on the E-score, which showed accurate predictive power.

Data Collection and Preprocessing
We downloaded RNA-seq BAM files, mRNA expression matrixes, and clinical features of acute myeloid leukemia (TCGA-LAML) from the TCGA data portal (https://portal.gdc. cancer.gov/) (Weinstein et al., 2013). The validation datasets were downloaded from the Beat AML (http://www.vizome.org/) (Tyner et al., 2018a) and Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/, GSE37642 (Li et al., 2013), GSE12417 (Metzeler et al., 2008), and GSE10358 (Tomasson et al., 2008)). To ensure the stability of the prediction, patients with survival time <1 month were excluded. The PROseq BAM files, Histone ChIP-seq, and P300 ChIP-seq data of the Kasumi-1 cell line were downloaded from the GEO database (GSE83660 and GSE100446) (Zhao et al., 2016;Xu et al., 2019). From the Encyclopedia of DNA Elements (ENCODE) portal (https://www.encodeproject.org/), we downloaded the Hi-C interactions data of K562 and GM12878 cell lines with the following identifiers: ENCFF996XEO and ENCFF355OWW (Inoue et al., 2017). The Histone ChIP-seq data of the mononuclear cells of the bone marrow were collected from the Blueprint project (IHEC portal: https://epigenomesportal.ca/ ihec/index.html, IHECRE00000277) (Adams et al., 2012;Stunnenberg et al., 2016). The ribosome-deleted RNA-seq data of two AML patients' BM samples were obtained from GSE87285 (Loke et al., 2017). All detailed information on datasets mentioned above is listed in supplementary table 1. filtered out eRNA loci overlapping with known annotated elements and 1 kb extension from both the transcription start site (TSS) and transcription end site (TES) to eliminate introns, exons, promoters, rRNA repeats, and other transcribed elements. The 3 kb of the middle locus of the enhancer was designated as the potential eRNA region . RNA-seq data of the TCGA-LAML cohort were mapped to eRNA regions by featureCounts and normalized to reads per million (RPM) (Liao et al., 2014). eRNA with RPM ≥1 in at least one sample was defined as an active eRNA Zhang et al., 2021). The enhancer loci from all the human tissues and cell lines in the database were included in the analysis. Liftover was used to convert genomic regions from various genomic versions to the hg19 genome version (Rhead et al., 2010). The annotation of enhancer loci was supported by R package annotatr (Cavalcante and Sartor, 2017). The histograms of promoters, enhancers, and random loci were created by HOMER (http://homer.ucsd.edu/homer/) (Heinz et al., 2010). Immunerelated eRNAs (ir-eRNAs) were defined as at least correlated with an immune cell type (absolute Spearman's correlation Rs > 0.3, p < 0.05). We used univariate Cox regressions (p < 0.01, hazard ratio (HR) > 1.1 or <0.9) to identify prognostic eRNA.

Identification of Enhancer-Gene Pair and Validation of Chromatin State Analysis
The pair of eRNAs and protein-coding genes was identified based on close distance (1 MB) and the corresponding expression level (absolute Spearman's correlation Rs > 0.3 and FDR <0.05) concurrently . We selected the defined eRNA loci, promoter of protein-coding genes, and random loci of the genome (0.5-4.5 kb away from the region of the defined enhancer and excluded any overlapping annotation) to exhibit the genome coverage. The Histone and P300 ChIP signal was visualized by Deeptools, and the HiC data were illustrated by HiCExplorer (Wolff et al., 2020).

Unsupervised Clustering Analysis
To investigate distinct patterns, we used the R package ConsensusClusterPlus for unsupervised clustering analysis (Wilkerson and Hayes, 2010). The number of iterations was set at 1,000 to ensure clustering stability. The other parameters, clusterAlg, distance, and linkage, were set as pam, Pearson, and ward.D.

Gene Set Variation Analysis, Gene Set Enrichment Analysis, and Gene Ontology Annotation
GSVA is a nonparametric, unsupervised method for assessing alternations in the activity of signaling pathways and biological processes using expression datasets. We compared signaling pathways across different ir-eRNA patterns using the R package GSVA (Hänzelmann et al., 2013). The gene set "c2. cp.kegg.v6.2. symbols" from the Molecular Signatures Database (MSigDB) was used as the predetermined dataset. Statistical significance was defined as an adjusted p value of less than 0.05. We used the R package clusterProfiler to annotate the ir-eRNA-related gene for Gene Ontology (GO) analysis by biological processes (BPs) (p < 0.05, FDR <0.05) and the differentially expressed genes (DEGs) for gene set enrichment analysis (GSEA) by the Kyoto Encyclopedia of Genes and Genomes (KEGG) with a default parameter.
Evaluating Cell Infiltration in the TME By generating a gene set using cell-specific markers, single-sample gene set enrichment analysis (ssGSEA) allows us to assess the number of immune cells in AML BM samples. A previous study yielded gene sets for labeling each kind of TME-infiltrated immune cells, such as activated dendritic cells, activated CD8 + T cells, macrophages, active CD8 T cells, and regulatory T cells (Barbie et al., 2009;Charoentong et al., 2017). The relative abundance of TME cells in each sample was represented by the estimated abundance score in silico. In addition, the ESTIMATE method was also used to compute immune and stromal scores, as well as tumor purity, corroborating the ssGSEA approach.

Identifying DEGs in Different ir-eRNA Patterns
Based on prognostic ir-eRNAs, AML patients were split into two different expression patterns. To determine pattern-related genes, we used the R package limma to compare the expression profiles between the two patterns. An empirical Bayesian algorithm using moderated t-tests was used to carry out this strategy. The adjusted p value of 0.05, determined with Benjamini-Hochberg correction, was set as the significant criteria for detecting DEGs (Ritchie et al., 2015).

Generating Gene Signature and E-Score
We developed a scoring system to assess the immune pattern of individual AML patients. First, we standardized the expression of DEGs recognized from different immune subtypes in all TCGA samples. Patients were divided into two subgroups by DEGs through unsupervised clustering analysis. Next, we used a univariate Cox regression model for prognostic analysis of each gene in the signature (p < 0.01, HR > 1.1 or <0.9). We calculated the correlation between different subtypes and signature genes (Pearson correlation analysis). The positively correlated genes were defined as gene set A, and the negatively correlated genes were defined as gene set B. We then performed principal component analysis (PCA) to construct pattern-related signatures. The benefit of this strategy was that it concentrated the scores on the dominant network of gene signatures while downweighting contributions from genes that did not correlate with other set members. Principal component 1 of each gene set was extracted as PC1. A and PC1.B to calculate E-score, which was defined in an approach similar to that of the Genomic Grade Index (GGI) ( Van 'T Veer et al., 2002;Zeng et al., 2019): where i is the expression of a gene related to each signature gene set.

Statistical Analysis
Spearman's rank correlation analysis evaluated the correlation coefficient between the expression of ir-eRNAs and TMEinfiltrated immune cells (Sokol et al., 2017). We utilized the R package survminer to identify the cutoff point between each E-score subgroup based on the connection between the E-score and the overall survival of the patients. To get the maximum rank statistic, all possible cutoff points were evaluated repeatedly. Patients were separated into high-and low-Escore groups using log-rank statistics to minimize batch effects for maximum selection (Kerschke et al., 2020). The Kaplan-Meier curves were used to create survival curves for prognostic analysis, and the significance of differences was determined using the log-rank test. The hazard ratios (HRs) and p values of ir-eRNAs and signature genes were calculated using a univariate Cox regression model. To find independent prognostic factors, we employed a multivariate Cox regression model. The R package forestplot was applied to visualize the results of univariate and multivariate Cox regression. Receiver operating characteristic (ROC) curves were used to assess the sensitivity and specificity of E-score, and the area under the curve (AUC) was determined using the R package survivalROC. The R package RCircos was used to illustrate the chromosomal location of the ir-eRNAs and target genes. The significance criterion for all statistical comparisons was 0.05, and all comparisons were two sided. R 4.0.2 was used to conduct all statistical analyses.

Availability of Data and Source Code
All expression profiles and clinical annotation files were downloaded from public databases. Our code is available on request.

The Landscape of eRNAs in Acute Myeloid Leukemia
The overview of the study is depicted in Figure 1. We started by gathering approximately 65,000 enhancer loci from the FANTOM5 project, a robust enhancer database including enhancers from most human tissues and cell lines (Andersson et al., 2014). Then, we identified 13,251 enhancer RNAs in the TCGA-LAML cohort, none of which overlapped with annotated genomic regions. Active enhancer loci were typically described as eRNA transcription and unique epigenomic modification patterns, such as high H3K4me1, high H3K27ac, low H3K4me3, and high P300 levels (Klemm et al., 2019). To verify the reliability of enhancer loci in AML, we then investigated the epigenomic modification patterns of 13,521 eRNA regions, which showed enhancer-like epigenomic patterns compared with random loci and protein-coding gene loci (Supplementary Figure S1). In addition, the distribution of RNA transcription abundance revealed obvious transcription peaks in a classic manner Frontiers in Molecular Biosciences | www.frontiersin.org May 2022 | Volume 9 | Article 877117 among enhancer regions, demonstrating that the expression of enhancer regions was not the transcription background noise (enhancer loci vs. random loci, Figure 2A). On average, the expression level of eRNAs was less than that of protein-coding genes (enhancer loci vs. random loci and transcription start site (TSS) of protein-coding gene vs. random loci, Figures  2A,B). Enhancers usually play roles as cis-elements (Klemm et al., 2019). Previous studies revealed that eRNAs were excellent markers for active enhancers, so eRNAs usually correlated with adjacent genes. To illustrate the characteristics of eRNAs, pairwise Pearson correlation analysis revealed that eRNAs were more correlated with neighboring genes in a distance-related way, which indicated that eRNA has a stronger association with the closest gene than distal ones ( Figure 2C). The closest coding genes near active eRNAs (RPM >1) exhibited a higher expression level than those near nonactive eRNAs (p < 0.001, Wilcoxon test; Figure 2D). Similar to previous studies, we observed that most eRNAs were mapped to interact with one adjacent gene (51.20%, 1,171/3,459, Figure 2E), and most genes were only matched to a single eRNA as well (59.74%, 2,664/4,459, Figure 2F). In conclusion, we demonstrated the reliability of the defined enhancer loci and presented the landscape of eRNA in AML. Our results and recent studies suggested that eRNA levels were intimately associated with the expression of nearby genes (Andersson et al., 2014;Chen and Liang, 2020).

Prognostic ir-eRNAs Correlated With the Tumor Microenvironment in the Bone Marrow
We then concentrated on eRNAs correlated with prognosis, which may play crucial roles in tumor progression and proliferation. We filtered out 258/13,251 prognostic eRNAs (Supplementary Table  S2) by univariate Cox regression. To explore the potential functions of prognostic eRNAs, co-expression analyses were performed. Correlated mRNAs in the 1 Mb region were identified as potential target genes (Rs > 0.3, Spearman correlation analysis). Functional enrichment analysis revealed that the bitter taste-related pathway and several vital oncogenic signaling pathways such as PI3K-Akt-mTOR and MAPK were enriched in GO enrichment analysis ( Figure 3A) (Nepstad et al., 2020;Kim et al., 1999). Surprisingly, prognostic eRNAs were most concentrated on immune-related pathways, such as myeloid leukocyte activation and granulocyte migration ( Figure 3A). Therefore, we suggested that AML blasts shape immune microenvironment patterns through eRNAs and related genes. We evaluated the immune cell abundance of each sample by ssGSEA and filtered 72 eRNAs correlated with immune cell abundance (Supplementary Figure S2), which were defined as ir-eRNAs. Due to the diverse expression of ir-eRNAs in different samples, we explored whether ir-eRNAs exhibited distinct expression patterns in AML. Unsupervised clustering analysis was performed, and two clusters were identified (Supplementary Figures S3A-E, Figure 3B). One cluster defined as an immuneresistant subtype showed a higher infiltration level of suppressive Frontiers in Molecular Biosciences | www.frontiersin.org May 2022 | Volume 9 | Article 877117 6 cells, including dendritic cells, macrophages, regulatory T cells, and myeloid-derived suppressive cells (MDSCs). The other cluster defined as the immune-active subtype showed a lower level of suppressive cells and a trend of a higher level of activated T cells without statistically significant at the 5% level ( Figure 3C). The stromal and immune scores determined by ESTIMATE were also higher in the immune-resistant subtype (stromal score, p < 0.001, Wilcoxon test; immune score, p < 0.001, Wilcoxon test; Figure 3D), while the immune-active subtype had a higher tumor purity (p < 0.001, Wilcoxon test; Figure 3E), in line with the work of Yan et al. (2019). Cell clustering revealed complicated correlation networks of 23 immune cells and their prognostic significance by the log-rank test ( Figure 3F). In addition, Figure 3G shows the detailed hazard ratio of each immune cell by univariate Cox regression analysis. Plasmacytoid dendritic cell and macrophage, previously reported as cancer-promoting cells, had hazard ratio >1 (HR, 2.06 (1.314-3.228), p = 0.001; HR, 1.694 (1.083-2.648), p = 0.023, respectively) (Miari et al., 2021;Xiao et al., 2021). The Kaplan-Meier curves validated the poor prognosis of several immune cells, and the representative cell types are shown in FIGURE 4 | S100-enhancer and S100-eRNA play a role in shaping the immunosuppressive microenvironment. (A) GSEA of immune-related pathways confirmed different pathway activation models between immune-resistant and immune-active subtypes. (B) The chromosomal locations of potential target genes of ir-eRNAs identified in our study. Red ones were S100-eRNA and five S100 proteins. (C) Higher expression level of the nine genes regulated by S100-enhancer and S100-eRNA in the eRNA immune-resistant subtype compared with the immune-active subtype (ns, p > 0.05, *p < 0.05; **p < 0.01; ***p < 0.001, Wilcoxon test). (D) Correlations between the S100 protein family and immune checkpoint-related genes (Spearman correlation test, *p < 0.05; **p < 0.01). (E) The survival relevance of the S100 protein family in pan-cancers (red, risky factors; blue, protective factors; and grey, p > 0.05). (F) Schematic illustration of the putative mechanism by which S100-enhancer regulates immune infiltration. A phrase-separated structure is shown.
Frontiers in Molecular Biosciences | www.frontiersin.org May 2022 | Volume 9 | Article 877117 7 Supplementary Figures S3F-P (log-rank test, cutoff determined by R package survminer). Furthermore, patients of the immuneresistant subtype had a significantly shorter overall survival time than patients of the immune-active subtype (log-rank test, p < 0.001, Figure 3H). In general, we discovered that two clinically significant immune subtypes characterized by distinct ir-eRNAs profiles existed in AML, suggesting that ir-eRNAs played vital roles in modulating the immune microenvironment of BM. The immune-resistant subtype was characterized by a high level of suppressive and cancer-promoting immune cells. The immune-active subtype was accompanied by a relatively high level of protective immune cells. Although the absolute abundances of activated T cells, including CD4 T cells and CD8 T cells, were not significant between the two types, fewer suppressive immune cells in the immune-active subtype provided an immune microenvironment with less immune evasion and T-cell dysfunction.

S100-Enhancer and S100-eRNA May Induce a Poor Immune Microenvironment in AML
Based on the abovementioned findings, we applied the GSVA algorithm to explore pathways that differed between two immune subtypes. The result showed that the immune-resistant subtype was highly associated with innate immune response pathways, while the immune-active subtype was related to the T-cell receptor signaling pathways (Supplementary Figure S4A). The GSEA of immune-related pathways further confirmed that the NF-κB pathway and other innate immune response pathways were highly enriched in immune-resistant subtypes ( Figure 4A). Surprisingly, the PDL1 expression and PD1 checkpoint pathway showed a close connection to the immune-resistant subtype, leading to the following analysis to predict the response to ICI therapy. To further explore the possible biological roles of each immune subtype, we identified 1,057 ir-eRNA-related differentially expressed genes (DEGs) between immuneresistant and immune-active subtypes by the limma package (Supplementary Table S3). To investigate whether the DEGs directly correlated with ir-eRNAs, we subsequently screened genes with correlation index >0.3 in the approximately 1 Mb adjacent regions of ir-eRNAs Zhang et al., 2021) and identified 148 potential target genes (Supplementary Table S4). The Venn diagram illustrated that 55 genes showed both different expression levels between two clusters and correlations with ir-eRNAs, which suggested that these genes might be directly activated by immune-related enhancers (Supplementary Figure S4B). Most of these genes were reported as immune-related genes in AML Brenner and Bruserud, 2018;Borrego, 2013). Among these genes, we focused on five genes (S100A4, S100A6, S100A8, S100A9, and S100A12) from the S100 protein family that was widely reported in AML. The five genes were correlated to a specific eRNA transcribed from a chromosomal region (chr1: 153488459-153488748) (supplementary table 4, Figure 4B). We defined this enhancer and eRNA as S100-enhancer and S100-eRNA, respectively. In addition, Supplementary Figure S4C indicates the compact correlation among the five genes. Figure 4C reveals that the five genes exhibited consistent higher expression in the immune-resistant subtype (S100A4 p < 0.001, S100A6 p = 0.001, S100A8 p = 0.014, S100A9 p = 0.006, S100A12 p = 0.056, Wilcoxon test). Brenner and Bruserud (2018) summarized that the S100 protein family promoted tumor growth by enhancing a suppressive immune microenvironment in AML. Thus, we examined the relationship between S100-eRNA and immune cells and found a strong correlation with suppressive immune cells, including MDSC, macrophage, and regulatory T cell, as was expected (Supplementary Figure S4D). Figure 4D demonstrates that the S100 protein family tightly correlated with the immune checkpoint-related gene expression, suggesting the potential efficacy in predicting response to ICI treatment. Because S100 proteins were positively correlated with the immune-resistant subtype, a relatively poor prognosis would be observed. Supplementary Figures S4E-I demonstrate the hypothesis. We then expanded the S100 protein family to pan-cancer survival analysis. S100 proteins represented general risky factors in the pan-cancer scale, except for S100A6/A8/A9 in Sarcoma (SARC) and S100A4/A8 in Head and Neck Squamous Cell Carcinoma (HNSC) (univariate regression analysis, Figure 4E). These results suggested that the S100enhancer and S100-eRNA generally facilitated a suppressive inflammatory microenvironment, resulting in a worse outcome via the S100 protein family and related pathway, as Figure 4F illustrates. The illustration of the genomic region revealed that the S100-enhancer and correlated genes were in the same transcription activation domain (TAD) in leukemia cells, which suggested S100-enhancer interacted with nearby genes and induced their transcription ( Figure 5). Additionally, PROseq demonstrated an active transcription event in the enhancer region, and histone modifications confirmed the epigenomic characteristics of S100-enhancer.

Immune Infiltration Patterns Indicate a Novel Clinical Classification and Potential Response to Immunotherapy
To explore the roles of immune microenvironment patterns in predicting clinical prognosis and potential response to immunotherapy, we created a gene signature based on prognostic DEGs between immune-resistant and immuneactive subtypes (246/1,057 prognostic DEGs, Supplementary  Table S5). Random Forest regression, a machine learningbased algorithm, was used to get representative genes (66/246 after Random Forest regression, Supplementary Table S5). We then performed an unsupervised clustering analysis with the 66 signature genes and obtained two gene clusters ( Supplementary  Figures 5A-D). The immune infiltration and survival analysis based on both ir-eRNAs and signature genes were similar, thus supporting the robustness of immune subtypes ( Supplementary  Figures 5E-G). The result of PCA revealed the well-defined characteristics of immune subtypes as well (66 signature genes, Figure 6A).
The abovementioned results demonstrated that ir-eRNA alteration had a significant impact in generating distinct TME patterns. These findings, however, were limited to the population Frontiers in Molecular Biosciences | www.frontiersin.org May 2022 | Volume 9 | Article 877117 8 level and could not properly predict the immune profiles of AML at the individual level. Given the alternative heterogeneity and diversity of ir-eRNA-mediated immune profiles, we developed a scoring system based on these phenotype-related genes to quantify the TME pattern of individual patients in AML through the PCA algorithm, which was named eRNA-derived score (E-score) (see Methods, Figure 6B). E-score well distinguished the two immune subtypes as was expected (p < 0.001, Wilcoxon test; Figures 6B,C). The high-E-score group had a better prognosis than the low-E-score group (log-rank test, p < 0.001, Figure 6D). Nearly all immune-active subtypes were stratified into the high-E-score group ( Figure 6C). FIGURE 5 | S100-enhancer and correlated genes in the same transcription activation domain (TAD). Upper panel: Hi-C analysis of the genomic region containing the S100-enhancer gene in the K562 cell line, compared with the control cell line. Arrow pointed out the enhancer location. Medium panel: UCSC gene annotations and PRO-seq analysis within ±0.5 Mb region of the S100-enhancer gene, including five correlated genes located in the same TAD (bold). Lower panel: PRO-seq analysis and enhancer-related epigenomic modification analysis (H3K4me1, H3K4me3, and H3K27ac) within the S100-enhancer gene region. S100-enhancer is highlighted in yellow.
Frontiers in Molecular Biosciences | www.frontiersin.org May 2022 | Volume 9 | Article 877117 9 Frontiers in Molecular Biosciences | www.frontiersin.org May 2022 | Volume 9 | Article 877117 AML was commonly classified by the French-American-British (FAB) classification and the morphology, immunology, cytogenetics, and molecular biology (MICM) classification via genetic, immunophenotypic, pathological, and clinical features. CALGB cytogenetics-based risk classification and ELN 2017 risk classification provided widely used predictive models based on the integration of cytogenetics and mutational status (Döhner et al., 2017;Tallman et al., 2019). We examined the efficacy of E-score in distinguishing various clinical characteristics by the TCGA cohort and another independent cohort, Beat AML. It was found that the high-and low-E-score groups showed distinct clinical characteristics ( Figure 6E). When favorable CALGB risk patients were almost concentrated in the high-E-score group, the low-E-score group contained a greater proportion of intermediate and poor CALGB risk patients ( Figure 6F). The M3 subtype was categorized exclusively into the high-E-score group, while M4 and M5 subtypes were most concentrated in the low-E-score group ( Figure 6G). Furthermore, PML-RARα translocations only existed in the high-E-score group ( Figure 6H). On the contrary, CBFB-MYH11, GATA2-MECOM, and MLLT3-KMT2A only occurred in the low-E-score group ( Figure 6H). Differences in other clinical features, including the response to chemotherapy, ELN2017 risk classification, and WHO classification, are shown in Supplementary Figures 5H-J. The group with both low E-score and intermediate/poor risk showed poor prognosis compared with the other groups, which indicated that a good proportion of patients with intermediate/poor risk showed relatively better outcomes ( Figure 6I). We have mentioned that the immune-resistant subtype was enriched in the PD−L1 expression and PD−1 checkpoint pathway ( Figure 4A). The immune-resistant TME with abundant suppressive immune microenvironment generally adopted a better response to immune checkpoint therapy, such as the anti-PD1/PDL1 agents, in solid tumors (Turan et al., 2018). We predicted the response to immunotherapy in the two groups. The Microsatellite Instability (MSI) score and T-cell dysfunction score of Tumor Immune Dysfunction and Exclusion (TIDE) were higher in the low-E-score group, suggesting a predictive better response to immunotherapy ( Figure 6J). We also detected other nonredundant immune checkpoints and checkpoint-related gene expression between the two clusters ( Figure 6K). Several checkpoint molecules, including IFNG, TNFRSF9, CTLA4, CD86, TIM-3, and PD-L2, had significantly increased expression in the low-E-score group, and their effects on T-cell dysfunction and immune escape have been verified in AML (Hobo et al., 2018;Kikushige, 2021). However, there was no significant difference, but only a trend, in the expression of PD1 and PDL1, implying that the TME modulation is complicated and understudied. On the other side, our result partly explained the reason why current clinical trials of anti-PD1 agents exhibited controversial results. A group of intermediate/poor-risk patients with low E-score might be the potential population benefiting from the ICI treatment, which requires clinical trials for further validation. Based on the aforementioned results, our E-score system firstly proposed a quantizable system to predict response to ICI therapy in AML.

E-Score Associated With Genetic Risk Stratification Serves as an Independent Prognostic Factor in AML
To discover whether E-score could provide satisfactory predictive efficacy in AML, we proved that E-score was an independent prognostic factor in 4 validation cohorts (Beat AML, GSE37642, GSE12417, and GSE10358). Consistent survival differences were observed in all the validation cohorts mentioned above (Beat AML cohort, log-rank test, p = 0.021; GSE37642, log-rank test, p = 0.020; GSE12417, log-rank test, p = 0.003; and GSE10358, logrank test, p = 0.019; Figures 7A-D). Both univariate cox regression and multivariate cox regression confirmed the independent predictive power of age, CAGLB risk classification, and E-score (TCGA-LAML cohort, Supplementary Figure 6A). Log-rank regression of these clinical features also revealed that age and CAGLB risk classification were independent predictors (TCGA-LAML cohort, Supplementary Figures 6B-F). E-score showed a high area under the curve (AUC) in a time sequence (AUC 1year = 0.741, AUC 3year = 0.730, AUC 5year = 0.810, Figure 7E). Additionally, compared with other clinical characteristics (gender, age, race, FAB subtype, and CAGLB risk classification), the E-score exhibited the best capability at predicting OS in the TCGA cohort ( Figure 7F). Given that the E-score showed an excellent predict power in OS, we constructed a nomogram integrating E-score with other clinical features (TCGA-LAML cohort as the training cohort, Figure 7G). The calibration curves of the nomogram for 1-year, 3-year, and 5-year survival probability exhibited excellent consistency with the ideal performance, indicating the high accuracy of our nomogram (TCGA-LAML cohort, Figure 7H). The nomogram discrimination as evaluated by C index was 0.714 (95% CI = 0.651-0.777, p < 0.0001), which was superior to that of the E-score alone (0.674, 95% CI = 0.608-0.740, p < 0.0001), CAGLB risk classification (0.594, 95% CI = 0.526-0.663, p < 0.0001), and age (0.661, 95% CI = 0.596-0.726, p < 0.0001) in the TCGA-LAML cohort. The nomogram was externally confirmed in the validation cohort, the Beat AML cohort. The C index was 0.684 (95% CI = 0.614-0.755, p < 0.0001; Beat AML cohort), indicating a significant discriminative ability. The calibration curve presented a good agreement for 1-year OS rates in the validation cohort (Supplementary Figures 6G). However, FIGURE 6 | Patient stratification based on the CALBG risk classification (favorable, intermediate, and poor) and E-score group (high and low) and corresponding prognosis (log-rank test p < 0.001). (J) Enrichment of immunotherapy-related markers in high-and low-E-score groups (ns, p > 0.05, *p < 0.05; **p < 0.01; ***p < 0.001, Wilcoxon test). (K) Expression of immune checkpoints and related genes in two E-score groups (ns, p > 0.05, *p < 0.05; **p < 0.01; ***p < 0.001, Wilcoxon test). eRNA.cluster.A and gene. cluster.A, immune-resistant subtype; eRNA.cluster.B and gene. cluster.B, immune-active subtype.
Frontiers in Molecular Biosciences | www.frontiersin.org May 2022 | Volume 9 | Article 877117 11 calibration curves exhibited underestimated 3-and 5-year OS rates ( Supplementary Figures 6H-I). In all, the E-score system had outstanding predict power in both OS and response to ICI therapy, which improved the current predictive models and made it potential for clinical practice.

CONCLUSION
This study identified immune-related prognostic eRNAs (ir-eRNAs) for AML. Clustering of ir-eRNAs revealed two distinct immune subtypes, the immune-resistant subtype and the immune-active subtype. S100-eRNA was found to lead to worse outcomes with the activation of the S100 protein family and subsequently induced a suppressive TME. A novel eRNA-derived scoring system (E-score) was developed and showed great stratification power for patients in terms of the immune microenvironment and clinical characteristics. Finally, we established a prognostic nomogram with satisfactory accuracy by integrating E-score and other clinical features. This study provided a comprehensive understanding of the impact of eRNAs on shaping the immune microenvironment in AML and developed eRNA-derived tools to predict prognosis.

DISCUSSION
In this study, we identified that AML bone marrow samples can be classified into immune-resistant and immune-active subtypes by ir-eRNAs. Then, we revealed the underlying rationality of the immune classification by an enhancer-gene pair. Finally, we developed a prognostic model based on immune scoring, E-score, for predicting overall survival and response to immunotherapy in AML patients.
Our study revealed that the high suppressive immune microenvironment with high counts of regulatory T cells and MDSCs correlated with poor prognosis. In this paper, we used eRNA, a non-coding RNA, for effective clustering and suggested that enhancer plays an important role in the construction of the BM microenvironment. We cautiously inferred the conclusion based on the fact that eRNAs are excellent markers for active enhancers that facilitated the transcription of target genes (Andersson et al., 2014). We demonstrated that a large amount of eRNAs was correlated with the abundance of immune cells in AML. Moreover, the putative genes in the vicinity of these ir-eRNAs have been experimentally validated to correlate with the invasion and progression of AML (Borrego, 2013;Brenner and Bruserud, 2018;Li et al., 2020). Bulk-level RNA-seq data of the current public database were the confounding profiles of the different cell types, which compromised the power of mRNAs to predict the molecular and clinicopathological features of tumor cells. The essential TME-promoting mRNAs expressed in AML blasts could be averaged by the background expression of other cell types in the TME. However, eRNA has been reported to associate with tumor traits, such as response to immunotherapy, more powerfully because of the strong cell lineage specificity (Chen and Liang, 2020). E-score derived from eRNA classification had a higher AUC (0.810 vs. 0.68) compared to other risk scores derived from immune cell abundance as a predictor of OS , implying that eRNA-based immune subtypes had better efficacy to predict the OS of AML patients.
The interaction between tumor cells and the microenvironment together promoted tumorigenesis and tumor progression. The S100 protein family provided a suppressive immune microenvironment in a variety of tumors, promoting tumorigenesis and metastasis (Brenner and Bruserud, 2018). Among the genes with a high correlation with S100-eRNA, our attention was focused on S100A8/A9, a heterodimer that is overexpressed in many cancers and associated with suppressive TME. S100A8/A9 is an activator of monocytes and macrophages and leads to neutrophil infiltration, correlating with a poor prognosis of patients. S100A8/A9 promotes tumor development and invasiveness by increasing leukemic cell growth through RAGE and secretion of pro-inflammatory cytokines through the TLR4-NFκB pathway (Hiratsuka et al., 2008). Kaikkonen et al. (2013) have reported that eRNA takes part in the activation of the TLR4-NFkB pathway in macrophages by promoting enhancer-promoter looping and regulating mediator recruitment, which increases our confidence in that the S100-eRNA-S100 protein family axis is potentially involved in shaping the immune microenvironment.
In this study, we first systematically concluded the prediction of response to immunotherapy linked to immunophenotyping in AML. Several clinical trials about checkpoint inhibitors therapy have been pursued in AML and got controversial conclusions (Ghosh et al., 2020). Identification of specific populations who might benefit from ICI therapy would be the ideal strategy to solve this problem. Our study firstly developed a method to predict the potential population that benefited from ICI therapy and implied that other ICIs, including TIM-3, might be more attractive targets. We collected general predictors of response to ICI therapy. T-cell exclusion score was excluded due to suspiciously overestimated MDSCs. In addition, the tumor mutation burden (TMB) was also omitted because TMB is generally quite low in AML. Further validation of our model requires further detailed survival data of the immunotherapy cohort. Our results suggest E-score may be a promising predictor for the choice of immunotherapy, especially for M4/M5 patients with high E-score and abundant suppressive immune infiltration.
It is important to note that limitations existed in our study. First, because of the similarities between leukemia cells and myeloid-derived cells such as MDSC and monocytes, a possible overestimation of MDSC and monocyte existed. We carefully based the conclusion on the abundance of these cells. Whether eRNA functions as an initiation factor of vital gene transcription needs further functional experiment because eRNA may be only a marker of enhancer active transcription instead of a functional factor. The enhancer-gene pairs were determined by co-expression analysis as in other studies , which was the best strategy for large clinical cohorts currently. More advantages of single-cell RNA-seq or single-cell ATAC-seq could improve the accuracy of prediction and support further analysis.
In conclusion, our study distinguished AML into immuneresistant and immune-active subtypes by ir-eRNA for the first time. The important roles of these eRNAs in the clinic require investigation, which may provide new insight into the pathogenesis of AML. The predictive model with high accuracy provides new tools for physicians. This study increases the understanding of the important roles of enhancers on the BM microenvironment and might provide novel therapeutic targets.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
MC and YZ were the principal investigators who designed and conceived the study and obtained financial support. ZJ and JL analyzed the data and prepared the dataset. ZJ and KD wrote the manuscript. All authors have read, revised, and approved the final manuscript.

FUNDING
This work was supported by the National Natural Science Foundation of China (NSFC) (81970106).


text:
Melatonin ameliorates lung cell inflammation and apoptosis caused by Klebsiella pneumoniae via AMP-activated protein kinase

Klebsiella pneumoniae is a Gram-negative bacterium and the causative agent of several life-threatening nosocomial infections, including pneumonia. K. pneumoniae induces acute lung injury and inflammation in humans that require immediate hospitalization and treatment. Therefore, attenuation of K. pneumoniae-induced inflammation is necessary for the survival of patients. This study investigated the mechanisms by which melatonin abrogated K. pneumoniae-induced inflammation and apoptosis of lung cell lines, HLF-1 and BEAS-2B. Our results showed that in vitro infection of HLF-1 and BEAS-2B cells by K. pneumoniae significantly induced inflammation and apoptosis increased elevated levels of IL-6, CXCL1, CXCL2, and caspase-9 mRNA. However, these effects were abrogated by melatonin treatment. Infection with K. pneumoniae significantly increased the expression of AMP-induced protein kinase (AMPK). Furthermore, AMPK silencing significantly abrogated the suppression of inflammation and apoptosis in melatonin-infected K. pneumoniae lung cells. Melatonin could alleviate K. pneumoniae infection-induced inflammation in three-dimensional lung spheroids. In conclusion, our study demonstrated that melatonin abrogated K. pneumoniae-induced inflammation and apoptosis in lung cells through AMPK. Our study demonstrated the potential of melatonin for therapy against K. pneumoniae infections including pneumonia.

Introduction
Klebsiella pneumoniae (K. pneumoniae) is a Gram-negative bacterium that causes life-threatening nosocomial infections, including pneumonia (Wyres et al. 2020). In humans, K. pneumoniae is mainly located in the gastrointestinal tract and some locate in the nasopharynx, through which bacteria can enter the bloodstream or other tissues and cause infection (Wyres et al. 2020). To date, four virulent factors have been identified, including pili, capsule, lipopolysaccharide (LPS), and iron carriers (Martin et al. 2018). Importantly, K. pneumoniae is assembled by adhesins, type 1 and type 3 pili, and it promotes bacterial adhesion to epithelial, immune cells, and abiotic surfaces (Martin et al. 2018). K. pneumoniae induces acute lung injury in humans and is associated with high mortality rates of nearly 50% (Wyres et al. 2020). It is difficult to treat pneumonia because K. pneumoniae is highly resistant to carbapenems, the last resort antibiotic against Gram-negative bacteria (Chong et al. 2018). K. pneumoniae acquires antimicrobial resistance (AMR) genes that cause multidrug resistance through de novo mutations or the acquisition of plasmids and transferable genetic elements (Vardakas et al. 2018). Therefore, treatment of K. pneumoniae infections require the development of drugs that bypass AMR genes.
Pathogenic bacteria induce acute or chronic inflammation of various human organs, including the intestine and lungs (Cao 2017). For example, invasive pathogens induce inflammation by disrupting the intestinal niche and gut equilibrium generated by the commensal bacterial communities (Abt and Pamer, 2014). The in vitro intestinal organoid model demonstrated that Salmonella infections-induced intestinal cell death via inflammation (Yin and Zhou 2018). Chronic inflammation by Helicobacter pylori causes peptic ulcer disease and gastric cancer (Lehours and Ferrero 2019). Inflammation could also induce serious cell injury through different mechanisms, including the mTOR/STAT3 pathway activation, autophagy, and induction of oxidative stress (Ryter 2021). Therefore, attenuation of inflammation is critical to alleviating disease caused by bacterial infections.
Melatonin (N-acetyl-5-methoxy-tryptamine) is a sleep hormone that was originally discovered in plants (Bhattacharya et al. 2019). In humans, melatonin influences several biological activities such as circadian rhythms, mood, sleep, body temperature, locomotor activity, food intake, retina physiology, sexual behavior, seasonal reproduction, and the immune system (Bhattacharya et al. 2019). Increasing evidence has shown that melatonin exerts important regulatory effects on pathogens, including viruses and bacteria (Bahrampour Juybari et al. 2020;Anderson and Reiter 2020). For example, melatonin also exhibits broad antiviral activity against COVID-19 (Bahrampour Juybari et al. 2020), the influenza virus (Anderson and Reiter 2020), and the Ebola virus (Anderson et al. 2015). Several studies have also shown the antibacterial properties of melatonin. For example, melatonin overcomes mobilized colistin resistance (MCR) of Gram-negative pathogens . Melatonin also demonstrates anti-inflammatory and cytoprotective functions. For example, in the mouse model, melatonin alleviated adipose inflammation by increasing alpha-ketoglutarate levels and diverting adipose-derived exosomes to macrophages (Liu et al. 2018). Melatonin has also been reported to regulate T-cell differentiation, interfere with T/B cell interaction, and attenuate the production of pro-inflammatory factors, achieve its antioxidant action via specific receptors to exert immunomodulatory effects (Yildirim et al. 2019). In bacterial meningitis (BM), melatonin was also found to exert several protective effects in BM through various mechanisms, including the immune response, antibacterial activity, protection of the integrity of the blood-brain barrier (BBB), free radical scavenging, anti-inflammatory activity, activation of signaling pathways, and effects on the gut microbiome (Gao et al. 2021). However, the effects of melatonin on inflammation and cell death induced by bacterial infections have not been reported. Therefore, it is attractive to investigate the effects of melatonin on bacteria-induced inflammation.
In the present study, we investigated the effects of melatonin on lung epithelial cell inflammation and cytotoxicity induced by infection with K. pneumoniae and the underlying mechanisms. Importantly, an advanced in vitro threedimensional (3D) lung spheroid model was used in the study, since it has a substantial advantage compared to conventional two-dimensional (2D) cell line models, including improved mimicking of a 3D structure in vitro. We believe that our study provides new information for understanding the pathogenesis of K. pneumoniae infections.

Preparation of bacteria
K. pneumoniae was obtained from the Biological Resource center of Tianjin First Central Hospital and grown as described previously (Anand et al. 2020). K. pneumoniae was cultured in the microbial culture room and grown in Luria Bertani (LB) broth with constant shaking at 37 °C or solid media containing LB and 1.5% agar in an incubator maintained at 5% CO 2 , pH 7.0-7.4, and 37 °C. The bacteria grown overnight in the exponential phase were harvested by centrifugation at 2700× g for 20 min and resuspended in PBS. Culturing of three-dimensional (3D) lung spheroids 3D lung spheroids were purchased from Innovation Biotechnology (Tianjin, China). Culture of 3D lung spheroids was used previously published method (Kumari et al. 2020). Briefly, a mixture of 3000 HLF-1 cells and 3000 BEAS-2B cells were plated on 96-well ultra-low attachment plates (Corning) in complete media (DMEM media containing 10% FBS, 1% penicillin-streptomycin, and 1% 2 mM L-glutamine) supplemented with 4% Matrigel™ (Corning, New York, USA). Media was refreshed every 4 days by carefully aspirating 50% of the well volume and replacing with fresh complete media containing 1% Matrigel™. The 3D spheroids were re-fed with fresh media every 2 or 3 days.

Bacterial infection
HLF-1 or BEAS-2B cells or 3D lung spheroids grown in DMEM medium at 5% CO 2 and 37 °C were serum-starved for 18 h and infected with the bacterial inoculum prepared in DMEM. Briefly, 80-90% confluent HLF-1 or BEAS-2B cells were seeded in 24-well plates at a density of 2 × 10 5 cells per well or 3D lung spheroids were seeded in a 4% Matrigel coated 24-well plate at a density of 50 spheroids per well and infected with K. pneumoniae at a multiplicity of infection (MOI) ranging from 5:1 to 10:1. The cells or 3D lung spheroids were then centrifuged for 4 min at 200× g at 22 °C and incubated between 2 and 24 h at 37 °C and 5% CO 2 in a humidified incubator. HLF-1 or BEAS-2B cells infected with K. pneumoniae or 3D lung spheroids were treated with 1 μM melatonin for 24 h.

CCK-8 assay of cell viability
Lung cell viability was analyzed using the CCK8 assay kit (Beyotime, Cat. No. C0037, Shanghai, China) according to the manufacturer's instructions. Briefly, 1 × 10 4 HLF-1 or BEAS-2B cells were seeded in 96-well culture plates, infected with K. pneumoniae, and treated with or without 1 μM melatonin for specified periods of time. The cells were then incubated for 4 h with 100 μL CCK-8 and the absorbance was measured at 450 nm. Cell viability (%) = OD (treated cells)/OD (control cells) × 100.

Reverse transcription-polymerase chain reaction (RT-qPCR)
Total RNA was extracted from lung cells using the RNeasy plus Mini Kit (Qiagen, Maryland, USA). Reverse transcription (RT) was performed with 500 ng of total RNA using the BeyoRT™ III cDNA synthesis kit (Cat. No. D7178M, Beyotime, Shanghai, China). qPCR was performed with the cDNA samples using the TB Green qPCR Kit (Takara) in the 7500 Fast Real-Time PCR System (Applied Biosystem, CA, USA). The qPCR reaction conditions included initial denaturation at 94 °C for 5 min followed by 40 cycles of  Table 1.

Flow cytometry
For analysis of apoptosis rates, after treatments, 3D lung spheroids were harvested to 15 centrifuge tubes, followed by spinning down at 100 g for 5 min. Then, supernatant was discarded, followed by adding 5 mL TrypLE™ Express Enzyme (1×, phenol red, Thermo Fisher Scientific GmbH, Waltham, Massachusetts, US) and incubating in water bath (°C) for 10 min. Then, dissociation was terminated by adding 5 mL FBS. Then, single cells were washed with PBS for three times, followed by being resuspended in 1× annexin binding buffer (10× , 0.1M hepes (pH 7.4), 1.4M NaCl, 25 mM CaCl 2 ) and incubated for 5 min at room temperature. Afterwards, cells were pelleted and stained with annexin-V-FITC (Miltenyi Biotech, Bergisch Gladbach, Germany) and propidium iodide (PI) solution (Miltenyi Biotech, Bergisch Gladbach, Germany) according to the manufacturer's manual. They were incubated for 15 min at room temperature before an additional 300 μl Annexin V Binding Buffer was added to each sample. Gating parameters were determined using Veh/PBS-treated cells as a negative control. Unstained cells were run to account for auto fluorescence. Stained cells were analyzed by flow cytometry using a FACSLyric system (BD Biosciences, Franklin Lakes, New Jersey, US). For instrument setup, the excitation wavelength was set to 488 nm, and emission wavelength was set to 530 for annexin-V, and the excitation wavelength was set to 488 nm, and emission wavelength was set to 585 for PI.

Statistical analysis
GraphPad Prism (version 9, GraphPad Software, San Diego, CA, US) was used for statistical analysis. Data are shown as mean ± SEM. One-way analysis of variance (ANOVA) was used to compare the data between multiple groups. P < 0.05 was considered statistically significant.

K. pneumoniae-induced inflammation in lung cell lines
We analyzed the inflammation induced by K. pneumoniae in the lung cell lines. Infection of the HLF-1 and BEAS-2B lung cell lines with K. pneumoniae (MOI = 10) significantly increased the mRNA levels of inflammation-related genes such as IL-6 ( Fig. 1A, D), CXCL1 (Fig. 1B, E), and CXCL2 ( Fig. 1C, F). This demonstrated that K. pneumoniae-induced inflammation in lung cells.

Infection with K. pneumoniae significantly reduced the viability of lung cells
Lung injury and inflammation are the common characteristics of lung diseases (Daniel et al. 2019). Therefore, we investigated the effects of K. pneumoniae infection on lung cell viability. K. pneumoniae infection (MOI = 5:1 and 10:1) significantly decreased the viability of HLF-1 and BEAS-2B cells ( Fig. 2A, C). Furthermore, HLF-1 and BEAS-2B cells infected with K. pneumoniae (MOI = 5:1 and 10:1, respectively) showed significantly higher levels of caspase-9 mRNA, a pro-apoptotic gene (Fig. 2B, D). This demonstrated that infection with K. pneumoniae-induced apoptosis in lung cells.

Melatonin-alleviated inflammation induced by K. pneumoniae in lung cells
Melatonin reduced inflammation in several cell types (Zhang et al. 2020). IL-6, CXCL1, and CXCL2 mRNA expression in HLF-1 and BEAS-2B cells infected with K. pneumoniae treated with 1 μM melatonin was significantly lower compared to the corresponding controls ( Fig. 3A-F). The dose and treatment time were referenced in previous studies (Lee et al. 2018;Rahman et al. 2005). This demonstrated that melatonin treatment significantly reduced the inflammatory response induced by K. pneumoniae infection in lung cell lines.

Melatonin treatment increased the viability of K. pneumoniae-infected lung cells
Next, we investigated the effects of melatonin on lung cell apoptosis induced by K. pneumoniae infection. Treatment of K. pneumoniae-infected lung cell lines with 1 μM melatonin significantly reduced caspase 9 expression levels ( Fig. 4A, C). Furthermore, melatonin significantly increased the viability of HLF-1 and BEAS-2B cells infected with K. pneumoniae (Fig. 4B, D). These results demonstrated that melatonin significantly increased the viability of the K. pneumoniae-infected lung cell lines.

AMPK expression was elevated in lung cells infected with K. pneumoniae
AMPK is a key regulator of cellular metabolism that plays a significant role in modulating bacterial infections ). Therefore, we investigated whether infection by K. pneumoniae modulated AMPK expression. K. pneumoniae (MOI = 10:1) significantly increased AMPK mRNA expression in HLF-1 and BEAS-2B cells at different time points (2, 4, and 8 h) (Fig. 5A, C). Furthermore, K. pneumoniae infection significantly increased AMPK protein levels in HLF-1 and BEAS-2B cells (Fig. 5B, D). This showed that K. pneumoniae infection significantly increased AMPK mRNA and protein levels in lung cells.

Melatonin-suppressed K. pneumoniae-induced inflammation in lung cells via AMPK
We then investigated the role of AMPK in the anti-inflammatory effects of melatonin in K. pneumoniae-infected lung cells. We silenced AMPK in HLF-1 and BEAS-2B cells by transfection of cells with siRNA against AMPK (AMPK KD) (Fig. 6A,  D). The qRT-PCR results indicated that the IL-6 and CXCL1 mRNA levels mRNA levels were significantly higher in AMPK-silenced lung cell lines treated with melatonin infected with K. pneumoniae compared to lung cell lines treated with melatonin infected with K. pneumoniae (Fig. 6B, C, E, F). This demonstrated that melatonin suppressed inflammation caused by K. pneumoniae infection in the lung cells via AMPK.

Melatonin-suppressed apoptosis in lung cells infected with K. pneumoniae by AMPK
Next, we analyzed the effects of melatonin on the viability of K. pneumoniae-infected lung cells and the corresponding mechanism. Treatment of HLF-1 and BEAS-2B cells infected with K. pneumoniae with 1 μM melatonin significantly reduced caspase9 mRNA expression levels (Fig. 7A, C) and increased cell viability (Fig. 7B, D) compared to HLF-1 and BEAS-2B cells infected with K. pneumoniae without melatonin treatment.
Interestingly, effects of melatonin on the promotion of caspase9 expression and reduction of cell viability by K. pneumoniae infection were found in both cell types (Fig. 7). This showed that melatonin alleviated apoptosis of K. pneumoniae-infected lung cells via AMPK.

Melatonin-alleviated inflammation caused by K. pneumoniae in 3D lung spheroids
3D in vitro models have been shown to more accurately mimic in vivo physiology (Yin et al. 2021;Clevers 2016). Thus, we evaluated inflammation due to K. pneumoniae infections and the effects of melatonin on K. pneumoniae infections-induced inflammation in 3D lung spheroids (Fig. 8A). K. pneumoniae infections increased IL-6 ( Fig. 8B), CXCL1 (Fig. 8C), and CXCL2 (Fig. 8D) mRNA expression in 3D lung spheroids. Importantly, 1 μM melatonin significantly reduced IL-6 ( Fig. 8E), CXCL1 (Fig. 8F), and CXCL2 (Fig. 8G) mRNA expression levels in 3D lung spheroids. To further confirm the effects of K. pneumoniae infections on lung cell apoptosis and the effects of melatonin on lung cell apoptosis, the flow cytometry assays were performed to confirm cell apoptosis induced by K. pneumoniae infections, while melatonin (1 μM) attenuated K. pneumoniae infection-induced cell apoptosis in 3D lung spheroids (Fig. 8H). This demonstrated that melatonin treatment significantly reduced the inflammatory response induced by K. pneumoniae infection in a 3D lung spheroid model.

Discussion
K. pneumoniae was first described by Carl Friedlander in 1882 as an encapsulated bacterium isolated from the lungs of patients who died from pneumonia. K. pneumoniae is the third leading cause of healthcare-associated infections in the United States (9.9%) behind Clostridium difficile and Staphylococcus aureus (Martin and Bachman 2018). In the present study, we demonstrated that infection with K. pneumoniaeinduced inflammation and cytotoxicity in two lung cell lines, HLF-1 and BEAS-2B, by increasing the expression levels of AMPK. This study also demonstrated that melatonin alleviated inflammation and cell death in lung cell lines in response to infection by K. pneumoniae. Furthermore, this study demonstrated that AMPK mediated the effects of melatonin in alleviating inflammation and cell death in K. pneumoniae-infected lung cell lines. We also confirmed that melatonin could alleviate inflammation induced by K. pneumonia infections in 3D lung spheroids.
The immune system plays a key role in protecting the human body from disease-causing pathogens such as bacteria and viruses. In addition, infections from foreign pathogens such as bacteria activate the immune response or inflammation. Pathogenic oral bacteria and related toxins induce inflammation at distant sites in the body (Martin and Bachman 2018;Han and Wang 2013). Wong et al (2014) demonstrated prostatic inflammation and fibrosis in a mouse model of chronic bacterial infection. In patients with Crohn's disease, abnormal ileal bacterial colonization is associated Fig. 8 Melatonin-alleviated inflammation caused by K. pneumoniae in 3D lung spheroids. A Morphology of 3D lung spheroids. Infection with K. pneumoniae increased the expression of IL-6 B, CXCL1 C and CXCL2 D mRNA in 3D lung spheroids. Melatonin (1 μM) significantly reduced IL-6 E, CXCL1 F and CXCL2 G mRNA expression in 3D lung spheroids; flow cytometry indicating that K. pneumoniae infections-induced cell apoptosis, while melatonin (1 μM) could attenuate K. pneumoniae infection-induced cell apoptosis in 3D lung spheroids H ◂ with altered innate immune responses to invasive bacteria (Barnich et al. 2013). Wu et al (2018) demonstrate that the oral pathogen Salmonella enterica induced intestinal inflammation to eliminate competitors to gain exclusive access to fructose-asparagine, an essential nutrient for their growth. In the lungs, innate immunity plays an important role in the elimination of several pathogens without adverse immunopathological effects and in the absence of adaptive immunity (Dulek et al. 2014). In the lungs, bacterial pathogens, including K. pneumoniae, can induce immune responses that result in increased expression of innate immune genes such as interleukin-17A (IL-17A), IL-17F, and IL-22 (Dulek et al. 2014). In the present study, K. pneumoniae-infected lung cells showed significantly higher mRNA and protein expression of pro-inflammatory genes such as IL6, CXCL1, and CXCL2 (Fig. 1). Hoogerwerf et al. (2012) demonstrated that K. pneumoniae infection was closely associated with the inflammatory response in the lung epithelium. Therefore, it is essential to determine the mechanisms underlying the inflammatory responses during K. pneumoniae infection of lung cells. Furthermore, K. pneumoniae infection is also associated with severe lung injury (Dulek et al. 2014). In the present study, we confirm that K. pneumoniae infectioninduced cellular apoptosis in lung cell lines (Fig. 2).
In humans, melatonin is produced by the pineal glands. However, melatonin is widely produced naturally in several living organisms, from bacteria to humans . Melatonin exerts a variety of biological functions, including circadian rhythms, sleep, immune responses, and others . Several studies have shown that melatonin protects against bacterial infections. Liu et al (2020) demonstrated that melatonin substantially increased the activity of colistin against mobilized colistin-resistant (MCR)-positive Gram-negative pathogens. Melatonin is a powerful antioxidant molecule that can reduce oxidative stress and up-regulate antioxidative defense mechanisms that can influence innate immunity against bacterial pathogens . Xu et al (2019) demonstrated that melatonin exerted its protective effects against polymicrobial sepsis in the mouse model by promoting the development of the neutrophil extracellular trap (NET), which contributed to antibacterial activity during polymicrobial infection. In the present study, we showed that melatonin suppressed inflammation and apoptosis in lung cells infected with K. pneumoniae (Figs. 3, 4, and 8).
AMPK is a highly conserved serine/threonine kinase that functions as a key metabolic sensor in eukaryotic cells and plays an important role in the adaptive response to exercise (Spaulding and Yan 2022). AMPK is associated with several biological processes, including tumorigenesis, virus infections, metabolism, and others (Prantner et al. 2017). AMPK also regulates host defenses against bacterial pathogens and is closely associated with innate immunity (Prantner et al. 2017). In this study, we demonstrated that melatonin alleviated inflammation and apoptosis related to K. pneumoniae infection in lung cells via AMPK (Figs. 6,7,8).
In conclusion, our study demonstrated that K. pneumoniae infection-induced inflammation and apoptosis in lung cell lines, which was associated with increased AMPK expression. Furthermore, melatonin alleviated K. pneumoniae infection-related inflammation and apoptosis of lung cell lines via AMPK. These results suggested that melatonin is a promising therapeutic agent for alleviating lung pathology and pneumonia related to K. pneumoniae infection.