The relationship between temperature (T) and species diversity is one of the most fundamental issues in ecology (Allen et al., 2002; Tittensor et al., 2010; Peters et al., 2016). Studies on the biotas of forests, grasslands, wetlands, continental shelves, the open ocean, and even the estuarine and coastal areas have all shown the significant correlations between T and diversity (Ray, 1991; Yasuhara and Danovaro, 2016; Henseler et al., 2019). Under the current circumstances of T change in most oceanic regions, their correlations have received more extensive attention from ecologists and climatologists (Currie et al., 2004; Mora et al., 2013; Yasuhara et al., 2020). However, how exactly species diversity is related to T remains a central yet controversial question in ecological research (Danovaro et al., 2004; O’Hara and Tittensor, 2010; Yasuhara and Danovaro, 2016), although many studies and discussions have surrounded the T–diversity relationship for many years (von Humboldt, 1871; Pianka, 1966; Brown, 2014).

Previous studies have identified three apparent T–diversity relationships. (1) Positive: As expected by ecological theories (Allen et al., 2002; Evans et al., 2005), most empirical T–diversity relationships were positive (Rutherford et al., 1999; Tittensor et al., 2010). (2) Negative: Negative relationship had been found on both temporal and spatial scales (Danovaro et al., 2004; O’Hara and Tittensor, 2010), which was supported by climate change science (Cheung et al., 2009; Mora et al., 2013). (3) Unimodal: Yasuhara synthesized two relationships above and stated that T and diversity were more likely to be a unimodal relationship (Yasuhara and Danovaro, 2016; Yasuhara et al., 2020).

The above studies of the T–diversity relationship might be biased, because (1) the temperature range of the datasets was not sufficiently broad (e.g., Danovaro et al., 2004; Hunt et al., 2005); (2) only one diversity metric (generally species richness, S) has been considered (Tittensor et al., 2010; Yasuhara and Danovaro, 2016); and (3) few studies had explored the T–diversity relationship from the perspective of species abundance distribution (SAD), which can better represent species diversity than many other metrics (Mouillot et al., 2000; McGill et al., 2007; Connolly et al., 2009). Thus, it is difficult to draw an overall conclusion on how species diversity relates to T.

The main purpose of this study is to seek a more comprehensive understanding of how species diversity changes with T. A useful model system for this purpose is the microfossils of planktonic foraminifera, because of their complete abundance records, large temperature coverage, and good taxonomic resolution (Rutherford et al., 1999; Siccha and Kucera, 2017; Yasuhara et al., 2017). Thus, the Brown University Foraminiferal Data Base (BFD) (Prell et al., 1999) is used here to explore the correlations between T and species diversity, which records the absolute abundances of 37 extant morpho-species of planktonic foraminifera at 1,265 core-top sites of the global oceans. Species diversity of foraminifera is estimated by four most commonly used diversity metrics and a new SAD parameter (Su, 2016).

To explore the T–diversity relationship from multiple perspectives, the four most common metrics are used to measure diversity, namely, Species richness (S), Shannon’s index (H’), Simpson’s index (D), and Pielou’s evenness (J) (Pielou, 1975). Additionally, a new SAD parameter (Su, 2016) is also used as a diversity surrogate.

SST ranges from −0.28 ± 0.01 to 29.60 ± 0.05°C. SST at the breakpoint of S, H’, D, and J are 21.08, 14.18, 9.24 and 7.61°C, respectively ( Table 1 ). In the first segment of the piecewise regression ( Figure 1 ), S, H’, and J are positively correlated with T (r_{S 1} = 0.87, r_{H’ 1} = 0.79, r_{J 1} = 0.75). D is negatively correlated with T (r_{D 1} = −0.78). Four metrics are significantly correlated with SST (P< 0.05). In the second segment of the piecewise regression ( Figure 1 ), S and H’ decrease significantly as SST increases (r_S2 = -0.31, r_H'2= -0.14, P< 0.05). D and J are not significantly correlated with T (P > 0.05) ( Table 1 ). As T increases, D increases and J decreases (r_{D 2} = 0.05, r_{J 2} = −0.06).

The R ² of the fractal model to the BFD data is from 0.632 to 1.000, and more than 86% is above 0.8 (see Supplementary Table 2 ). SST at the breakpoint of p is 8.10°C. When SST< 8.10°C, p decreases pronouncedly as SST increases (r_{p 1} = −0.69, P< 0.05). When SST > 8.10°C, p decreases slightly with SST (r_{p 2} = −0.11, P< 0.05) and tends to be flat ( Figure 2 ). The SST–p relationships are significantly negative over the entire SST range (see the non-line regression in Figure 2 , P< 0.05).

Although the bootstrap simulation is used, the results are found to be mostly consistent (see Supplementary Table 1 ). Thus, only the relationships between T and each metric without the bootstrap simulation are presented ( Table 1 and Figures 1 , 2 ).

The correlation between temperature (T) and species diversity has fascinated and occupied scientists for centuries (von Humboldt, 1871; Pianka, 1966; Brown, 2014). Ecologists generally suggested that increasing T promoted higher rates of speciation or provided more tolerable living conditions leading to greater diversity (Allen et al., 2002; Currie et al., 2004; Tittensor et al., 2010). However, climatologists often implied that climate warming would have negative biological consequences triggering a decline in diversity (Cheung et al., 2009; Mora et al., 2013). Thus, how diversity relates to T has remained controversial (Yasuhara and Danovaro, 2016). Since coastal T may change rapidly in the coming decades (Zhu et al., 2021), a better understanding of the T–diversity relationship is one of the central problems in current research (Cheung et al., 2009; Brown, 2014; Yasuhara et al., 2020).

Yasuhara stated that “when considered over a sufficiently broad T range, the T–diversity relationship appears to be unimodal” (Yasuhara and Danovaro, 2016). This idea stems from a unimodal model, noting that fewer species can physiologically tolerate conditions in extremely cold or warm places than in other sites. Thus, S will show a peak at intermediate T (Currie et al., 2004). The well-known unimodal depth–diversity relationship was considered to support the unimodal model (Yasuhara and Danovaro, 2016; Yasuhara et al., 2020).

In this study, S, H’, and J show similar unimodal relationships with T ( Figures 1A, B, D ). Their relationships are positive below T at the breakpoint (r_{S 1} = 0.87, r_{H’ 1} = 0.79, r_{J 1} = 0.75) and negative over T at the breakpoint (r_{S 2} = −0.31, r_{H’ 2} = −0.14, r_{J 2} = −0.06). The T–D relationship is opposite to that of the other three metrics in two segments of the piecewise regression (r_{D 1} = −0.78, r_{D 2} = 0.05). However, D represents the species dominance in a community. The lower the D, the higher the diversity (Pielou, 1975). Thus, four metrics are not completely positive or negative with T. Species diversity measured by four metrics increases first and then decreases with increasing T. This is basically consistent with the expectation of the unimodal model (Yasuhara and Danovaro, 2016).

As noted before, the unimodal model is mostly derived based on the relationship between T and S (Yasuhara and Danovaro, 2016). However, S only represents one aspect of species diversity. It has a number of other components, such as dominance (D) and evenness (J) (Hill, 1973; Pielou, 1975; Willig and Presley, 2013). A better understanding of the T–diversity relationship needs to take these components into account (Pielou, 1975; Stevens and Willig, 2002; Willig and Presley, 2013). According to the piecewise regression results ( Table 1 ), T at the breakpoints of D and J are apparently lower than that of S. This means that dominance has started to increase and evenness has started to decrease, when S is still increasing with T. Although the present results based on a single database may not be very universal, they at least imply that S cannot describe some changes that have occurred in species diversity. Therefore, this study suggests that evenness and dominance should be considered when investigating the T–diversity relationship, as they cannot be reflected by S. Additionally, it is worth noting that the change of D and J over T at the breakpoint is less significant than that of S (P > 0.05). This means that the unimodality of dominance and evenness is not as obvious as that of S. Similar results have also been found in the study of New World bat (Stevens and Willig, 2002). Stevens and Willig characterized three aspects of diversity for bat communities and suggested that the evenness and dominance did not change significantly at the warmer lower latitudes (Stevens and Willig, 2002). Thus, whether dominance and evenness can be treated as the typical unimodal relationship also need to be further investigated.

Finally, many ecologists have suggested that the traditional diversity indices (e.g., H’, D, and J) only represent a synthetic characterization of relative abundance in a community (Pielou, 1975; Stevens and Willig, 2002; Willig and Presley, 2013), which may lead to an excessive simplification of diversity (Mouillot et al., 2000; McGill et al., 2007). Since SAD can reflect the relative abundance for each different species encountered within a community (Tokeshi, 1993; Mouillot et al., 2000; McGill et al., 2007; Ulrich et al., 2010), it is a more comprehensive description of diversity (Mouillot et al., 2000; McGill et al., 2007; Connolly et al., 2009). According to Eq. 4, p is the fractal parameter that determines SAD. The T–p relationship shows that the direction of variation of p is the same in two segments of the piecewise regression (r_{p 1} = −0.69, r_{p 2} = −0.11, P< 0.05) ( Figure 2 ) and their relationship is significant (P< 0.05, see the non-line regression in Figure 2 ). This is totally different from the other metrics (S, H’, D, and J). Accordingly, it is argued here that the fractal p has an asymptotic relationship with T rather than the typical unimodal relationship.

This study assumes that p with increasing T will gradually approach the minimum 1. Firstly, when p = 1, SAD (N_r /N ₁) is 1:1/2:1/3… (Eq. 4), which is consistent with Zipf’s law (Zipf, 1949). Secondly, in previous studies, Su suggested that Zipf’s law (1:1/2:1/3…) might be a general pattern of SAD (Su, 2018), which is supported by nearly 20,000 community samples. Thirdly, Su (2018) proposed two hypotheses to elucidate this pattern. H₁: Species diversity was determined by the entropy that increased with the energy transformation. H₂: The total assimilated energy of the community (E_T ) was finite. These two hypotheses provide a direction for the interpretation of the T–p relationship.

(1) N_T /N ₁ (Eq. 6) is an effective number of species diversity, which is related to Rényi’s entropy (Rényi, 1961; Hill, 1973). According to H₁, increasing entropy with T leads to the rise of N_T /N ₁ (Sethna, 2006). (2) According to H₂, N_T /N ₁ is finite as N_T is usually equivalent to E_T (Hutchinson, 1959; Brown, 2014). On the one hand, increasing N_T /N ₁ with T will cause p to decrease (Eq. 6). On the other hand, the finiteness of N_T /N ₁ determines that p ought to be higher than 1 (N_T /N ₁ converges when p > 1, Eq. 6), which means that the theoretical minimum of p is 1. Therefore, if H₁ and H₂ hold true, their combined effect will make p approach the theoretical minimum 1 with increasing T. Since a lower p means a higher diversity (Su, 2016); this study conjectures that species diversity with increasing T will rise and gradually approach the theoretical maximum, which is Zipf’s law (1:1/2:1/3…).

It should be emphasized that the results of this study depend on the accuracy and representation of the dataset. In fact, although BFD has been used in many studies of diversity (Rutherford et al., 1999; Tittensor et al., 2010), there are some important biases that need to be noted here. Firstly, the “pachyderma-dutertrei intergrade” category is not recognized in BFD (King and Howard, 2003). Secondly, calcium carbonate dissolution in deep-sea sediments of the Pacific and Indian Oceans may decrease the species richness of core-top assemblages (Rutherford et al., 1999). Finally, the representation of planktonic foraminifera to other organisms (especially terrestrial taxa) is still limited (Rutherford et al., 1999; Siccha and Kucera, 2017; Yasuhara et al., 2017). Thus, more details of the T–diversity relationship in planktonic foraminifera, and their pattern in other taxa, still need to be further verified.

When the correlation between T and diversity was mentioned in previous studies, most of the discussions revolved around the potential impact of T on S (Allen et al., 2002; Tittensor et al., 2010; Peters et al., 2016; Yasuhara and Danovaro, 2016). However, this study indicates that (1) the T–diversity relationship is not wholly captured by S and the other aspects of diversity (especially SAD) should be considered; (2) S and H’ support the typical unimodal relationship with T, while dominance and evenness do not have significant unimodality; (3) p with increasing T may decrease and gradually approach the theoretical minimum 1 (N_r /N ₁ is 1:1/2:1/3…). These points are the biggest difference between previous studies and this paper.

In the future, the relationship between T and the other aspects of diversity (e.g., evenness and dominance) should be discussed more over a broader T range. This can provide feasible guidance for dealing with the potential responses of diversity to global T changes. Additionally, the relationship between T and SAD should receive more attention and investigation. It is particularly noteworthy that the theoretical minimum of fractal p (or the theoretical maximum of diversity) can be supported by more empirical data. Finally, analyzing the T–diversity relationship from multiple aspects and exploring their underlying mechanisms are expected to enrich the theories of community ecology and biogeography.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material . The code supporting the results in the paper have been archived in Figshare (https://doi.org/10.6084/m9.figshare.14933301.v10).

This study does not involve related issues. Ethics approval was not required for this study.

QS conceived and designed this study. JG developed the mathematical models and performed statistical analyses. JG and QS wrote the manuscript. All authors contributed to the article and approved the submitted version.

This work was supported by the National Natural Science Foundation of China (Grant Nos. 42071137 and 41676113). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.