Diversity time-period and diversity-time-area relationships exemplified by the human microbiome

We extend the ecological laws of species-time relationship (STR) and species-time-area relationship (STAR) to general diversity time-period relationship (DTR) and diversity-time-area relationship (DTAR), and test the extensions with the human vaginal microbiome datasets by building 1460 DTR/DTAR models. Our extensions were inspired by the observation that Hill numbers, well regarded as the most appropriate measure of alpha-diversity and also particularly suitable for multiplicative beta-diversity partitioning, are actually in the units of effective species, and therefore, should be able to substitute for species in the STR and STAR. We found that the traditional power law (PL) model is only applicable for DTR at diversity order zero (i.e., species richness); at higher diversity orders (q = 1–4), the power law with exponent cutoff (PLEC) and power law with inverse exponent cutoff (PLIEC) are more appropriate. In particular, PLEC has an advantage over PLIEC in predicting maximal accumulation diversity (MAD) over time. In fact, with the DTR extensions, we can construct DTR and MAD profiles. To the best of our knowledge, this is the first comprehensive investigation of the DTR/DTAR in human microbiome. Methodologically, our DTR/DTAR profiles can characterize general diversity scaling beyond species richness, covering both alpha- and beta-diversity regimes across different diversity orders.

Even with the rapidly expanding revolutionary advances in microbial community ecology, to the best of our knowledge, testing of STR with the human microbiome data has not been systematically investigated except for perspective general discussions (Ma et al., Oliver et al., Shade et al.) 21,22,36 . In this article, we report our comprehensive test and extensions of the traditional STR with recently available dataset from a longitudinal study of the human vaginal microbiome (Gajer et al.) 37 .
Methodologically, we extend the traditional STR and STAR models to the general diversity-time relationship (DTR) and diversity-time-area relationship (DTAR) models, respectively, by using the Hill Numbers as a general diversity metric system, including both alpha-and beta-diversity. The choice of Hill numbers not only makes our study more comprehensive than existing studies of SAR/STR, but also makes it possible to extend the traditional STR to both alpha-and beta-diversities beyond species richness (i.e., the zero-order diversity in terms of the Hill numbers). We take advantage of the recent consensuses that Hill numbers offer the best measure of alpha-diversity and that multiplicative beta-diversity partitioning of the Hill numbers is more appropriate than the alternative additive partition (Ellison, Chao et al., Jost, Chiu et al., Chao & Jost) [38][39][40][41][42][43][44] . As a side note, as explained in following section of methods, the time (T) we discuss in this article should be termed time period or period, because what is usually accumulated is a period of time, rather than time points. That is, more appropriate terms should be diversity time period or diversity period. In consider of the traditional usage of species time in ecology, we use the term diversity time and diversity-time period interchangeably. Accordingly, the acronym DTR may be interpreted as either diversity time relationship or diversity time-period relationship in this article.
To the best of our knowledge, this is the first extension of the STR/STAR to general diversity-time scaling (DTR/DTAR) beyond species-richness level in terms of the Hill numbers. Previously, a few authors, notably Helmus & Ives 45 , Mazel et al. 46 have successfully extended the SAR to phylogenetic and functional diversities. Their extensions not only verified the applicability of SAR models beyond traditional species richness, but also found important applications in identifying more comprehensive conservation hotspots and predicting the impacts of habitat loss. However, we are not aware of any similar extensions to STR in the existing literature. We expect that our methodological extensions of STRs/STARs to general Hill numbers based DTRs/DTARs should not only enrich the theoretical modeling of the diversity scaling in terms of more comprehensive diversity profiles, but also offer important novel insights for their ecological applications. In addition, the Hill numbers, as the metric for a diversity profile, have advantages lacked by other existing diversity measures. For instances, the Hill numbers effectively capture the essential properties (such as rarity vs. commonness) of species abundance distribution (SAD) by computing the entropy at different orders (non-linearity levels) ( [38][39][40][41][42][43][44] . These advantages are naturally carried over to our general DTR/DTAR extensions. Theoretically, the proposed DTR/DTAR can characterize the temporal dynamics/spatial-temporal dynamics of community diversity. Furthermore, with our derived formula, the MAD (maximal accumulation diversity) over time can be predicted. Practically, Oliver et al. 21 emphasized that future importance of studying the microbial STR includes to distinguish between anthropogenic perturbations and underlying natural dynamics, and to offer ecological insights for clinical benefit in fighting bacterial infections. Their point makes particular sense in the case of the human vaginal microbiome (HVM) because the DTR can potentially be utilized to assess the bacterial temporal turnover or diversity dynamics of the HVM.

The Methods-Extensions from STR/STAR to DTR/DTAR
We use the following definitions and computational procedures to design, conduct, and interpret DTR and DTAR modeling analyses.
Definitions of the alpha-and beta-diversity measured in the Hill numbers. We adopt the Hill numbers to measure both alpha-and beta-diversities, and adopt multiplicative partition of the Hill numbers as beta diversities.
The Hill numbers, originally introduced as an evenness index in economics by Hill 47 who was apparently inspired by Renyi's 48 general entropy, did not receive the attention it deserves in ecology until recent years (e.g., Jost, Chao et al.) [39][40][41][42] : where S is the number of species, q is the order number of the Hill numbers. The Hill number is undefined for q = 1, but its limit as q approaches 1 can be derived as: The parameter q determines the sensitivity of the Hill number to the relative frequencies of species abundances. When q = 0, the species abundances do not count at all and 0 D = S, i.e., species richness. When q = 1, 1 D equal the exponential of Shannon entropy, and is interpreted as the number of typical or common species in the community. When q = 2, 2 D equal the reciprocal of Simpson index, i.e., The general interpretation of diversity of order q is that the community contains q D = x equally abundant species, or the so-termed effective number of species (or species equivalents such as OTUs). When the species are equally abundant, the relative abundance does not weigh in the computation of the Hill numbers, and consequently, the Hill numbers are determined by the number of effective species completely. Note that, although the assumption of equal relative abundance is rather rarely satisfied and the computation of Hill numbers is usually dependent on both the number of species and their abundances, the principle underlying the equivalence between the effective number of species and the Hill numbers still holds regardless of the assumption. It is this equivalence principle that inspired our extensions from STR to DTR (this article) and from SAR to DAR (Ma) 49 .
The diversity of a community can be measured with a series of Hill numbers, possibly plotted on a single graph as a 'continuous' function of the parameter q. Chao et al. 39 termed the series of plots "community diversity profile" that characterizes the species abundance distribution (SAD) of a community and offers complete information on its diversity. As explained later, we extend the concept of diversity profile to the DTR profile and MAD (maximal accumulation diversity) profile.
Recent advances (e.g., Jost, Ellison, Chao, Chao et al., Chiu et al., Chao & Jost) [38][39][40][41][42][43][44] have made a convincing case that Hill numbers and multiplicative beta-partition offer to date the most generally consistent and appropriate, yet simple solution to measure biodiversity. Jost 44 demonstrated that the partition of Hill numbers into independent alpha-diversity (within community) and beta-diversity (between communities) is necessarily multiplicative, i.e., This beta-diversity derived from the above partition takes the value of 1 if all communities are identical, and the value of N (the number of communities) when all the communities are completely different from each other (no shared species). In general, the beta-diversity defined in this multiplicative partition is, with Jost 44 words, "the effective number of completely distinct communities. " In this article, we compute diversities until q = 4, i.e., to the fourth order, which includes traditional species richness (q = 0), the exponential of Shannon index (q = 1), the reciprocal of Simpson index (q = 2) and two additional sets of diversity indexes for q = 3 and 4, respectively.
The DTR (Diversity Time Relationship) models. We were inspired by the observation that Hill numbers, well regarded as the most appropriate measure of the alpha-diversity, are actually in the units of species; indeed, they are referred to as the effective number of species or as species equivalents (such as OTUs), as explained in the previous section. Furthermore, Hill numbers capture the essential distribution properties of species abundances (rarity vs. commonness) by computing the entropy at different orders (non-linearity levels); they are particularly suitable for multiplicative beta-diversity partitioning, which is superior to its alternative-additive partitioning. Based on these inspirations, we postulate that the models of traditional STR should be extendable to model general DTR.
In consideration of the debates on the functional forms of SAR/STR, we selected three models, i.e., power law (PL) model as well as its two extensions, the PL with exponential cutoff (PLEC) and the PL with inverse exponential cutoff (PLIEC), from more than 20 SAR/STR models existing in the literature ( [11][12][13]19,50,51 to build the DTR models. We assume that the basic power law (PL) function widely adopted in the STR scaling study, known as power-law species scaling law can be extended to model the general DTR: where q D is diversity in the form of q-order Hill numbers, T is time, and c & w are parameters. Strictly speaking, T should be time period or period, because what is usually accumulated is a period of time, rather than time points. In consider of the traditional usage of species time in ecology, we use the term diversity time. Nevertheless, obviously more appropriate terms should be diversity time period or diversity period. We also extend two additional models, which were originally introduced to model SAR by Plotkin et al. 52 and Ulrich & Buszko 53 , to investigate the DTR. The first is the power law with exponential cutoff (PLEC), where d is a third parameter, and exp(dT) is the exponential decay term that eventually overwhelms the power law behavior at very large value of T. The justification for adding the exponential decay term is because diversity accumulation may saturate and even decline over time, and therefore, there should be a taper-off item to reflect the reality.
Another function similar to (6) but with Sigmoid shape, rather than convex as (6) is the following power-law function with an inverse exponential cutoff (PLIEC): q w PLEC has an advantage over the other two models in predicting the maximal accumulation diversity (MAD), derived below. The two models we here refer to as the PLEC and PLIEC were usually referred to as the first and second persistence model, or persistence function-I and persistence function-II, respectively (e.g., Tjørve).
We use the following log-linear transformed Eqns (8), (9), (10) to estimate the model parameters of Eqns (5), (6), (7), respectively: where all the symbols represent the exactly same variables (parameters) as those in Eqns (5), (6), (7). These linear models are computationally simple to fit, and we use both linear correlation coefficient (R) and p-value to judge the goodness of the model fitting, either of which should be sufficient to judge the suitability of the models to data. The taper-off parameter (d) in both the PLEC and PLIEC not only addresses a critique to the traditional power law for overestimating diversity (He & Hubbell) 16 , but also preserves the biological interpretations of the scaling parameter (slope z of SAR or w of STR) since d is primarily a revision to the other biologically less important parameter c (Tjørve) 13 . Furthermore, PLEC allows for the estimation of the maximal accumulation diversity (MAD) of the human microbiome over time times, as detailed below.
The DTAR (Diversity Time Area Relationship) models. Rosenzweig  Here we further extend it to general diversity time area relationship (DTAR). Both the extensions adopt the following bivariate power law (BPL) model: Its log-transformed linear form is: where q D is the diversity measured with the Hill numbers of q-order, A is the area, T is the time, z and w are scaling exponents corresponding to area and time, respectively, c is a parameter that is largely determined by sampling design but with limited biological meanings. Different from most studies in macro-ecology, where there is often a natural spatial sequence (or arrangement) among the communities sampled, the community samples in our HVM data are somewhat 'random' since there is not a natural order among the 32 individual subjects. To avoid the potential bias from an arbitrary order of the community samples, we totally permutate the orders of all 32 subjects under investigation, and then randomly select (without replacement) 50 orders of the communities generated from the permutation operation. That is, rather than building only one DTAR model for the dataset, we build 50 such models and take the average of their parameters as the final DTAR parameters. In the time dimension of each individual subject, we simply follow the natural (calendar) sequence of the time-series data.
Therefore, to build a DTAR model, one of the 50 randomly sampled sequences is used to determine the accumulation order of the area (individuals), and the diversity profile is then computed for each of the time-series point of each individual subject. If D ij (i = 1, …, 32, j = 1, …T i ; T i is the number of times that the i-th subject were sampled over time series) represents the diversity of the i-th individual at time j, its area (A) and time (T) for fitting the DTAR models [Eqns (11) and (12)] would be A = i, and T = j. In this study, to avoid potential bias caused by different number of samples for different individuals (T i ), we fix T to 25 times. Consequently, the total numbers of samples utilized to construct the DTAR models are fixed to 800 (32*25) times. Finally, the averages of the model parameters from the 50 times of DTAR fittings are adopted as the model parameters of the DTAR for the set of community samples under investigation.
Sampling and accumulation schemes for the DTR and DTAR models. Since we need to accumulate not only time (period) for the DTR and DTAR, but also area (individual) for the DTAR, we must devise the accumulation schemes for both area and time. We adopt the Scheiner 54,55 Type-III sampling scheme, noncontiguous quadrats grid. Strictly speaking, quadrats (the same size) are arranged in a regular but noncontiguous grid. Since the quadrat (individual subject), i.e., our body (habitats for our microbiome) is highly mobile, figuring out a definite spatial coordinate relationship, not to mention a grid of quadrats, is hardly possible. Nevertheless, Type-III is obviously the most natural choice we can make.
Carey 56 extended Scheiner's 54 typology to species time relationships (STR): these are nested, completely nested, and island approaches to construct STR models depending on their definitions of time, whether it is defined as an interval, a flow or a combination of both. In our study, the length of sampling time (period) T is defined as the total time (period) T from the start of the initial survey. The scheme is the nested construction, and time is defined as both an interval and a flow in terms of Carey's classification.
After selecting the sampling schemes, we need to specify the scheme to accumulate diversity based on the diversity formulae [Eqns (1-4)] listed above. Although the accumulation of species in the traditional STR is well defined and there is little ambiguity on how the species counts are computed once the scheme for time accumulation is decided, the computation of the diversity accumulation is still largely an uncharted area, and there may be more than one scheme to accumulate diversity across space or time, especially for the accumulation of beta-diversity.
To devise what we believe to be the most appropriate and also natural methods to accumulate diversity, we propose the following three principles. The first is to use the Hill numbers, or what Jost 44 termed the true diversity; the second is to follow the essence of SAR and STR, as captured by the word "accumulation" or "aggregate, " i.e., species (diversity) are accumulated for the accumulated times (periods); the third is that the diversity scaling SCiEntifiC RePoRTS | (2018) 8:7214 | DOI:10.1038/s41598-018-24881-3 model should be useful for predicting diversity at different accumulated times periods. We consider these three principles as axioms in traditional SAR or STR and we believe that any extension from STR to DTR should not violate them.
To accumulate alpha-diversity across times, assuming there are N community samples, at each accumulation step i = 1, 2, …, N-1, we simply add up all the OTU lists (rows) up to the accumulation step i in the standard OTU table of the communities, corresponding to the first i community samples, until the last sample (community N) is added when the accumulation is completed. For each of the aggregate (accumulated) community, we compute its alpha-diversity with Eqns (1) or (2) (for q = 1) with the added-up OTU lists for the accumulation step. The resulted pairs of accumulated alpha-diversities and times are regressed to fit the DTR model.
To accumulate beta-diversities across times, we start the computational procedure with two local communities (samples) by using the formulae specified by Eqns (1)-(4), from which the first beta-diversity value is computed for the two initial communities. For each newly added community (sample) at the accumulation step i (i = 3, 4, … N), we simply run the same computation procedure with a pair of communities consisting of the previously pooled i-1 communities as one regional community and the newly added i-th community as the other. With each newly added local community, we obtain a new Hill numbers of beta-diversity, until all N communities are accumulated for their beta-diversity. The series of beta-diversities are regressed with their respective times accumulated at each step. Obviously, this accumulation scheme calculates the accumulated beta-diversity of N communities (of individuals in the case of human microbiome), i.e., the difference between the pooled N-1 communities and the last community. It is also the maximum difference among N communities in terms of beta-diversity.
Statistical distribution of the DTR/DTAR parameters. We further investigate the statistical distributions of the DTR/DTAR scaling parameters [w, ln(c), d, T max , D max ] by fitting two contrastingly different statistical distributions: the normal distribution and power-law distribution. The former depicts a largely symmetric distribution of the scaling parameters, and the latter depicts an asymmetric (long-tail) probability distribution that has some unique properties not possessed by the normal distribution. For example, the power-law distribution usually suggests heavy heterogeneity or skewed data points. It has the so-termed "no-average" property, which means that the average of the power law distribution cannot represent majority of the data points because of the highly skewed long-tail. Such information should be particular valuable for further characterizing the individualities of the human microbiome, as well as for personalized medicine (Ma et al.) 57,58 .
Since the information on normal distribution can be readily found in standard statistics textbook (e.g., Gotelli & Ellison) 59 , we only list some basic information about the power law distribution below. Power law distribution has a probability density function as follows: where x is the random variable, x min is the minimum value of x, K > 1 and is the exponent or the scaling parameter of the power law distribution. K can be considered as a measure of asymmetricity (skewness) of the heterogeneity in the power law distribution. A comprehensive discussion on the power law distribution can be found in Clauset et al. 60 .
Predicting the MAD (maximal accumulation diversity) with DTR models. PLEC [Eqn. (6)] has an advantage over the other two models in predicting the maximal accumulation diversity. Ma 49 derived the maximum of PLEC, which can be adapted for DTR as follows: The necessary condition for Eqn. (6) to achieve maximum is its derivative equals zero, i.e., Solving Eqn. (14), when max q D may have a maximum in the following form: Theoretically, Eqns (15) and (16) can be utilized to predict the maximal accumulation diversity (MAD) of the human microbiome, whether it is alpha-or beta-diversity. Data Accessibility. The longitudinal cohort dataset of HVM used to test DTR/DATR is available at: http:// stm.sciencemag.org/content/4/132/132ra52.full.

Results
Before presenting the test results for our DTR/DTAR extensions, here we first briefly introduce the 16 s rRNA sequence datasets we used to perform the tests. The datasets we utilized are from a longitudinal study on the HVM (human vaginal microbiome) by Gajer et al. 37 with 32 healthy women at reproductive age offered to date the longest and also the most comprehensive longitudinal study of microbiome dynamics, conducted with 16s-rRNA sequencing technology. To reduce the influence of the noises in sequence data, we filtered out the samples that have a total number of lesser than 100 16s-rRNA reads because samples with too small number of reads are often due to the failure in amplicon amplification and may introduce significant bias in binning out the OTUs from the limited number of reads.

Diversity-Time Relationship (DTR) modeling.
We construct three DTR models (PL, PLEC, and PLIEC) for each of the 32 subjects, with their time-series data of their vaginal microbial communities, respectively. A total of 960 DTR models were built. To save page space, Tables 1 and 2 only list the partial results of alpha-diversity DTR or ADTR (for q = 0, 1) and beta-diversity DTR or BDTR (for q = 0), respectively. The online Supplementary Information (Suppl. Tables S1A and B) report the full results for diversity order q = 0-4, and a detailed examination of the modeling results is also presented there. Here, we summarize some important findings as explained below.
Performance and parameter range of the DTR models. Regarding the alpha-diversity DTR or ADTR (Table 1  and Table S1A), the traditional PL model did perform universally well (p-values < 0.001) when q = 0, i.e., the traditional STR models. Although the PLEC and PLIEC performed slightly better than the PL model when q = 0, we select the PL model in consideration of the parsimony principle. However, systematic failures occurred when PL model was fitted to the higher order diversities (q = 1, 2, 3, 4). Instead, the two extended PL models, PLEC and PLIEC, especially the former, perform very well (p < 0.01) with the higher order ADTRs. Therefore, we conclude that, although traditional PL model is not applicable to the DTR at the higher diversity orders, either PLEC or PLIEC can be utilized to model the DTR for higher diversity orders (q = 1-4), with PLEC performing best. Between PLEC and PLIEC, we preferred PLEC because it not only allows for fluctuated diversity with time, but also can be utilized to predict maximal accumulation diversity.
As to the beta-diversity DTR or BDTR (Table 2 and Suppl. Table S1B), the performances of the three models (PL, PLEC, and PLIEC) are similar with the previous described ADTR modeling. Similar to the ADTR, the traditional PL model indeed encountered difficulties at higher order diversities, but its performance is generally better than it did with the ADTR. The PL model successfully fitted to the zero-order BDTR of all 32 subjects (p < 0.0001), and to the first-order BDTR of 32 subjects with an average p-value = 0.024. It indeed had some failures at the higher order (q = 2 to 4). As in the case of alpha-diversity, both PLEC and PLIEC fitted to the BDTR at all diversity orders significantly better than the traditional PL model. We therefore conclude that both PLEC and PLIEC can be utilized to describe the BDTR at higher diversity orders, with a preference to PLEC given that it can predict maximal accumulation diversity. We adopted two modified PL models, PLEC and PLIEC, originally introduced to SAR modeling by Plotkin et al. and Ulrich & Buszko 52,53 , respectively (also see Tjørve) 13 . The former is termed persistence function-I and has maximum value, and the latter is termed persistence function-II. Essentially, PLEC and PLIEC can be considered as extensions to parameter c, rather than to w, i.e., exp( / ), respectively. Therefore, w is assumed to have the similar interpretation as its counterpart in the basic PL. PLEC and PLIEC, however, both behave very differently. The PLEC model asymptotically approaches cx w as x becomes small, whereas the PLIEC asymptotically approaches cx w as x becomes large. They were designed to remedy the potentially unlimited accumulation of species when the time approaches to infinity by introducing a taper-off exponent that may even produce asymptote. Indeed, both the models performed extremely well, usually outperforming the basic PL model, except in a handful case when the sample size was too small (3 or 4). As demonstrated previously, in the higher order diversity of the DTR, PLEC and PLIEC performed particularly well, while the basic PL failed.
The PLEC model has an advantage over the other two models in predicting the maximal diversity accumulation. In addition, we do concur with Plotkin et al. and Ulrich & Buszko 52,53 that a taper-off parameter may indeed be necessary. The issue may not be serious for DAR, because both the Earth and the diversity on Earth are finite anyway. But time is potentially infinite; hence the traditional PL would predict an ever increasing (if w > 0) or ever decreasing (if w < 0) diversity, which could hardly be realistic. The issue can be particularly troubling in the case of human vaginal microbial communities because of the periodic nature of menses, which may cause periodic changes of the vaginal microbiota. From this perspective, the failure of traditional PL in modeling the DTR of vaginal microbiota should be anticipated; accordingly, the success of PLEC and PLIEC with taper-off effect should also be expected.
The alpha DTR or ADTR has the average exponent w = (0.406, 0.073, 0.021, 0.010, 0.005) corresponding to the diversity order q = 0-4 (Table 1, Suppl. Table S1A); the beta DTR or BDTR has the average exponent w = (0.407, 0.018, 0.073, 0.074, 0.081), corresponding to the diversity order q = 0-4 (Table 2, Suppl. Table S1B). Note that the parameter (w) of ADTR displays a monotonic decreasing relationship with diversity order (q) (Fig. 1), and w of BDTR displays a valley-shape relationship (Fig. 3). It is also noted that w values listed above are the average of 32 individual subjects, and w for a specific individual subject may be negative.
It is also noted that the average PL exponent (w) of the zero-order BDTR model is equal to 0.407 for the HVM. This w-value is very close to its counterpart in the temporal scaling of alpha-diversity ADTR (w = 0.406), and both nonetheless have very different ln(c) values. The average ln(c) across all 32 subjects is 0.078 and 2.898 for BDTR and ADTR, respectively. Therefore, although the values of the scaling parameter (w) in both alpha-and beta-diversity are almost equal, the predicted diversities by the scaling models can be very different due to the significant difference in parameter c. The same magnitude of w in both alpha and beta DTR suggests that the rates of diversity accumulation at the zero-order (species-richness level or traditional STR) are nearly the same.
White et al. 61 survey revealed that the temporal scaling exponent (w) is 0.10-0.51 with a mean of w = 0.29 through the comparison of 984 community type series of eukaryotic organisms. Oliver et al. 21   which corresponds to the traditional STR, does fall in the ranges reported in the literature. Since our study, to the best of our knowledge, is the first investigation of the DTR beyond richness level (q = 0), we establish the first report on the alpha ADTR ranges for diversity order q = 1-4 as 0.005-0.073 for the HVMs. Similarly, since we are not aware of any report on the beta-diversity scaling of BDTR, we establish the zero-order BDTR of the HVM as 0.407, and the range of the BTDR for order q = 1-4 as 0.018-0.081. The possible negative exponent (w) in PL and its extended versions (PLEC & PLIEC) reveal an important difference between DAR (Ma) 49 and DTR: the spatial scaling is usually positive, but the temporal scaling can be negative. That is, diversity could decrease with time, which makes sense even intuitively.
Visualization of the DTR models. Figure 1 shows the 3D-graph of the PL exponent (w) of 32 alpha-DTR or ADTR models for 32 subjects at each diversity orders (q = 0-4), and Fig. 2 shows the beta-diversity counterpart, i.e., PL w of BDTR. Figures 3 and 4 show the PLEC counterparts, i.e., PLEC w of ADTR and PLEC w of BDTR, respectively. Figure 5 shows the 3D-graph alpha-MAD or maximal accumulation alpha diversity, predicted by the ADTR models for each of the 32 subjects at each diversity order (q = 0-4), and Fig. 6 shows the beta-diversity counterpart, i.e., beta-MAD or maximal accumulation beta diversity, predicted by the BDTR models.
Both Figs 1 and 2 show that negative PL w-values exist widely in DTR; especially with alpha diversity-time scaling (ADTR), w was drawn within the range between [-0.523 and + 0.656] (Fig. 1), a potential equilibrium around w = 0 (no change in diversity). With beta-diversity scaling (BDTR), w was drawn in the range between [−0.137 and +0.569] (Fig. 2), and the negative amplitude of accumulation (fluctuation) is smaller than the positive amplitude. Figures 1 and 2 also exhibit rather strong heterogeneities of w among 32 individual subjects, and the heterogeneity appears more significant at lower order diversities (especially for q = 0, i.e., the traditional STR).
The PLEC exponents (w) show similar pattern ranging from negative to positive, with the PLEC alpha w of ADTR ranging [0.935, 1.045], and the PLEC beta w of BDTR ranging [−0.125, 0.823]. As demonstrated previously, PLEC is more appropriate for DTR modeling since PL may fail to describe the DTR at higher diversity orders. One observation is that the range of alpha-diversity scaling exponent (w) is largely symmetrical around zero, but that of beta-diversity w is asymmetrical around zero, with a larger positive fluctuation.
The predicted maximal accumulation diversities or MADs (Figs 5 and 6) demonstrate high heterogeneity, varying significantly from individual to individual, especially at lower diversity orders. The alpha-MAD across individuals and diversity orders is 170 (Fig. 5), and the beta counterpart is 10.26 (Fig. 6). The latter predicts that the maximal beta diversity fluctuation in the 32-healthy cohort can be equivalent to 10 distinctly different communities according to the interpretation of multiplicative beta-diversity (Chao et al.) 40,41 .
The heterogeneities of the DTR parameters among individual subjects. The heterogeneities of the PL parameters [w, ln(c)] for DTR and those of the PLEC & PLIEC parameters [w, ln(c), d] were also captured quantitatively by fitting these parameters to two standard statistical distributions, i.e., the normal distribution and power law distribution [Eqn. (13)], which are often utilized to represent the typical symmetrical and asymmetrical heterogeneities-respectively. In Suppl. Table S3, we listed the distribution-fitting results, i.e., the p-values for testing the statistical distribution of each DTR parameter [w, ln(c), d] for five diversity orders (q = 0-4) and both diversity types (alpha and beta). Approximately 60% parameter cases followed the normal distribution and 97% cases followed the power law distribution. However, the predicted maximal accumulation diversities in all but one case (subject) followed the power law distribution. This suggests that, while the distributions of some parameters are symmetrical, the distributions of most parameters are highly skewed. Furthermore, the distribution of the maximal accumulation diversities among individual subjects is again highly skewed. It is noted that the statistical             (11,12)], and the full results were exhibited in Suppl. Table S2A and B. To avoid the potential bias from arbitrarily ordering the sequence of individual subjects in accruing individuals (area) and corresponding diversities, 50 runs of random samplings from the total permutations of the 32 subjects were conducted to build 50 DTAR models separately, and the average parameters from the 50 runs were adopted as the final DTAR model parameters. The average BPL model parameters from Suppl. Tables S2A and B are summarized in the following Table 3. From Table 3 and Suppl. Tables 2A and B (detailed results from the 50 runs), we report the following findings as explained below.
The performance of BPL in modeling the alpha DTAR and its parameter ranges. Regarding the alpha-diversity scaling of DTAR, the BPL [Eqns (11), (12)] fits to the DTAR data of the 32-healthy cohort highly significant with an average p-value = 0.01 across 50 runs (all but one has p < 0.001). The average temporal-scaling parameter w is (0.242, −0.010, −0.029, −0.029, −0.029) for q = 0 th -order to 4 th -order of diversities, respectively, exhibiting a decline trend with the increase of the diversity order (q) in general, becoming stabilized after order (q > 1), and negative scaling (w) beyond the zero-order diversity scaling (species richness). The average spatial-scaling parameter z is (0.461, 0.292, 0.209, 0.183, 0.171) for q = 0 th -order to 4 th -order of diversities, respectively, exhibiting a decline trend with the increase of the diversity order (q), but keeping all positive values. The approximately twice size of the spatial-scaling parameter (z) over temporal-scaling parameter (w) at the diversity order q = 0, suggests that, although both have positive scaling effects at the species-richness level (0 th -order Hill numbers), space ('area' or individual) is more significant than time in accumulating diversity (species). This certainly makes sense since it simply says that the inter-subject heterogeneity (variability) is more significant than the temporal heterogeneity (variability) of an individual subject in contributing the accumulation of species-richness across the members of a cohort. For the higher order diversities, the difference is even more categorical, with time having negative scaling effects at higher order beyond the basic species-richness level. This finding casts doubt on Preston 1 ergodic conjecture; that is, space and time at least in ecological time, do not seem to be equivalent in terms of the accumulation of diversity.
The performance of BPL in modeling the beta DTAR and its parameter ranges. Regarding the beta-diversity scaling of DTAR, the same BPL, utilized for describing the alpha-diversity DTAR as described above, also fits to the beta-diversity DTAR of the 32-healthy cohort dataset extremely significant with all p-values < 0.001 for 250 tested models. The average temporal-scaling parameter w is (−0.151, −0.100, −0.076, −0.078, −0.083) for q = 0 th -order to 4 th -order of diversities, respectively, exhibiting a mostly increasing trend with the increase of the diversity order (q) but all negative (even for species richness) scaling w. This mostly increasing trend of temporal beta-diversity scaling is also opposite with the trend demonstrated by the alpha-diversity scaling summarized above. The average spatial-scaling parameter z of DTAR model is (0.404, 0.240, 0.201, 0.205, 0.217) for q = 0 th -order to 4 th -order of beta-diversities, respectively. The scaling trend exhibits a mixed non-monotonic trend (valley-shaped) with the increase of the diversity order (q). Note that all scaling parameters are positive with the spatial-scaling parameter (z). The results from the beta-diversity DTAR further suggest that space and time can be different in terms of the accumulation of beta-diversity; that is, Preston 1 conjecture may not be true at the ecological time scale we tested. Indeed, the results in beta-diversity scaling indicate that even at species-richness level, time has a negative scaling parameter, and is different from space (area) in shaping the species accumulation. Figure 7 shows the graph of the average BPL exponents (z & w) of the alpha DTAR models across 50 runs (each run produced a DTAR model) for each of the 5 diversity orders (q = 0-4). Figure 8 shows their beta-diversity counterparts. The above-described properties such as the ranges of BPL parameters (z & w) as well as their change trends with diversity order (q) all can be visually checked in Figs 7 and 8. The distribution of DTAR parameters. Suppl. Table S4 demonstrates that the power law statistical distribution fits to the BPL parameters z, w, & ln(c) extremely significant, although the normal distribution also succeeded in 15 out of 60 cases, mostly for PLIEC or higher order beta-diversities (q = 2-4).
In conclusion, the scaling effects of both area and time can be adequately captured by the BPL, with area seemingly having a larger effect than time on diversity accumulation. Area seems to always affect diversity accumulation positively (z > 0). In contrast, beyond the species-richness level, time may have a negative effect on diversity accumulation (w may be less than 0). Especially with the beta-diversity scaling, the effect of time can be negative even for species-richness accumulation. Therefore, Preston's 1 conjecture may not be true beyond species richness at the ecological time scale. It is generally recognized that the bacteria diversity of HVM is likely to experience periodic decline and restoration due to menses (Ma et al.) 36 . Therefore, a negatively temporal-scaling exponent is certainly possible. In addition, we suggest investigating the general DTAR further when more extensive joint spatial-temporal data sets are available.

Discussion
Multiple mechanisms have been proposed to explain the STRs, including random sampling process at short time scales and ecological processes such as climatic variability, succession changes, meta-population dynamics at longer time scales ( 11,13,14,19,50,61,62 . However, direct experimental evidence to support these hypotheses is still scarce. Here, we only discuss one possible mechanism-the self-similarity or scale invariance, determined by the mathematical function of the DTR, i.e., the power law (Ma) 63 . Since the self-similarity property has been extensively studied in the context of the SAR, the following discussion aims to extend the results from the SAR to the STR and DTR.
From the basic power-law of the DTR, Eqn. (5) 14 , we can make the following derivations: The exponent w is here the proportionality constant in the rate-equation form of the power law. The left side is the diversity increase per increase in time (unit), i.e., the tangent of the PL curve [Eqn. (5)], the right side is the w multiplied by the diversity per unit time, termed temporal diversity density. The above equation can be rewritten as the following equation: which further reveals that w is the ratio of diversity accumulation rate to time-increase rate, or the relative diversity-accumulate rate.
The second parameter (c) is not independent of the temporal sampling unit used, which makes the comparison of c between different case studies difficult. According to the self-similarity hypothesis, c is the diversity in one unit of time, that is, D 0 = c when T = 1. It is the diversity within one unit of time, but not per unit of time.
Self-similarity is also known as scale-invariant, which refers to the following mathematical property of the power law: w w that is, scaling the argument T (time) by a constant factor α, is equivalent to scaling its function proportionally by a constant factor α . w Therefore, all power laws with a particular scaling exponent w are equivalent up to constant factors because each is a scaled version of the others. The scale-invariance is also responsible for the linear relationship after log-transformation of power law [ = + D c w T ln( ) ln( ) ln( )], and the resulted straight line on