Ecology

Study site
The study site (N 38° 49′, E 115° 26′) is located in the Guanzhuang Village in Baoding City of Hebei Province, China, in the humid temperate and monsoon climatic zone with the average annual air temperature of 13 °C, annual rainfall of 500 mm, and frost-free period is 210 days. Although the experiment was a one-year, the distribution of precipitation (488.50 mm) and temperature (13.45 °C) during the experimental period (2014–2015) were close to the the latest 10-year averages (2005–2015) (500.19 mm and 13.61 °C) (Fig. 1). Determine the basic nutrient indexes of the 0–20 cm surface soil in the test plot. The soil type is silty loam, consists with 22.55% sand, 71.09% silt and 6.36% clay. Analysis of soil basic characteristics showed that it has a pH of 8.3, and its content of organic matter, total N, available phosphorus (P) and available potassium (K) was 11.27 g kg−1, 1.47 g kg−1, 25.49 mg kg−1 and 127.43 mg kg−1, respectively.
Figure 1

Monthly precipitation and average temperature during the experimental year (2014–2015) and the mean values in the last ten years (2005–2015) in the test area. Data was from meteorological station.

Full size image

Experiment materials
The planting mode of the experimental site was a winter wheat-summer maize rotation, the winter wheat variety was ‘Jinnong 6’ that verage thousand weight was 47.6 g. The summer maize variety was ‘Zhengdan 958′ that average thousand weight 330 g.
Test fertilizers include inorganic fertilizer, organic fertilizer, soil conditioner, compound bateria, amino acid liquid fertilizer and nutrient agent. Inorganic N, P and K in the tested fertilizers were provided by urea (N 46%), superphosphate (P2O5,16%) and potassium chloride (K2O, 54%), respectively, as well as in the form of zinc and humic acid urea which is mainly a combination of N with humic acid (N 46% and HA 1.2%). The organic fertilizer used in the experiment was mainly decomposed chicken manure. But the N content in the chicken manure is 1.32% in wheat season and 4.48% in maize sason. Soil conditioner mainly containing calcium (Ca) and magnesium (Mn), Compound bacteria could fix N potentially, promote root growth, decompose cellulose lignin and thus rapidly to degrade. The number of living bacteria reached 2 billion per gram. Amino acid liquid fertilizer and nutrient agent sprayed according to crop growth to provide the required amino acids and trace elements for plant growth.
Experiment design
Field experiment consisted of five treatments with 3 replicates. The experiment uses a completely randomized block setting, the plot size was 79.2 m2 (13.2 m × 6 m). Before the experiment, no crops were planted in the area and it was idle for more than one year. The five treatments were: CK (zero N), FN (farmers’ traditional inorganic N rate, through mass surveys on actual production), RN (recommended inorganic-N rate, according to the experimental results of many scholars, combined with the local soil N supply, crop straw returning in the previous season as well as wheat or maize N demand for target yield in the current season)22,23, HAN (zinc and humic acid urea, the N supply same as RN), RN40% + HOM (40% inorganic N rate of RN (RN40%) with homemade organic matters (HOM). HOM was an organic control measure, it including organic fertilizer, soil conditioner compound bacteria, amino acid liquid fertilizer and nutrient agent. these constituents and amount according to Shu et al24 (Table 1).
Table 1 Rates of pure N and organic matters in different fertilization treatments.
Full size table

For wheat, N fertilizer was broadcast for ratio of 4:3 (basal to topdressing) in RN40% + HOM, whereas for the rest N treatments the ratio was 1:1. During maize planting, N ratio (basal to topdressing) for all treatments was 2:3.The N, P and K fertilizers for wheat were applied in the form of urea, single superphosphate and potassium chloride, respectively. The amount of N fertilizer applied in different treatments of different crops is different, specific application amount reference Table 1. Except for treatment RN40% + HOM, all treatments have the same amount of single superphosphate and potassium chloride. Single superphosphate (120 kg P2O5 ha−1)and potassium chloride (150 kg K2O ha−1)were used in winter wheat season. For maize, single superphosphate (90 kg P2O5 ha−1) and potassium chloride (150 kg K2O ha−1) were used. P and K fertilizers were applied once before sowing. For the doses of P and K in RN40% + HOM brought by organic fertilizer were firstly assessed (48.4 kg P2O5 ha−1 and 149.3 kg K2O ha−1 for winter wheat; 87.2 P2O5 ha−1 and 19.1 kg K2O ha−1 for maize), remaining amounts were supplemented with chemical P and K fertilizers.
Wheat at a rate of 187.5 kg ha−1 with a row space of 15 cm, was sown on 12 October 2014 and harvested at 7 June 2015. Then, at the same wheat plot, Maize of 37.5 kg ha−1 with a row space of 57 cm, was sown on 18 June 2015 and harvested at 5 October 2015.
N2O sampling and measurements
N2O gas was collected using a closed static chamber from sowing to harvest of wheat and maize25. The sampling box was divided into two parts and made by PVC material: a box body and a base. The upper part of the box body was provided with a gas sampling port sealed with a rubber plug, and a thermometer probe was arranged inside the box body to monitor the soil surface temperature. The box body is 15 cm high and the bottom diameter is 25 cm. The base was ring shaped, and buried into the soil. Gas collection was performed from 9:00 to10:00 am. A 30 mL of air sample was collected at 0, 8, 16 min after closure26. The air samples were taken once at an interval of 7 days in general and subsequently continuous 5 days following fertilization or precipitation. Continue to collect gas samples from the beginning of the experiment. No gas samples are collected during the freezing of wheat field soil in February and March every year. At the same time, the air temperature was measured by a thermometer and the soil moisture in the 0–5 cm depth was measured by soil moisture tester (TK3-BASIC). Gas concentrations were analyzed by using a gas chromatography (Agilent 7890 A, USA), fitted with a 4 mm by 3 m stainless steel column packed with Porapack Q and N2 was used as the carrier gas. The column and the detector temperatures were set at 70 °C and 300 °C, respectively. The standard N2O was supplied from National Center of Standard Measurement.
N2O flux was calculated using the following equation (Wang et al.27).

$${text{F}} = rho times {text{H}} times T_{0} frac{{left( {c_{2} /T_{2} – c_{1} /T_{1} } right)}}{Delta t}$$
(1)

where F is N2O emission flux, ρ = m/22.414, ρ is the density of gas in airtight box, m is molecular weight, H is the height of the static chamber, T0 is 273 K, c1 and c2 are the gas concentration in time of t1and t2, respectively, T1and T2 are gas temperatures, ∆t = t2 − t1 ,where t2 and t1are times.
Cumlative N2O emissions were from the growth season was calculated by the equation:

$${text{T = }}sum {left[ {{{left( {{text{F}}_{{text{i + 1}}} {text{ + F}}_{{text{i}}} } right)} mathord{left/ {vphantom {{left( {{text{F}}_{{text{i + 1}}} {text{ + F}}_{{text{i}}} } right)} {2}}} right. kern-nulldelimiterspace} {2}}} right]} times left( {{text{D}}_{{text{i + 1}}} – {text{D}}_{{text{i}}} } right) times {{{24}} mathord{left/ {vphantom {{{24}} {{1000} times {667} times {15} times {10}^{{ – {6}}} }}} right. kern-nulldelimiterspace} {{1000} times {667} times {15} times {10}^{{ – {6}}} }}$$
(2)

where T is the total amount of N2O emissions from the growth stage (kg N ha−1), Fi and Fi+1 denote the N2O flux of the i and i + 1 sub-sampling (μg N m−2 h−1); Di and Di+1 represent sampling days (d)26.
N2O emission coefficient (EF) was estimated with equation28:

$$EF(% ) = left[ {left( {{text{Cumulative}};{text{N}}_{{2}} {text{O}};{text{emissions}};{text{from}};{text{fertilized}};{text{plots}} – {text{control}};{text{plots}}} right)/{text{N}};{text{fertilizer}};{text{rate}}} right] times 100$$
(3)

Soil sampling and measurements
At wheat and maize maturity, soil samples were collected from depths of 0–20, 20–40 and 40–60 cm with a hand probe from three places in central rows of each plot and mixed together. Fresh soil samples were sieved through a 2 mm, extracted with 1 mol L−1 KCl and a soil-solution ratio of 1:10, and analyzed for inorganic N (mainly including NH4+–N and NO3−–N) contents with continuous flow analysis technique(AA3-HR, Germany)22. Soil moisture and density of each soil layer were measured simultaneously, and the soil residual N in 0–60 cm was calculated. Other soil sample was air-dried and sieved, organic matter and total N content were measured by agrochemical analysis method29.
Plant harvest
For wheat, plants with double rows (1 m length) in each plot were harvested and 20 spikes were selected to count the numbers of effective spikes. All the harvested plant samples were separated into straw (including stem, leaves and remaining of ears) and grains, and the grain yield was calculated to 12.5% moisture content (PM-8188, Japan). Three samples were chosen from each plot and weighted to get the average 1000-grain weight.
For maize, two representative plants in each plot were harvested and separated into straw (including stems, leaves, tassels, husks, cobs) and grains in the central rows. Moreover, 20 ears were continuously selected to thresh and measured grain yield. Grain yield was calculated to 14% moisture content (PM-8188, Japan).
All harvested wheat and maize samples were dried, weighed, ground into powder to measure the total N content using H2SO4–H2O2 Kjeldahl digestion method29.
N balance and N efficiencies
Total N input was comprised of N fertilizer, the initial inorganic N in soil before planting (including both NO3−–N and NH4+–N), pre-crop N straw return (no straw was returned when sowing wheat and the N uptake in maize stage was calculated from pre-wheat straw), N deposition from dry and wet atmosphere and mineralized N in soil. Atmospheric N deposition was derived from Research result by Liu et al.30. N output was comprised of crop uptake, post-harvest residual soil N and apparent N loss. This study calculated soil N to a depth of 0–60 cm. Mineralized soil N, apparent N loss, Nitrogen production effiency (NPE) , Nitrogen agronomic effiency (NAE) and Nitrogen use efficiency (NUE) were calculated as follows:

$$begin{aligned} {text{Mineralized}};{text{N}}left( {{text{kg}},{text{ha}}^{ – 1} } right) & = {text{Crop}};{text{N}};{text{uptake}};{text{in}};{text{CK}} + {text{Post – harvest}};{text{residual}};{text{soil}};{text{N}};{text{in}};{text{CK}}{-}{text{Pre – planting}};{text{soil}};{text{N}};{text{in}};{text{CK}} \ & quad {-}{text{N}};{text{deposition}};{text{from}};{text{atmosphere}};{text{in}};{text{CK}} – {text{Pre – crop}};{text{straw}};{text{return}};{text{N}};{text{in}};{text{CK}} \ end{aligned}$$
(4)

$${text{Apparent}};{text{N}};{text{loss}}left( {{text{kg}},{text{ha}}^{ – 1} } right) = {text{Total}};{text{N}};{text{input}} – {text{crop}};{text{N}};{text{Uptake}} – {text{post – harvest}};{text{residual}};{text{soil}};{text{N}}$$
(5)

$${text{NPE}}left( {{text{kg}},{text{kg}}^{ – 1} } right) = {text{Plant}};{text{yield}}/{text{N}};{text{fertlizer}};{text{rate}}$$
(6)

$${text{NAE}}left( {{text{kg}},{text{kg}}^{ – 1} } right) = left[ {left( {{text{Plant}};{text{yield}};{text{with}};{text{N}};{text{application}} – {text{plant}};{text{yield}};{text{without}};{text{N}};{text{fertlizer}}} right)/{text{N}};{text{fertlizer}};{text{rate}}} right]$$
(7)

$${text{NUE}}left( % right) = left[ {left( {{text{N}};{text{content}};{text{in}};{text{plant}};{text{with}};{text{N}};{text{fertilizer}}{-}{text{N}};{text{content}};{text{in}};{text{plant}};{text{without}};{text{N}};{text{fertilizer}}} right)/{text{N}};{text{fertlizer}};{text{rate}}} right] times 100$$
(8)

Net income analyses
Prices of fertilizers and grains as well as other costs in Chinese Yuan (RMB: 1 USD = 6.71 RMB in the experiment year) were based on local prices. Net income was calculated by the equation:

$${text{Net}};{text{income}} = {text{Output}};{text{value}}{-}{text{fertilizer}};{text{cost}}{-}{text{other}};{text{field}};{text{management}};{text{costs}}$$
(9)

$${text{Output}};{text{value}} = {text{Grain}};{text{yield}} times {text{grain}};{text{price}}$$
(10)

where Fertilizer costs were composed of the prices of inorganic N (3.9 RMB kg−1), P2O5 (5.65 RMB kg−1), K2O (6.5 RMB kg−1), pure N in zinc and humic acid urea (5.0 RMB kg−1), decomposed chicken manure (0.5 RMB kg−1), soil conditioner (2.8 RMB kg−1), compound bacteria, amino acid liquid fertilizer and nutrient agent together (30 RMB kg−1). Other field management costs included seed, labor for fertilization, irrigation, mechanical sowing, etc. Grain prices of wheat and maize during the experiment were 2.2 and 1.8 RMB kg−1, respectively.
Statistical analysis
This research adopted SPSS Statistics 20.0 software (SPSS Inc., Chicago, IL, USA) to date analysis. Through least significant differences (LSD) method, the statistically significant differences were calculated. The differences level was prominent when P  More

Jaw shape variation and diet in small mammals
Using 2D geometric morphometrics (Fig. 2a), we found that jaw shape is a good proxy for diet among small extant mammals. In Fig. 3, taxa with negative PC1 scores have shorter jaws, and taxa with positive PC1 scores have longer jaws; taxa with positive PC2 scores have taller ascending rami and taxa with negative PC2 scores have shorter ascending rami. Among extant mammals, most dietary categories (excluding omnivores) can be distinguished along PC1 (Fig. 3a): herbivores plot at the negative end of PC1, insectivores towards the positive end, and carnivores in between. These categories are also statistically different from each other (Table 2), showing that jaw shape can distinguish between most major dietary types. However, our data cannot distinguish between carnivores and omnivores.
Fig. 2: Data acquired from the jaws of Mesozoic and extant small mammals.

a Jaw landmarking regimen used in this study. Modified from ref. 12. In orange: six fixed landmarks; in blue: 58 sliding semi landmarks. b Moment arm measurements taken in this study. Modified from ref. 19.

Full size image

Fig. 3: Scatter plots of the principal component analysis (PCA) results (PC1 vs. PC2).

a Extant taxa, b extinct taxa. Convex hulls shown for extant insectivores (yellow), carnivores (red), omnivores (purple) and herbivores (blue). Icon colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). See Table 1 for taxon names.

Full size image

Table 2 Summary of the Procrustes ANOVA (Type II, Conditional SS) performed for jaw shape data as a function of dietary group.
Full size table

Data on the jaw shape of Mesozoic mammals were projected onto the extant taxa morphospace (Fig. 3b). In order to determine whether jaw shape could be used as a dietary proxy in Mesozoic mammals, we obtained previous independent determinations of likely diets, which variously employed dental morphology, tooth wear facets and body size (e.g., see refs. 1,7,12,14,24,25,26,27,28,29,30,31,32). We saw a very good correspondence between previous proposed diets for Mesozoic mammals and their position on the morphospace. See Supplementary Fig. 6 for a principal components analysis scatter plot which includes multituberculates and haramiyidans; these taxa were excluded from our study because the vast majority of them have jaw shapes dissimilar to the other extinct and extant mammals in our sample (i.e., allotherians have shorter jaws and thus more negative PC1 scores).
Stem mammals
Most stem mammals plot within the morphospace of extant insectivores and have positive PC1 scores. One exception is Sinoconodon (taxon #2, Fig. 3), which plots within the morphospace of extant carnivores; Sinoconodon is considered a carnivore based on dental morphology5. Haramiyavia (#1) is thought to have been a plant-dominated omnivore23 based on dental morphology, but here it plots within the morphospace of extant insectivores. Both morganucodontans in this study, Morganucodon (#3) and Dinnetherium (#4), have similar PC1 scores to extant insectivores, echoing the findings of Gill et al.14.
Molar morphology indicates omnivorous or faunivorous diets for docodontans; here they mostly plot within the morphospace of extant insectivores, with the exception of Haldanodon (#6) and Docofossor (#7). Agilodocodon (#9) was previously considered a plant-dominated omnivore, with exudativorous dental features which indicated a diet mainly composed of plant sap33; more recently, Wible and Burrows34 challenged this hypothesis and suggested that the teeth of Agilodocodon most closely resemble those of extant insectivores. Here, Agilodocodon plots firmly within the morphospace of extant insectivores, close to the insectivorous dusky antechinus (Antechinus swainsonii, #61) and the elephant shrews (Elephantulus rufescens [#114] and E. brachyrhynchus [#115]), which are insect-dominated omnivores.
According to Ji et al.28 the swimming docodontan, Castorocauda (#5), has dental features indicative of feeding on aquatic invertebrates and small vertebrates, like fish. Castorocauda is often depicted as being carnivorous and, particularly, piscivorous7,28,33. The jaw shape of Castorocauda is similar to that of modern day insectivores; this docodontan might have been feeding on “soft” aquatic invertebrates (Fig. 3). The other Mesozoic mammal in our sample proposed to have been semi-aquatic, Teinolophos (#13), plots in a similar area of the morphospace to Castorocauda. Our extant sample also includes a semi-aquatic carnivore, the water opossum (Chironectes minimus, #69), which plots in the middle of the carnivore morphospace, far away from Castorocauda and Teinolophos.
Docofossor (#7) has skeletal features indicative of a fossorial lifestyle and a dentition similar to those of extant mammals foraging underground, such as moles, solenodons, and tenrecs35. This docodontan has previously been considered an insectivore7. Here, Docofossor plots within the morphospace of extant carnivores; however, it plots close to the burrowing Hispaniolan solenodon (Solenodon paradoxus, #109), which has an insectivorous diet. Among the extant insectivores in our sample, the burrowing vermivores (e.g., the hairy-tailed mole, Parascalops breweri [#108], and the Hispaniolan solenodon) have more negative PC1 scores than other insectivores (similar to that of Docofossor), and their PC1 values are more similar to those of carnivores.
The dental morphology of Haldanodon (#6) is indicative of an insectivorous diet. Here, it plots within the carnivore morphospace (very near extant herbivores), because of its tall coronoid process and comparatively shorter jaw. Docodon (#8) likely ate insects and other small invertebrates27 and, based on its diminutive size36, Microdocodon (#10) was probably insectivorous. Both of these docodontans plot within the insectivore morphospace.
Non-therian crown mammals
The jaw shape of non-therian crown mammals varies widely, plotting mostly within the morphospace of insectivores and carnivores. Fruitafossor (#11), a fossorial mammal with teeth similar to extant armadillos, has been considered an omnivore eating insects, small invertebrates and some plants26. Here, it plots within the insectivore morphospace, closely to the insectivorous and fossorial hairy-tailed mole (Parascalops breweri, #108), and shares similar PC1 scores with other fossorial taxa, such as Docofossor (#7) and the Hispaniolan solenodon (#109).
Extant monotremes eat insects and other small invertebrates. It has been proposed that the Early Cretaceous monotreme Teinolophos (#13) had a semiaquatic lifestyle (on the basis of its enlarged mandibular canal37) and ate in a similar manner to the insectivorous Kuehneotherium38. Here Teinolophos, and the australosphenidan Henosferus (#12), have PC1 scores similar to insectivores and omnivores.
The eutriconodontans are a very diverse group of insectivores and carnivores which had a wide range of body sizes, including some of the largest Mesozoic mammals known1. Here all eutriconodontans fall within or very close to the extant carnivore morphospace. In particular, Triconodon (#16) and Argentoconodon (#19) plot within the carnivore morphospace, Trioracodon (#17) and Volaticotherium (#18) plot between the carnivore and insectivore morphospaces, and Yanoconodon (#15) plots within the insectivore morphospace. Both gobiconodontids, Gobiconodon (#20) and Repenomamus (#21), have more negative PC1 scores and plot closer to the herbivore morphospace, but still remain within or close to the carnivore morphospace. Triconodon, Trioracodon, Gobiconodon, and Repenomamus are all considered carnivores based on craniodental morphology and body size1,7,31; additionally, there is direct evidence for the carnivorous diet of Repenomamus from fossilized stomach contents4. Yanoconodon and Volaticotherium are considered insectivores7.
“Symmetrodontans” like Spalacotherium (#22), Zhangheotherium (#24) and Maotherium (#25) have often been considered insectivores based on their craniodental morphology1,7 (note “symmetrodontans” likely do not represent a monophyletic group, but are often grouped together based on their tooth morphology1). Here, all “symmetrodontans” plot within the insectivore morphospace. Dryolestids are also commonly considered insectivorous1,29. Here, Crusafontia (#26) plots between the morphospace of extant carnivores and insectivores, while Amblotherium (#27) plots within the insectivore morphospace.
Vincelestes (#29) has previously been considered a carnivore on the basis of jaw shape12. Here, it plots near the morphospaces of both omnivores and herbivores. Bonaparte24 considered the incisor wear of Vincelestes reminiscent of Cenozoic carnivores, and Rougier25 considered its jaw morphology indicative of a forceful bite enabling the incorporation of tough plant matter into a primarily carnivorous/insectivorous diet.
Therian crown mammals
Extant marsupials have a large diversity of diets, including herbivory, but the extinct metatherians in our sample are considered to have been limited in diet to insectivory and carnivory (note that there are some putatively herbivorous/omnivorous extinct metatherians, like Glasbius and polydolopimorphians39,40). Their jaw shape is very similar to that of extant carnivores and insectivores (Fig. 3). Dental morphology indicates that Eodelphis (#32) and Deltatheridium (#30) were carnivorous, Didelphodon (#31) durophagous or molluscivorous31,32, and Alphadon (#33) is considered to have been insectivorous, on the basis of its jaw shape and body size12. Dental microwear indicates a broad diet consisting of vertebrates, plants, and hard-shelled invertebrates for Didelphodon; biomechanical analyses of its skull and jaw points towards a durophagous diet15,16. Biomechanical analyses of the resistance to bending and torsion of Eodelphis jaws, points to a durophagous diet in Eodelphis cutleri and non-durophagous faunivory for Eodelphis browni16. Here, Eodelphis, Deltatheridium and Didelphodon plot closely to the extant carnivores, while Alphadon plots closely to the extant insectivores.
Extant placentals also have a wide range of diets, but many of the extinct eutherians in this study (i.e., Sinodelphys [#34], Juramaia [#35], Eomaia [#36], Kennalestes [#40], Barunlestes [#44], and Kulbeckia [#43]) are considered insectivorous7,12. Here, we corroborate this hypothesis (Fig. 3): all extinct eutherians plot within the insectivore morphospace, with the exception of Asioryctes (#38) which plots in the insectivore/carnivore morphospace, and Juramaia and Sinodelphys, which plot just outside the insectivore morphospace.
Using jaw shape to infer diet in Mesozoic mammals
We performed a phylogenetic flexible discriminant analysis (phylo FDA) following Motani and Schmitz41 to determine the posterior probability of the Mesozoic taxa belonging to one of three dietary categories: insectivore, carnivore, or herbivore (we omitted omnivores as they are not well discriminated in Fig. 3). We used the first seven PC scores (of the PCA of Procrustes coordinates of jaw shape), which together accounted for 81.39% of the variance. The results of the analysis can be seen in Fig. 4 and the posterior probability values can be seen in Supplementary Data 1. We used the extant taxa of known diets as the training dataset for the discriminant analysis: these taxa were classified correctly 89.19% of the time. For the most part, we see a good separation between dietary groups among extant mammals (Fig. 4a), with some exceptions: the primarily herbivorous olingo (Bassaricyon gabbii, #94) plots with the carnivores (although mainly frugivorous, it can consume small vertebrates), and a couple of insectivores plot very near the carnivores (i.e., the little brown bat [Myotis lucifugus, #104] and the Hispaniolan solenodon [Solenodon paradoxus, #109]). These three taxa, alongside the carnivorous greater bulldog bat (Noctilio leporinus, #101), were the only extant taxa misclassified by the discriminant analysis.
Fig. 4: Phylogenetic flexible discriminant analysis results, showing discriminant axis 1 (DA1) and two (DA2), of all taxa in this study.

Extinct taxa are color coded based on their posterior probability of belonging to one of the established dietary categories. Convex hulls show the position of the extant taxa in the plot and are color coded based on their dietary categories.

Full size image

The Mesozoic mammals included in our sample have largely been considered faunivorous and the results of the phylo FDA (Fig. 4b) corroborate this hypothesis. The majority of them are classified as insectivorous, including most stem mammals, australophenidans, “symmetrodontans” and eutherians, among others. Among the eutriconodontans, Argentoconodon, Gobiconodon, Repenomamus, and Trioracodon, are classified as carnivores, Triconodon and Yanoconodon are classified as insectivores, but with moderate support (posterior probabilities: 48% and 52%, respectively), and Phascolotherium and Volaticotherium are more confidently classified as insectivores (posterior probabilities: 60% and 73%, respectively). Among the metatherians, Didelphodon and Eodelphis are classified as carnivores, while Alphadon and Deltatheridium are classified as insectivores with moderate support (posterior probabilities: 54% and 52%, respectively). The stem mammals, Haramiyavia, Sinoconodon, and Docofossor are all confidently classified as carnivores (posterior probabilities over 80%), and the crown mammals Crusafontia and Kennalestes are also classified as carnivores, but with moderate support (posterior probabilities: 54% and 52%, respectively). Two taxa in the analysis are classified as herbivores, because of their relatively tall ascending rami: Vincelestes (#29) and Haldanodon (#6). The dental morphology of Vincelestes points to a primarily faunivorous diet24, but it has been previously noted that its jaw morphology is indicative of a forceful bite; Rougier25 suggested that this jaw morphology might have enabled Vincelestes to incorporate tough plant matter into its diet, but it might also be indicative of durophagy. The dental morphology27 and body size of Haldanodon point towards an insectivorous diet; in this analysis, the posterior probability of Haldanodon being a herbivore is not high (only 40.3%). The evidence thus far suggests Haldanodon had a faunivorous diet; its jaw morphology might be indicative of the incorporation of tougher food sources into its diet.
Mechanical advantage of the jaws of small mammals
We obtained mechanical advantage (MA) data to test whether extant mammals of different dietary groups have distinct MA values (Table 3). The mechanical advantage measurements were standardized across all jaws to account for differences in jaw morphology (e.g., presence or absence of the angular process) (Fig. 2b); the outlever was measured at the anterior end of the jaw and at the first lower molar (m1). When measuring mechanical advantage at the jaw tip and considering extant taxa only, we find statistically significant differences in the mechanical advantage of the masseter (MAM) values in all pairwise dietary combinations except for carnivore-insectivore (Table 3). The mechanical advantage of the temporalis (MAT) is statistically distinct only between herbivores and insectivores, and carnivores and insectivores (Table 3). Herbivores and carnivores do not have statistically distinct MAT values. This may differ in a sample of larger ( > 5 kg) therians. When measuring the outlever at the m1, we find statistically significant differences in all pairwise comparisons of MAM between dietary groups, except for herbivore–omnivore and carnivore–insectivore. When considering MAT, we only find significant differences between omnivores and carnivores, insectivores and herbivores, and insectivores and carnivores.
Table 3 Pairwise p values (uncorrected significance) of one-way PERMANOVAs of the mechanical advantage values of the masseter (MAM) and temporalis (MAT) obtained in this study on extant taxa of known dietary preferences only (permutation N = 9999).
Full size table

Figure 5 shows the mechanical advantage of the masseter (left) and temporalis (right), measured at the jaw tip, in a phylogenetic context (see also Supplementary Fig. 7 for individual taxon names). Phylogeny appears to have a large influence on the mechanical advantage and diet of the jaws of small mammals. Most Mesozoic taxa have low (blue) to intermediate (green) MAM values. Most stem mammals have intermediate (green) to high (red) MAM values and non-therian crown mammals have low MAM values, with the exception of Fruitafossor and Vincelestes (the latter has the highest MAM value of all taxa, both extinct and extant). Most eutherians, both extinct and extant, have intermediate to low MAM values, with the exception of the relatively high values (yellow to orange) seen in elephant shrews (order Macroscelidea) and the four-toed hedgehog (order Eulipotyphla, Atelerix albiventris). Some members of the orders Carnivora (including canids and euplerids) and Afrosoricida have some of the lowest MAM values. Metatherians have MAM values ranging from low to intermediate (in the orders Dasyuromorphia and Didelphimorphia, as well as in the Mesozoic metatherians) to some of the highest in the order Diprotodontia (e.g., the sugar glider [Petaurus breviceps], the woylie [Bettongia penicillata], the cuscus [Phalanger orientalis]).
Fig. 5: Mechanical advantage values of the masseter and temporalis when biting at the jaw tip visualized in the context of the phylogeny used in this study.

See Supplementary Fig. 7 for individual taxon names.

Full size image

Most taxa have intermediate MAT values (Fig. 5 and Supplementary Fig. 7). Very low MAT values are seen in the extinct non-therian crown mammals Teinolophos and Zhangheotherium and a few extant taxa, including marsupials like the Western barred bandicoot (Perameles bougainville) and the numbat (Myrmecobius fasciatus), and placentals such as the striped treeshrew (Tupaia dorsalis) and the short-snouted elephant shrew (Elephantulus brachyrhynchus). The highest MAT values belong to members of the order Carnivora, including skunks (Mephitis macroura and Conepatus humboldtii), the least weasel (Mustela nivalis) and the tayra (Eira barbara). Some diprotodontians like the common ringtail possum (Pseudocheirus peregrinus) and the sugar glider (Petaurus breviceps) also have relatively high MAT values. Some extinct taxa also have relatively high MAT values, including the stem mammal Docofossor, and the non-therian crown mammals, Triconodon and Vincelestes.
Figures 6 and 7 present a visualisation of the mechanical advantage of the masseter and the temporalis (x axis, outlever measured at the jaw tip) and the PC1 scores of Fig. 3 (y axis, mainly describes jaw length) because, as previously mentioned, this is the axis in which dietary categories among extant mammals are best discriminated. In the y axis of Figs. 6a and 7a, herbivores have short jaws, carnivores have short to intermediate-length jaws and insectivores have intermediate-length to long jaws. In Fig. 6a, insectivores and carnivores have low mechanical advantage values of the masseter (i.e., when biting: less forcefulness, more speed), and herbivores have higher mechanical advantage values (i.e., when biting: more forcefulness, less speed). In Fig. 7a, insectivores have lower mechanical advantage values of the temporalis, while carnivores and herbivores have higher mechanical advantage values. Note that most carnivores have intermediate MAT values, but some mustelids (i.e., the least weasel [Mustela nivalis, #99], the American badger [Taxidea taxus, #96], and the North American river otter [Lontra canadensis, #98]), have the highest MAT values among extant mammals. Also note that, among insectivores, those with the highest MAT values are burrowing vermivores (i.e., the short-tailed shrew tenrec [Microgale brevicaudata, #111], the hairy-tailed mole [Parascalops breweri, #108], and the Hispaniolan solenodon [Solenodon paradoxus, #109]). By using a combination of their MAM and MAT values (as well as their jaw length), we can distinguish dietary categories among extant mammals. We decided to omit omnivores from these figures because, as seen in Fig. 3, they cannot be distinguished from other dietary groups on the basis of jaw shape.
Fig. 6: Scatter plot of the mechanical advantage of the masseter (x axis) vs. PC1 scores from Fig. 3 (y axis), which mainly describes jaw length.

a Extant taxa, b extinct taxa. Colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). Ovals indicate where extant taxa of known dietary categories plot, as in part a.

Full size image

Fig. 7: Scatter plot of the mechanical advantage of the temporalis (x axis) vs. PC1 scores from Fig. 3 (y axis), which mainly describes jaw length.

a Extant taxa, b extinct taxa. Colors indicate known dietary categories of extant mammals and suggested dietary categories for Mesozoic mammals (obtained from the literature). Ovals indicate where extant taxa of known dietary categories plot, as in part a.

Full size image

We also obtained additional mechanical advantage measurements, in which the outlever was measured at the first lower molar (m1), rather than the jaw tip (Supplementary Figs. 8, 10, 11, and 13). We made this alternative measurement because Grossnickle17 found that the length of the posterior portion of the jaw (measured from the jaw joint to the m1) is a strong predictor of diet in mammals. Compared to the mechanical advantage (MA) measurements at the jaw tip (Figs. 6a and 7a), we see a less clear distinction between dietary groups among extant mammals. There is considerable overlap between dietary groups in Supplementary Fig. 10 (jaw length~MAM). In Supplementary Fig. 11 (jaw length~MAT), there is a better separation between dietary groups.
Based on previous likely determinations of diet of Mesozoic mammals (see Supplementary Data 1 for the full list of sources), most taxa plot where it is “expected” of them, with some exceptions (Figs. 6b and 7b): 1) about half of the stem mammals (i.e., Haramiyavia, Sinoconodon, Morganucodon, Haldanodon, and Docofossor), most of which are thought to have been faunivorous, have higher MAM values than modern insectivores and carnivores, and 2) the docodontan Castorocauda has MAM and MAT values consistent with an insectivorous diet, as opposed to the carnivorous diet proposed for this taxon7,28,33. Most Mesozoic mammals have mechanical advantage values similar to modern insectivores, a few taxa are similar to carnivores (e.g., Sinoconodon, Triconodon, Trioracodon, Argentoconodon, Gobiconodon, Repenomamus, Deltatheridium, Didelphodon, and Eodelphis), and some are more similar to herbivores (e.g., Vincelestes and Fruitafossor). More

Model equations
We introduce a spatiotemporal model to describe the bacteria-phage interaction in the upper part of the soil with the depth H (we consider (H=1) m) in a typical agricultural field. Here we consider a 1D model where all abiotic and biotic components depend on time t and vertical coordinate h. The biotic component of the model consists of 4 compartments: phage-free bacteria (S) susceptible to infection by the phage, bacteria infected by the phage in its lysogenic ((I_1)) and lytic ((I_2)) states, and free phages (P). The total density of the host bacterial populations N is defined as (N = S + I_1 + I_2). The schematic diagram illustrating bacteria-phage interactions is similar to that of Egilmez and co-authors16. The local species interactions are described based on the classical modelling approach6,19. Our spatiotemporal model is of reaction-diffusion type and is described by the following equations

$$begin{aligned} begin{aligned} frac{partial S(t,h)}{partial t}&= D_b frac{partial ^2 S(t,h)}{partial h^2} +alpha (T) S(t,h) Big [1-frac{N(t,h)}{C(h)}Big ] – K_S S(t,h)P(t,h), \ frac{partial I_1(t,h)}{partial t}&= D_b frac{partial ^2 I_1(t,h)}{partial h^2} + {overline{alpha }}(T) I_1(t,h) Big [1- frac{N(t,h)}{C(h)}Big ] + K_1(T) S(t,h) P(t,h) – lambda _1(T) I_1(t,h), \ frac{partial I_2(t,h)}{partial t}&= D_b frac{partial ^2 I_2(t,h)}{partial h^2} + K_2(T) S(t,h) P(t,h) + lambda _1(T) I_1(t,h) – lambda _2 I_2(t,h), \ frac{partial P(t,h)}{partial t}&= D_P frac{partial ^2 P(t,h)}{partial h^2} -K N(t,h) P(t,h) – mu P(t,h) + b lambda _2 I_2(t,h). end{aligned} end{aligned}$$
(1)

In the above model, we parameterise the growth of susceptible bacteria via a standard logistic growth function6, where (alpha) is the maximal per capita growth rate and C is the carrying capacity of the environment; we assume that C(h) varies with depth. Infection of S by phages P at low temperatures results in lysogeny which is described by a mass action term (K_s S(t,h) P(t,h)). The growth of lysogenic bacteria (I_1) is described by a logistic function as in the case of S; however, with a different maximal growth rate ({overline{alpha }} (T)) as detailed in the next subsection. At high temperatures, the transition from the lysogenic to the lytic cycle of infection occurs: this is described by the term (lambda _1 (T) I_1(t,h)). Infection by the phage via the lytic cycle is modelled by the term (K_2 (T)S(t)P(t)). The death rate of infected bacteria due to lysis is modelled by (lambda _2 (T) I_2). The lysis of a bacterium results in the release of b new phages, the the burst size6. In the equation for P, KN(t)P(t) stands for the loss of phage due to binding to any type of bacteria (for simplicity, we assume that there is no saturation in binding). Finally, (mu P(t,h)) is the natural mortality or deactivation of phages.
According to this framework, the vertical displacement of the phage and bacteria are modelled by a diffusion term (first term in each equation), where (D_b) and (D_P) are the diffusion coefficients of bacteria and phage, respectively. The variation of the temperature T across the soil is described by the heat equation

$$begin{aligned} frac{partial T(t,h)}{partial t} = D_h frac{partial ^2 T(t,h)}{partial h^2}, end{aligned}$$
(2)

where (D_h) is the diffusion coefficient of heat transfer (see more detail in the next section). Models (1)–(2) should be supplied with appropriate boundary conditions. We assume that the model has the zero-flux boundary condition for all biotic components (bacteria and phage) at (h=0) and (h=H). For the temperature, we consider Dirichlet boundary conditions such that (T(t,0)= T_s (t)) and (T(t,H)= T_H), where (T_s (t)) is the surface temperature and (T_H) is a constant temperature in deeper soil layers.
Parameterisation of equation terms
Next we describe the functional forms of the dependence of model parameters on the temperature and the depth. Following the previous study16, we assume that the maximal bacterial growth rates (alpha (T)) and ({overline{alpha }}(T)) are described by

$$begin{aligned} alpha (T)= & {} exp left (-frac{(T-T_0)^2}{2sigma ^2}right )alpha _{text {max}}, end{aligned}$$
(3)

$$begin{aligned} {overline{alpha }}(T)= & {} alpha (T) left [1-frac{T^n}{T_1^n + T^n}right ] = alpha _{text {max}}exp {left (-frac{(T-T_0)^2}{2sigma ^2}right )} left [1-frac{T^n}{T_1^n + T^n}right ], end{aligned}$$
(4)

where (T_0=38.2 ^circ text {C}) is the optimal temperature; (T_1=34.8 ^circ text {C}) is the temperature corresponding to the switch between the lytic and the lysogenic cycles; (alpha _{text{max}}=23 text {day}^{-1}) is the maximal possible growth, (sigma =9.1 ^circ text {C}) describes the decay of growth with temperature T16,20.
In the equation for ({overline{alpha }}(T)), we assume that at a high temperature normal cell division of (I_1) stops since there is a transition to a lytic state in bacteria. In the soil bacteria grow anaerobically or microaerophillically, and the growth rates of B. pseudomallei under such conditions are yet to be studied. For simplicity they are assumed to be the same as under aerobic conditions. Realistic values of the above parameters are listed in Table 1. Note that in the model both (alpha (T)) and ({overline{alpha }}(T)) are in fact effective growth rates of the bacterial populations, i.e. they incorporate the replication of cells and as well as their mortality.
Table 1 Parameters used in the model along with their units and ranges.
Full size table

The overall adsorption rate of the phage K is estimated as (2 times 10^{-7} text {ml}^{-1} text { day}^{-1}) from Egilmez et al.16. The adsorption constants (K_1 (T)), (K_2 (T)) and the transition rate from lysogenic to lytic cycle (lambda _1(T)) depend on temperature as follows16:

$$begin{aligned} K_1(T)&= left(1-frac{T^n}{T_1^n + T^n}right) K_S, nonumber \ K_2(T)&= frac{T^n}{T_1^n + T^n} K_S, nonumber \ lambda _1(T)&= frac{T^n}{T_1^n + T^n} {lambda _1}_{text {max}} , end{aligned}$$
(5)

where (K_S) is the maximal phage adsorption constant ((K_S=epsilon K) where (epsilon =0.3) is the adsorption efficiency) and (lambda _{1_text {max}}=23 text {day}^{-1}) is the maximal transition rate which is assumed to be equal to the maximal growth rate of the bacteria16. The lysis rate of bacteria (lambda _2=20 text {day}^{-1}) (depending on 50 min latency time13) and the burst size (b = 100) in the model are assumed to be constant16. The temperature dependence of (alpha (T)), ({overline{alpha }}(T)), (K_1 (T)), (K_2 (T)) and (lambda _1(T)) are shown in Fig. 2. The mortality rate of phages (mu) is high near the surface due to ultraviolet radiation, but the role of ultraviolet radiation becomes negligible starting from a depth of a few centimetres because sunlight cannot penetrate the soil. For the above reason, we can assume (mu =3 text {day}^{-1}) to be constant.
Figure 2

(a) Temperature dependence of the adsorption constants (K_i) ((i=1,2)) of the phage (measured in (text {ml}^{-1} text {day}^{-1})). (b) Growth rates of susceptible (alpha (T)) and lysogenic ({overline{alpha }}(T)) bacteria and the transition rate (lambda _1(T)) from the lysogenic cycle to the lytic cycle (measured in (text {day}^{-1})). The corresponding analytical expressions for the temperature dependence are given by (3)–(5).

Full size image

The carrying capacity C of the bacteria varies with the depth of the soil, according to empirical observations21,22,23. This can be explained by the fact that the humus, oxygen, nitrogen contents, or/and water content in the soil generally decrease with depth24. We use a combined approach to parameterise C(h) based on the available empirical data. We assume that in the absence of phages, the bacteria achieve numbers close to the carrying capacity at a given depth. Firstly, we parameterise the dependence of the overall bacterial load on depth in paddy soils in Southern Asia using the existing data22. Then we re-scale the obtained curve based on the available observations of B. pseudomallei at a depth (h=30 text {cm})25,26. We approximate C(h) using the following simple Gaussian-type curve

$$begin{aligned} C(h)=(C_text {surf} -C_0)exp (-B h^2)+C_0, end{aligned}$$
(6)

where (C_text {surf}) gives the maximal number of bacteria near the surface (h), B determines how fast the bacterial abundance decreases with depth, (C_0) is background bacterial density which takes into account the fact that bacteria can survive even at large depths (e.g. (h=100 text {cm})). Based on our estimates (see supplementary material SM1 for more detail), we will use the following parameter values as defaults: (C_text {surf} = 1 times 10^6) (text {cell/ml}), (B=7.5 times 10^{-4}) (1/{text{cm}}^2), (C_0=10^4) (text {cell/ml}) . One can easily see that C(h) has a maximum at the surface and monotonically decreases with depth. We assume that the carrying capacity of the environment is not influenced by seasonal variation.
The coefficient (D_h) in the equation for the temperature distribution can be estimated as follows. Generally, (D_h) is related to (rho _s), (C_{rho s}) and (k_s) which are the bulk density, specific heat and thermal conductivity in soil, respectively, i.e. (D_h=k_s/(rho _s C_{rho s)}). We use the estimates for (rho _s), (C_{rho s}) and (k_s) from27 which gives (rho _s=110.52 text {kg}/text{m}^3), (C_{rho s} = 1130) (text {J/kg K}) and (k_s = 0.0967) (text {W/m K}) and, for the diffusion coefficient (D_h=7.7 times 10^{-8}) (text {m}^2 text{s}^{-1}). The variation of (T_s)—the surface temperature—is obtained from the historical weather report for the surface16. The bottom boundary temperature (T_H) at (h=H=1 text {m}) is considered to be (22 ^circ text{C}). The initial value of the temperature distribution (T_s (0)) is assumed to be linear, but this assumption does not affect long-term temperature dynamics.
The paddy fields in which we model the bacteria-phage interactions are flooded lands, where the soil is either mud or muddy water. Many factors can affect vertical dispersal of bacteria and phages in such soil. For instance, rain water can carry bacteria and phage up or down in the soil, which can be mathematically modelled by adding an advection term; however, for simplicity we ignore such effects in this paper. We also assume the phage and bacteria vertical diffusion coefficients to be constant; however, it is rather hard to provide accurate estimates for (D_p) and (D_b). In water, the diffusion coefficient of bacteria and phages can be estimated as (3.6times 10 ^{-10} text {m}^2 text{s}^{-1}= 0.3 text {cm}^2 text{day}^{-1}) and (2.8 times 10^{-12} text {m}^2 text{s}^{-1}= 0.002 text {cm}^2 text{day}^{-1}), respectively28, but the diffusivity in soil should be smaller than this. As such, these values should be considered as upper limits for (D_P) and (D_b), with the actual coefficients being orders of magnitude smaller. We undertook simulations with different combinations of diffusion coefficients in this range, and found that the patterns of vertical distribution do not largely depend on the diffusion coefficients provided (D_P< 10^{-3} text {cm}^2 text{day}^{-1}) and (D_b < 10^{-2} text {cm}^2 text{day}^{-1}), due to the strong external forcing of the system by temperature (see “Results” section for details). In our numerical simulations, we use both explicit and implicit numerical schemes. We take a 0.1 cm spatial step size to get a proper resolution. We separately compute the heat equation to define T(t) with a smaller time resolution and then apply the temperature obtained to model bacteria-phage interactions for a larger time resolution (for example, (Delta t cong 7 times 10^{-5}) day). We compute the average densities of the species (both in terms of spatial and temporal averaging) using a numerical right Riemann sum. The accuracy of our numerical simulation was verified by reducing both time and space steps and comparing the results obtained. We use daily and seasonal variation of temperatures (for the period of 2013–2016) in the provinces of Nakhon Phanom and Sa Kaeo in Thailand to parameterise the model (http://www.worldweatheronline.com). The unit of the densities of bacteria and phages are cells/ml. The summary of model parameters as well their values are provided in Table 1. More

Elements of TA
Identifying biogeographic trends within volatile data required several innovations in the MASTIF model20, building from multivariate state-space methods in previous applications41,52. Standard modeling options, such as generalized linear models and their derivatives, do not accommodate key features of the masting processes. First, multiple data types are not independent. Maturation status is binary with detection error, CCs are non-negative integers, also with detection error, and STs require a transport model (dispersal) linking traps to trees, and identification error in seed identification. Of course, a tree observed to bear seed, now or in the past, is known to be mature now. However, failure to observe seed does not mean that an individual is immature because there are detection errors and failed crop years41,64.
Second, seed production is quasiperiodic within an individual (serial dependence), quasi-synchronous between individuals (“mast years”), [and] there is dependence on environmental variation, and massive variation within and between trees41,53,65. Autoregressive error structures (AR(p) for p lag terms) impose a rigid assumption of dependence that is not consistent with quasiperiodic variation that can drift between dominant cycles within the same individual over time43. It does not allow for individual differences in mast periodicity.
Third, climate variables that affect fecundity operate both through interannual anomalies over time and as [a] geographic variation. The masting literature deals almost exclusively with the former, but our application must identify the latter: the potentially smooth variation of climate effects across regions must be extracted from the many individual time series, each dominated by local “noise.”
Finally, model fitting is controlled by the size classes that dominate a given site and thus is insensitive to size classes that are poorly represented. Large trees are relatively rare in eastern forests, making it hard to identify potential declines in large, old individuals41,53. Conversely, the shade-intolerant species that dominate second-growth forests often lack the smaller size classes needed to estimate maturation and early stages where fecundity may be increasing rapidly.
Several of the foregoing challenges are resolved in the MASTIF model by introducing latent states for individual maturation status and tree-year seed production. The dependent data types (maturation status, CCs, STs) become conditionally independent in the hierarchical MASTIF model (e.g., ref. 66). The serial dependence is handled as a conditional hidden Markov process for maturation that combines with CCs and STs by way of stochastic (latent) conditional fecundity. Maturation status and conditional fecundity must be estimated jointly, that is, not with separate models. The latent maturation/fecundity treatment avoids imposing a specific AR(p) structure. In the MASTIF model there is a posterior covariance in maturation/fecundity across all tree-year estimates that need not adhere to any specific assumption20. The dependence across individuals and years is automatic and available from the posterior distribution.
Separating the spatial from temporal components of climate effects is possible here, not only because the entire network is analyzed together but also because predictors in the model include both climate norms for the individual sites and interannual anomalies across sites35,52. TA depends on both of these components.
Extracting the trends from volatile data further benefits from random individual effects for each tree and the combination of climate anomalies and year effects over time. A substantial literature focuses on specific combinations of climate variables that best explain year-to-year fecundity variation, including combinations of temperature, moisture, and water balance during specific seasons over current and previous years19,25,41. Results vary for each study, presumably due to the differences in sites, species, size classes, duration, data type, and modeling assumptions. For TA, the goal is to accommodate the local interannual variation to optimize identification of trends in space and time. Thus, we include a small selection of important climate anomalies (spring minimum T of the current year, summer T of the current and previous year, and moisture D of the current and previous year). The climate anomalies considered here do not include every variable combination that could be important for all size classes of every species on every site. For this reason, we combine climate anomalies with year effects. Year effects in the model are fixed effects within an ecoregion and random between ecoregions (ecoregions are shown in Fig. 2 and listed in Supplementary Data 2). They are fixed within an ecoregion because they are not interpreted as exchangeable and drawn at random from a large population of possible years. They are random between ecoregions due to the uneven distribution of sites (Supplementary Data 1)20.
To optimize inference on size effects, the sampling of coefficients in posterior simulation is implemented as a weighted regression. This means that the contribution of tree diameter to fecundity is inversely proportional to the abundance of that size class in the data. This approach has the effect of balancing the contributions of abundant and rare sizes. Identifying size effects further benefits from the introduction of opportunistic field sampling, which can target the large individuals that are typically absent from field study plots.
MASTIF data network
Data included in the analysis come from published and unpublished sources and offer one or both of the two data types, CCs and STs (Supplementary Data 1). Both data types inform tree-year fecundity; they are plotted by year in Fig. 6.
Fig. 6: Distribution of observation trees by year in the North American region of the MASTIF network.

Sites are listed by ecoregion in the Supplementary Data 2.

Full size image

CCs in the MASTIF network are obtained by one of three methods. Most common are counts with binoculars that are recorded with an estimate of the fraction of the crop that was observed. A second CC method makes use of seeds collected per ground surface area relative to the crown area. This method is used where conspecific crowns are isolated and wind dispersal is limited. The crop fraction is the ratio of ground area for traps relative to the projected crown area. Examples include HNHR67 and BCEF68.
A third CC method is based on evidence for past cone production that is preserved on trees. This has been used for Abies balsamea at western Quebec sites69, Pinus ponderosa in the Rocky Mountains70, and for Pinus edulis at SW sites27.
ST data include observations on individual trees that combine with seed counts from traps. Because individual studies can report different subcategories of seeds, and few conduct rigorous tests of viability, we had to combine them using the closest description to the concept of “viable”. For example, we do not include empty conifer seeds. A dispersion model provides estimates of seeds derived from each tree. ST and CC studies are listed in Supplementary Data 1. The likelihoods for CCs and STs are detailed in ref. 20. Individually and in combination, the two data types provide estimates, with full uncertainty, for fecundity across all tree-years.
Fitted species had multiple years of observations from multiple sites, which included 211,146 trees and 2,566,594 tree-years from 123 species. Sites are shown in Fig. 2 of the main text by ecoregion, they are named in Fig. 1 and summarized in Supplementary Data 1. For TA the fits were applied to 7,723,671 trees on inventory plots. Mean estimates for the genus were used for inventory trees belonging to species for which there were not confident fits in the MASTIF model, which amounted to 7.2% of inventory trees. Detailed site information is available at the website MASTIF.
Covariates
Covariates in the model include as main effects tree diameter, tree canopy class (shading), and the climate variables in Fig. 1 of the main text and described in Table 1. A quadratic diameter term in the MASTIF model allows for changes in diameter response with size52. Shade classes follow the USDA Forest Inventory and Analysis (FIA)/National Ecological Observation Network (NEON) scheme that ranges from a fully exposed canopy that does not interact with canopies of other trees to fully shaded in the understory. Shading provides information on competition that has proved highly significant for fecundity in previous analyses41,52.
Table 1 Predictors in the model, not all of which are important for all species.
Full size table

To distinguish between the effects of spatial variation versus interannual variability, spring T and moisture D are included in the model as site means and site anomalies35. Spring minimum T affect phenology and frost risk during flowering and early fruit initiation. Summer mean T (June–August) is included both as a linear and quadratic term. Mean summer T is linked to thermal energy availability during the growing season, with the quadratic term allowing for potential suppression due to extreme heat. Moisture D (cumulative monthly PET-P (potential evapotranspiration[-] minus precipitation) for January–August) is included as a site mean and an annual anomaly. Moisture D is important for carbon assimilation and fruit development during summer in the eastern continent and, additionally, from the preceding winter in the western continent. For species that develop over spring and summer, anomalies incorporate the current and previous year. We did not include longer lags in covariates. For species that disperse seed in spring (Ulmus spp. and some members of Acer), only the previous year was used. Temperature anomalies were included for spring, but not summer, simply to reduce the number of times that temperature variables enter the model, and these two variables tended to be correlated at many sites.
Climate covariates were derived from gridded climate products and combined with local climate monitoring where it is available. Terraclimate71 provides monthly resolution, but it is spatially coarse. For both norms and trends, we used the period from 1990 to 2019 because global temperatures have been increasing consistently since the 1980s, and this span broadly overlaps with fecundity data (Fig. 6). CHELSA72 data are downscaled to a 1 km grid, but it does not extend to 2019. Our three-component climate scaling used regression to project CHELSA forward using Terraclimate, followed by downscaling to 1 km with CHELSA, with further downscaling to local climate data. Even where local climate data exist, they often do not span the full duration of field studies, making the link to gridded climate data important. Local climate data were especially important for mountainous sites in the Appalachians, Rockies, Sierra Nevada, and Cascades.
Of the full list of variables, a subset was retained, depending on species (some have narrow geographic ranges) and deviance information criteria of the fitted model (Supplementary Data 2). Year effects in the model allow for the interannual variation that is not absorbed by anomalies20.
Model fitting and TA
As mentioned above, model fitting applied the hierarchical Bayes model of ref. 20 to the combination of time series and opportunistic observations summarized in Fig. 1. Posterior simulation was completed with Markov chain Monte Carlo based on direct sampling, Metropolis, and Hamiltonian Markov chain. Model fitting used 211,146 trees and 2,566,594 tree-years from 123 species (Supplementary Data 2). Only species with multiple observation years were included.
The climate variable referenced as C in Eq. (1) of the main text is, in fact, a vector of climate variables described in the previous section, spring minimum T, summer mean T, and moisture D (Table 1). The anomalies and year effects in the fitted model contribute to the trends not explained by biogeographic variation as γ in Eq. (1). For main effects in the model, the partial derivatives are fitted coefficients, an example being the response to spring minimum temperature (partial f/partial {T}_{mathrm{sp}}={beta }_{{T}_{mathrm{sp}}}). For predictors involved in interactions, the partial derivatives are combinations of fitted coefficients and variables. For example, the response to moisture D, which interacts with tree size, is (partial [F], f/partial {D}={beta }_{{D}} + beta_{GD}G). The response to diameter G, which is quadratic and interacts with D, is (partial f/partial G={beta }_{G}+2{beta }_{{G}^{2}}G ,+{beta }_{GD}D).
Trend decomposition applied the fitted model to every tree in forest inventories from the USDA FIA program, the Canada’s National Forest Inventory, the NEON, and our MASTIF collaboration. Each tree in these inventories has a species and diameter. For trees that lack a canopy class, regression was used to predict it from distances and tree diameters based on inventories that include both location and canopy class, including NEON, FIA, and the MASTIF network. Although inventories differ in the minimum diameter they record, few trees are reproductive at diameters below the lower diameter limits in these surveys, so the effect on fecundity estimates is negligible. For the indirect effects of climate coming through tree growth rates, the same covariates were fitted to growth as previously defined for fecundity, using the change in diameter observed over multiple inventories. A Tobit model was used to accommodate the fact that a second measurement can be smaller than an earlier measurement. The Tobit thus treats negative growth as censored at zero. TA to inventory plots used 7,717,677 trees. Because not all species in the inventory data are included in the MASTIF network, mean fecundity parameters for the genus were used for unfitted species. Species fitted in the MASTIF network accounted for >90% of trees in inventory plots (Supplementary Data 2).
From the predictive distributions for every tree in the inventory data, we evaluated predictive mean trends aggregated to species and plot in Fig. 2b. We extracted specific terms that comprise the components in Fig. 4 and aggregated them too to the plot averages.
General form for TA
Equation 1 simplifies the model to highlight direct and indirect effects. Again, climate variables and tree size represent only a subset of the predictors in the model that are collected in a design vector ({{bf{x}}}_{t}=[{x}_{1,t},ldots ,{x}_{Q,t}]^{prime}), where the q = 1, …, Q predictors include shading from local competition, individual size, and climate and habitat variables (Table 1). On the proportionate scale, Eq. (1) can be written in terms of all predictors, including main effects and interactions, as

$$frac{{mathrm{d}}f}{{mathrm{d}}t}=mathop{sum }limits_{q=1}^{Q}left(frac{partial f}{partial {x}_{q}}+sum _{q^{prime} in {I}_{q}}frac{partial f}{partial ({x}_{q}{x}_{q^{prime} })}{x}_{q^{prime} }right)frac{{mathrm{d}}{x}_{q}}{{mathrm{d}}t}+gamma$$
(2)

where Iq are variables that interact with xq. In this application, interactions include tree diameter with moisture deficit and diameter squared. Each term in the summation consists of a main effect of xq and interactions that are multiplied by the rate of change in variable xq. For the simple case of only two predictors, Eq. (2) is recognizable as Eq. (1) of the main text, where x1, x2 have been substituted for variables G and C. In our application, predictors include additional climate and shading (Table 1).
Recognizing that environmental variables affect not only fecundity but also growth rate, we extract the size effect, that is, the xq that is G, and incorporate these indirect effects (through growth) by expanding g = dG/dt in Eq. (1) of the main text as

$$g=mathop{sum }limits_{q=1}^{Q}left(frac{partial g}{partial {x}_{q}}+mathop{sum}limits _{q^{prime} in {I}_{q}}frac{partial g}{partial ({x}_{q}{x}_{q^{prime} })}{x}_{q^{prime} }right){x}_{q}+nu$$
(3)

where ν is the component of growth that is not accommodated by other terms. This expression allows us to evaluate the full effect of climate variables, including those coming indirectly through growth.
Connecting fitted coefficients in MASTIF to TA
This section connects the continuous, deterministic Eq. (1) to the MASTIF model20 with the interpretation of responses, direct effects, and full effects of Fig. 5. To summarize key elements of the fitted model20, consider a tree i at site j that grows to reproductive maturity and then produces seed depending on its size, local competitive environment, and climate. We wish to estimate the effects of its changing environment and condition on fecundity using a model that includes spatial variation in predictors that are tracked longitudinally over years t. Fecundity changes through maturation probability ρij(t), which increases as trees increase in size, and through conditional fecundity ψij(t), the annual seed production of a mature tree. Let zij(t) = 1 be the event that a randomly selected tree i is mature in year t. Then, ρij(t) is the corresponding probability that the tree is mature, E[zij(t)] = ρij(t)(ρ is not to be confused with the probability that a tree that is now immature will make the transition to the mature state in an interval dt = 1. That is a different quantity detailed in the Supplement to ref. 41). Fecundity has expected value Fij(t) = ρij(t)ψij(t). On a proportionate (log) scale,

$${f}_{ij}(t)={mathrm{log}},{F}_{ij}(t)={mathrm{log}},{rho }_{ij}(t)+{mathrm{log}},{psi }_{it}(t)$$
(4)

The corresponding rate equation is

$$frac{{mathrm{d}}f}{{mathrm{d}}t}=frac{{mathrm{d}},{mathrm{log}},rho }{{mathrm{d}}t}+frac{{mathrm{d}},{mathrm{log}},psi }{{mathrm{d}}t}$$
(5)

The discretized and stochasticized version of Eq. (1) is

$$frac{{mathrm{d}}{F}_{ij}}{{mathrm{d}}t} = , frac{{F}_{ij,t+{mathrm{d}}t}-{F}_{ij,t}}{{mathrm{d}}t}+{epsilon }_{ij,t}\ = , {{Delta }}{F}_{ij,t}+{epsilon }_{ij,t}$$
(6)

where dt = 1 and ϵij,t is the integration error. When applied to a dynamic process model, this term further absorbs process error (see above), which is critical here to allow for conditional independence where observations are serially dependent. In simplest terms, ϵ is model miss-specification that allows for dependence in data.
The MASTIF model that provides estimates for TA is detailed in ref. 20. Elements of central interest for TA are

$${F}_{ij,t} = , {z}_{ij,t}{psi }_{ij,t}\ left[{z}_{ij,t}=1right] sim , {{Bernoulli}}left({rho }_{ij,t}right)\ {rho }_{ij,t} = , {{Phi }}({{boldsymbol{mu }}}_{ij,t})\ mathrm{log},{psi }_{ij,t} = ,{{bf{x}}}_{ij,t}^{prime}{boldsymbol{beta }}+{h}_{t}left(Tright)+{epsilon }_{ij,t}$$

where μij,t = α0 + αGGij,t describes the increase in maturation probability with size, Φ(⋅) is the standard normal distribution function (a probit), ϵij,t ~ N(0, σ2), and ht(T) can include year effects, h(T) = κt, or lagged effects, (h(T)=mathop{sum }nolimits_{k = 1}^{p}{kappa }_{k}{psi }_{ij,t-k}), that contribute to γ in Eq. (1) of the main text. If year effects are used, then γ includes the trend in year effects. (The generative version of this model writes individual states at t conditional on t − 1 and is given in the Supplement to ref. 20.). If an AR(p) model is used, then γ = κ1 (provided data are not detrended). Random individual effects in the fitted model are marginalized for prediction of trees that were not fitted, meaning that σ2 is the sum of model residual and random-effects variance. Again, the length-Q design vector xij,t includes individual attributes (e.g., diameter Gij,t), local competitive environment, and climate (Table 1). There is a corresponding coefficient vector β.
Moving to a difference equation (rate of change) for conditional log fecundity,

$${{Delta }}{f}_{ij,t}={{Delta }}mathrm{log},{rho }_{ij,t}+{{Delta }}mathrm{log},{psi }_{ij,t}$$

where

$${{Delta }}mathrm{log},{psi }_{ij,t} =mathrm{log},{psi }_{ij,t+1}-mathrm{log},{psi }_{ij,t}\ ={{Delta }}{{bf{x}}}_{ij,t}^{prime}{boldsymbol{beta }}+{gamma }_{ij,t}+{nu }_{ij,t}\ {{Delta }}{{bf{x}}}_{ij,t} ={{bf{x}}}_{i,t}-{{bf{x}}}_{ij,t-1}\ {nu }_{ij,t} sim N(0,2{sigma }^{2})$$

The variance in the last line is the variance of the difference Δϵij,t.
Elements of basic theory in Eq. (1) of the main text are linked to data through the modeling framework as

$${{Delta }}{f}_{ij,t}= +{beta }_{{T}_{sp}}{{Delta }}{T}_{sp,j}\ +left({beta }_{T}+2{beta }_{{T}^{2}}{T}_{j}right){{Delta }}{T}_{j}\ +left({beta }_{D}+{beta }_{GD}{G}_{ij,t}right){{Delta }}{D}_{j}\ +left({alpha }_{G}frac{phi ({{boldsymbol{mu }}}_{ij,t})}{{{Phi }}({{boldsymbol{mu }}}_{ij,t})}+{beta }_{G}+2{beta }_{{G}^{2}}{G}_{ij,t}+{beta }_{GD}{D}_{j}right){{Delta }}{G}_{ij}\ +{gamma }_{ij,t}+{nu }_{ij,t}$$
(7)

where ϕ(⋅) is the standard normal density function that comes from the rate of progress toward maturation. Again, the anomalies do not appear in this expression for trends because trends in the anomalies and year effects enter through γ.
The first four lines in Eq. (7) are, respectively, the effects of trends in spring minimum temperatures ΔTsp,j, summer mean temperature ΔTj, moisture deficit ΔDj, and size ΔGij, where the latter comes from growth on inventory plots. The contribution of maturation to change in fecundity is the first term in the fourth line, αGϕ(μij,t)/Φ(μij,t). A map of this term in Fig. 7b shows the strong contribution to fecundity in the East due to the young (Fig. 7a) and/or small (Fig. 4b) trees there. The sum of these terms dominates the patterns in Fig. 3c.
Fig. 7: Size and maturation effects on fecundity.

a Stand age variable in FIA data and b positive effect of maturation for increasing fecundity in the eastern continent. In the West, maturation does not contribute to rising fecundity because large trees are predominantly [mature] larger.

Full size image

All terms in Eq. (7) have units of mean change in proportionate fecundity, and these are mapped in figures of the main text. We focus on proportionate fecundity because it reflects the full effect of climate as opposed to total fecundity, which would often be dominated by one or a few trees of a single species. However, from proportionate fecundity we can obtain change in fecundity as ΔFij,t = Fij,t × Δfij. Stand-level effects on fecundity change at site j can be obtained from individual change as

$${{Delta }}{F}_{j}=mathop{sum }limits_{i=1}^{{n}_{j}}{{Delta }}{F}_{ij}=mathop{sum }limits_{i=1}^{{n}_{j}}{F}_{ij}{{Delta }}{f}_{ij,t}$$

Again, maps in Fig. 5 show mean proportionate effects for all trees on an inventory plot.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Van der Laan, R. & Fricke, R. Eschmeyer’s Catalog of Fishes Family Group Names. http://www.calacademy.org/scientists/catalog-of-fishes-family-group-names (2020).
2.
Greenwood, P. H. African cichlids and evolutionary theories. In Evolution of Fish Species Flock (eds Echelle, A. A. & Kornfield, I.) 141–154 (University of Maine at Orono Press, Orono, 1984).
Google Scholar 

3.
Muschick, M., Indermaur, A. & Salzburger, W. Convergent evolution within an adaptive radiation of cichlid fishes. Curr. Biol. 22, 2362–2368 (2012).
CAS  PubMed  Article  Google Scholar 

4.
Wagner, C. E., Harmon, L. J. & Seehausen, O. Ecological opportunity and sexual selection together predict adaptive radiation. Nature 487, 366–369 (2012).
ADS  CAS  PubMed  Article  Google Scholar 

5.
Tiercelin, J.-J. & Mondeguer, A. The geology of the Tanganyika trough. In Lake Tanganyika and its Life (ed. Coulter, G. W.) 7–48 (Oxford University Press, Oxford, 1991).
Google Scholar 

6.
Irisarri, I. et al. Phylogenomics uncovers early hybridization and adaptive loci shaping the radiation of Lake Tanganyika cichlid fishes. Nat. Commun. 9, 3159 (2018).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

7.
Salzburger, W., Meyer, A., Baric, S., Verheyen, E. & Sturmbauer, C. Phylogeny of the Lake Tanganyika Cichlid species flock and its relationship to the Central and East African Haplochromine Cichlid Fish Faunas. Syst. Biol. 51, 113–135 (2002).
PubMed  Article  Google Scholar 

8.
Salzburger, W., Mack, T., Verheyen, E. & Meyer, A. Out of Tanganyika: genesis, explosive speciation, key-innovations and phylogeography of the haplochromine cichlid fishes. BMC Evol. Biol. 5, 17 (2005).
PubMed  PubMed Central  Article  Google Scholar 

9.
Koblmüller, S. et al. Age and spread of the haplochromine cichlid fishes in Africa. Mol. Phylogenet. Evol. 49, 153–169 (2008).
PubMed  Article  CAS  Google Scholar 

10.
Sturmbauer, C., Salzburger, W., Duftner, N., Schelly, R. & Koblmüller, S. Evolutionary history of the Lake Tanganyika cichlid tribe Lamprologini (Teleostei: Perciformes) derived from mitochondrial and nuclear DNA data. Mol. Phylogenet. Evol. 57, 266–284 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

11.
Sturmbauer, C., Levinton, J. S. & Christy, J. Molecular phylogeny analysis of fiddler crabs: test of the hypothesis of increasing behavioral complexity in evolution. Proc. Natl. Acad. Sci. U. S. A. 93, 10855–10857 (1996).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

12.
Joyce, D. A. et al. An extant cichlid fish radiation emerged in an extinct Pleistocene lake. Nature 435, 90–95 (2005).
ADS  CAS  PubMed  Article  Google Scholar 

13.
Katongo, C., Koblmüller, S., Duftner, N., Mumba, L. & Sturmbauer, C. Evolutionary history and biogeographic affinities of the serranochromine cichlids in Zambian rivers. Mol. Phylogenet. Evol. 45, 326–338 (2007).
CAS  PubMed  Article  Google Scholar 

14.
Sturmbauer, C., Koblmüller, S., Sefc, K. M. & Duftner, N. Phylogeographic history of the genus Tropheus, a lineage of rock-dwelling cichlid fishes endemic to Lake Tanganyika. Hydrobiologia 542, 335–366 (2005).
Article  Google Scholar 

15.
Meier, J. I. et al. Ancient hybridization fuels rapid cichlid fish adaptive radiations. Nat. Commun. 8, 14363 (2017).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

16.
Svardal, H. et al. Ancestral hybridization facilitated species diversification in the Lake Malawi Cichlid fish adaptive radiation. Mol. Biol. Evol. 37, 1100–1113 (2020).
PubMed  Article  Google Scholar 

17.
Kullander, S. O. & Roberts, T. R. Out of Tanganyika: endemic lake fishes inhabit rapids of the Lukuga River. Ichthyol. Explor. Freshw. 22, 355–376 (2011).
Google Scholar 

18.
West-Eberhard, M.-J. Developmental Plasticity and Evolution (Oxford University Press, Oxford, 2003).
Google Scholar 

19.
Rossiter, A. The Cichlid fish assemblages of Lake Tanganyika: ecology, behaviour and evolution of its species flocks. In Advances in Ecological Research (eds Begon, M. & Fitter, A. H.) 187–252 (Academic Press Ltd., London, 1995).
Google Scholar 

20.
Malinsky, M. et al. Whole-genome sequences of Malawi cichlids reveal multiple radiations interconnected by gene flow. Nat. Ecol. Evol. 2, 1940–1955 (2018).
PubMed  PubMed Central  Article  Google Scholar 

21.
Brawand, D. et al. The genomic substrate for adaptive radiation in African cichlid fish. Nature 513, 375–381 (2014).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

22.
Liem, K. F. Evolutionary strategies and morphological innovations: Cichlid Pharyngeal Jaws. Syst Biol. 22, 425–441 (1973).
Google Scholar 

23.
Carleton, K. L., Dalton, B. E., Escobar-Camacho, D. & Nandamuri, S. P. Proximate and ultimate causes of variable visual sensitivities: Insights from cichlid fish radiations. Genesis 54, 299–325 (2016).
PubMed  PubMed Central  Article  Google Scholar 

24.
Maan, M. E. & Sefc, K. M. Colour variation in cichlid fish: Developmental mechanisms, selective pressures and evolutionary consequences. Semin. Cell. Dev. Biol. 24, 516–528 (2013).
PubMed  PubMed Central  Article  Google Scholar 

25.
Salzburger, W. Understanding explosive diversification through cichlid fish genomics. Nat. Rev. Genet. 19, 705–717 (2018).
CAS  PubMed  Article  Google Scholar 

26.
Malinsky, M. Andinoacara coeruleopunctatus Genome Browser Gateway. http://em-x1.gurdon.cam.ac.uk/cgi-bin/hgGateway?hgsid=6400&clade=vertebrate&org=A.+coeruleopunctatus&db=0 (2015).

27.
Conte, M. A. et al. Chromosome-scale assemblies reveal the structural evolution of African cichlid genomes. GigaScience 8, giz030 (2019).
PubMed  PubMed Central  Article  CAS  Google Scholar 

28.
Thibaud-Nissen, F. et al. P8008 the NCBI eukaryotic genome annotation pipeline. J. Anim. Sci. 94, 184 (2016).
Article  Google Scholar 

29.
Zerbino, D. R. et al. Ensembl 2018. Nucleic Acids Res. 46, D754–D761 (2018).
CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
Conte,M.A., Gammerdinger,W.J., Bartie,K.L., Penman,D.J. & Kocher,T.D. A high quality assembly of the Nile Tilapia (Oreochromis niloticus) genome reveals the structure of two sex determination regions. bioRxiv https://doi.org/10.1101/099564 (2017).

31.
Vij, S. et al. Chromosomal-level assembly of the Asian Seabass genome using long sequence reads and multi-layered scaffolding. PLoS Genet. 12, e1005954 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

32.
Smit, A. F. A., Hubley, R. & Green, P. RepeatMasker Open-4.0. http://www.repeatmasker.org (2015).

33.
Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 29, 24–26 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

34.
Simão, F. A., Waterhouse, R. M., Ioannidis, P., Kriventseva, E. V. & Zdobnov, E. M. BUSCO: Assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics 31, 3210–3212 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

35.
Parra, G., Bradnam, K. & Korf, I. CEGMA: A pipeline to accurately annotate core genes in eukaryotic genomes. Bioinformatics 23, 1061–1067 (2007).
CAS  PubMed  Article  Google Scholar 

36.
Dohmen, E., Kremer, L. P. M., Bornberg-Bauer, E. & Kemena, C. DOGMA: Domain-based transcriptome and proteome quality assessment. Bioinformatics 32, 2577–2581 (2016).
CAS  PubMed  Article  Google Scholar 

37.
Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff. Fly 6, 80–92 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

38.
Hunt, M. et al. REAPR: a universal tool for genome assembly evaluation. Genome Biol. 14, R47 (2013).
PubMed  PubMed Central  Article  Google Scholar 

39.
Asalone, K. C. et al. Regional sequence expansion or collapse in heterozygous genome assemblies. PLoS Comput. Biol. 16, e1008104 (2020).
CAS  PubMed  PubMed Central  Article  Google Scholar 

40.
Conte, M. A. & Kocher, T. D. An improved genome reference for the African cichlid Metriaclima zebra. BMC Genomics 16, 724 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

41.
Finn, R. D. et al. The Pfam protein families database. Nucleic Acids Res. 38, D211–D222 (2010).
CAS  PubMed  Article  Google Scholar 

42.
McKenna, A. et al. The genome analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

43.
Rausch, T. et al. DELLY: Structural variant discovery by integrated paired-end and split-read analysis. Bioinformatics 28, i333–i339 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

44.
Liu, Y. et al. Comparison of multiple algorithms to reliably detect structural variants in pears. BMC Genomics 21, 61 (2020).
PubMed  PubMed Central  Article  Google Scholar 

45.
Supernat, A., Vidarsson, O. V., Steen, V. M. & Stokowy, T. Comparison of three variant callers for human whole genome sequencing. Sci. Rep. 8, 17851 (2018).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

46.
McCarthy, D. J. et al. Choice of transcripts and software has a large effect on variant annotation. Genome Med. 6, 26 (2014).
PubMed  PubMed Central  Article  Google Scholar 

47.
Gunter, H. M., Schneider, R. F., Karner, I., Sturmbauer, C. & Meyer, A. Molecular investigation of genetic assimilation during the rapid adaptive radiations of East African cichlid fishes. Mol. Ecol. 26, 6634–6653 (2017).
CAS  PubMed  Article  Google Scholar 

48.
Navon, D. et al. Hedgehog signaling is necessary and sufficient to mediate craniofacial plasticity in teleosts. Proc. Natl. Acad. Sci. U. S. A. 117, 19321–19327 (2020).
CAS  PubMed  PubMed Central  Article  Google Scholar 

49.
Boyle, E. A., Li, Y. I. & Pritchard, J. K. An expanded view of complex traits: From polygenic to omnigenic. Cell 169, 1177–1186 (2017).
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Adhikari, K. et al. A genome-wide association scan implicates DCHS2, RUNX2, GLI3, PAX1 and EDAR in human facial variation. Nat. Commun. 7, 11616 (2016).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

51.
Liu, F. et al. A genome-wide association study identifies five loci influencing facial morphology in Europeans. PLoS Genet. 8, e1002932 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

52.
Claes, P. et al. Genome-wide mapping of global-to-local genetic effects on human facial shape. Nat. Genet. 50, 414–423 (2018).
CAS  PubMed  PubMed Central  Article  Google Scholar 

53.
Lupo, G., Harris, W. A. & Lewis, K. E. Mechanisms of ventral patterning in the vertebrate nervous system. Nat. Rev. Neurosci. 7, 103–114 (2006).
CAS  PubMed  Article  Google Scholar 

54.
Dworkin, S., Boglev, Y., Owens, H. & Goldie, S. J. The role of sonic hedgehog in craniofacial patterning, morphogenesis and cranial neural crest survival. J. Dev. Biol. 4, 24 (2016).
PubMed Central  Article  PubMed  Google Scholar 

55.
Szabo-Rogers, H. L., Smithers, L. E., Yakob, W. & Liu, K. J. New directions in craniofacial morphogenesis. Dev. Biol. 341, 84–94 (2010).
CAS  PubMed  Article  Google Scholar 

56.
Zhou, H., Kim, S., Ishii, S. & Boyer, T. G. Mediator modulates Gli3-dependent Sonic hedgehog signaling. Mol. Cell Biol. 26, 8667–8682 (2006).
CAS  PubMed  PubMed Central  Article  Google Scholar 

57.
Vilhais-Neto, G. C. et al. Rere controls retinoic acid signalling and somite bilateral symmetry. Nature 463, 953–957 (2010).
ADS  CAS  PubMed  Article  Google Scholar 

58.
Clouthier, D. E., Garcia, E. & Schilling, T. F. Regulation of facial morphogenesis by endothelin signaling: Insights from mice and fish. Am. J. Med. Genet. A 152A, 2962–2973 (2010).
PubMed  PubMed Central  Article  Google Scholar 

59.
Fischer, C. et al. Complete mitochondrial DNA sequences of the Threadfin Cichlid (Petrochromis trewavasae) and the Blunthead Cichlid (Tropheus moorii) and patterns of mitochondrial genome evolution in cichlid fishes. PLoS ONE 8, e67048 (2013).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

60.
Andrews, S. FastQC A Quality Control tool for High Throughput Sequence Data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (2016).

61.
Marçais, G. & Kingsford, C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics 27, 764–770 (2011).
PubMed  PubMed Central  Article  CAS  Google Scholar 

62.
Davis, M. P. A., van Dongen, S., Abreu-Goodger, C., Bartonicek, N. & Enright, A. J. Kraken: A set of tools for quality control and analysis of high-throughput sequence data. Methods 63, 41–49 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

63.
Wingett, S. W. & Andrews, S. FastQ Screen: A tool for multi-genome mapping and quality control. F1000Res. 7, 1338 (2018).
PubMed  PubMed Central  Article  Google Scholar 

64.
Schmieder, R. & Edwards, R. Fast identification and removal of sequence contamination from genomic and metagenomic datasets. PLoS ONE 6, e17288 (2011).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

65.
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12 (2011).
Article  Google Scholar 

66.
Buffalo, V. Scythe. https://github.com/vsbuffalo/scythe (2014).

67.
CLCbio Assembly Cell. https://www.quiagenbioinformatics.com/products/clc-assembly-cell (2015).

68.
Bushnell, B., Rood, J. & Singer, E. BBMerge—Accurate paired shotgun read merging via overlap. PLoS ONE 12, e0185056 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

69.
Xu, H. et al. FastUniq: A fast de novo duplicates removal tool for paired short reads. PLoS ONE 7, e52249 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

70.
Leggett, R. M., Clavijo, B. J., Clissold, L., Clark, M. D. & Caccamo, M. NextClip: An analysis and read preparation tool for Nextera Long Mate Pair libraries. Bioinformatics 30, 566–568 (2014).
CAS  PubMed  Article  Google Scholar 

71.
Barnett, D. W., Garrison, E. K., Quinlan, A. R., Strömberg, M. P. & Marth, G. T. BamTools: a C++ API and toolkit for analyzing and managing BAM files. Bioinformatics 27, 1691–1692 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

72.
Li, H. et al. The sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).
PubMed  PubMed Central  Article  CAS  Google Scholar 

73.
Broad Institute Picard Tools. https://github.com/broadinstitute/picard (2016).

74.
Hackl, T., Hedrich, R., Schultz, J. & Förster, F. proovread: large-scale high-accuracy PacBio correction through iterative short read consensus. Bioinformatics 30, 3004–3011 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

75.
Zimin, A. V. et al. The MaSuRCA genome assembler. Bioinformatics 29, 2669–2677 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

76.
Le, H. S., Schulz, M. H., McCauley, B. M., Hinman, V. F. & Bar-Joseph, Z. Probabilistic error correction for RNA sequencing. Nucleic Acids Res. 41, e109 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

77.
Song, L. & Florea, L. Rcorrector: efficient and accurate error correction for Illumina RNA-seq reads. GigaScience 4, 48 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

78.
Liu, Y., Schröder, J. & Schmidt, B. Musket: A multistage k-mer spectrum-based error corrector for Illumina sequence data. Bioinformatics 29, 308–315 (2013).
CAS  PubMed  Article  Google Scholar 

79.
Liu,B. et al. Estimation of genomic characteristics by analyzing k-mer frequency in de novo genome projects. arXiv:1308.2012 (2013).

80.
Denisov, G. et al. Consensus generation and variant detection by Celera Assembler. Bioinformatics 24, 1035–1040 (2008).
CAS  PubMed  Article  Google Scholar 

81.
Kajitani, R. et al. Efficient de novo assembly of highly heterozygous genomes from whole-genome shotgun short reads. Genome Res. 24, 1384–1395 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

82.
Pryszcz, L. P. & Gabaldón, T. Redundans: An assembly pipeline for highly heterozygous genomes. Nucleic Acids Res. 44, e113 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

83.
Boetzer, M., Henkel, C. V., Jansen, H. J., Butler, D. & Pirovano, W. Scaffolding pre-assembled contigs using SSPACE. Bioinformatics 27, 578–579 (2011).
CAS  PubMed  Article  Google Scholar 

84.
Luo, R. et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. GigaScience 1, 18 (2012).
PubMed  PubMed Central  Article  Google Scholar 

85.
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25, 1754–1760 (2009).
CAS  PubMed  PubMed Central  Article  Google Scholar 

86.
Frith, M. C., Wan, R. & Horton, P. Incorporating sequence quality data into alignment improves DNA read mapping. Nucleic Acids Res. 38, e100 (2010).
PubMed  PubMed Central  Article  CAS  Google Scholar 

87.
English, A. C. et al. Mind the Gap: Upgrading genomes with pacific biosciences RS long-read sequencing technology. PLoS ONE 7, e47768 (2012).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

88.
Chaisson, M. J. & Tesler, G. Mapping single molecule sequencing reads using basic local alignment with successive refinement (BLASR): application and theory. BMC Bioinform. 13, 238 (2012).
CAS  Article  Google Scholar 

89.
Wences, A. H. & Schatz, M. C. Metassembler: Merging and optimizing de novo genome assemblies. Genome Biol. 16, 207 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

90.
Gurevich, A., Saveliev, V., Vyahhi, N. & Tesler, G. QUAST: Quality assessment tool for genome assemblies. Bioinformatics 29, 1072–1075 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

91.
Kosugi, S., Hirakawa, H. & Tabata, S. GMcloser: closing gaps in assemblies accurately with a likelihood-based selection of contig or long-read alignments. Bioinformatics 31, 3733–3741 (2015).
CAS  PubMed  Google Scholar 

92.
Kurtz, S. et al. Versatile and open software for comparing large genomes. Genome Biol. 5, R12 (2004).
PubMed  PubMed Central  Article  Google Scholar 

93.
Camacho, C. et al. BLAST+: Architecture and applications. BMC Bioinformatics 10, 421 (2009).
PubMed  PubMed Central  Article  CAS  Google Scholar 

94.
Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Meth. 9, 357–359 (2012).
CAS  Article  Google Scholar 

95.
Paulino, D. et al. Sealer: A scalable gap-closing application for finishing draft genomes. BMC Bioinform. 16, 230 (2015).
Article  Google Scholar 

96.
Simpson, J. T. et al. ABySS: A parallel assembler for short read sequence data. Genome Res. 19, 1117–1123 (2009).
CAS  PubMed  PubMed Central  Article  Google Scholar 

97.
Ponstingl, H. & Ning, Z. SMALT. https://www.sanger.ac.uk/science/tools/smalt-0 (2018).

98.
Birney, E., Clamp, M. & Durbin, R. GeneWise and genomewise. Genome Res. 14, 988–995 (2004).
CAS  PubMed  PubMed Central  Article  Google Scholar 

99.
Finn, R. D., Clements, J. & Eddy, S. R. HMMER web server: Interactive sequence similarity searching. Nucleic Acids Res. 39, W29–W37 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

100.
Stanke, M. & Morgenstern, B. Augustus: A web server for gene prediction in eukaryotes that allows user-defined constraints. Nucleic Acids Res. 33, W465–W467 (2005).
CAS  PubMed  PubMed Central  Article  Google Scholar 

101.
Grabherr, M. G. et al. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat. Biotechnol. 29, 644–652 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

102.
Haas, B. J. et al. De novo transcript sequence reconstruction from RNA-Seq: reference generation and analysis with Trinity. Nat. Protoc. 8, 1494–1512 (2013).
CAS  Article  Google Scholar 

103.
Haas, B. J. et al. Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies. Nucleic Acids Res. 31, 5654–5666 (2003).
CAS  PubMed  PubMed Central  Article  Google Scholar 

104.
Wu, T. D. & Watanabe, C. K. GMAP: A genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics 21, 1859–1875 (2005).
CAS  PubMed  Article  Google Scholar 

105.
Kent, W. J. BLAT—The BLAST-like alignment tool. Genome Res. 12, 656–664 (2002).
CAS  PubMed  PubMed Central  Article  Google Scholar 

106.
Oracle Inc. MySQL. https://www.mysql.com (2016).

107.
Cantarel, B. L. et al. MAKER: An easy-to-use annotation pipeline designed for emerging model organism genomes. Genome Res. 18, 188–196 (2008).
CAS  PubMed  PubMed Central  Article  Google Scholar 

108.
Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
CAS  Article  Google Scholar 

109.
Lomsadze, A., Ter-Hovhannisyan, V., Chernoff, Y. O. & Borodovsky, M. Gene identification in novel eukaryotic genomes by self-training algorithm. Nucleic Acids Res. 33, 6494–6506 (2005).
CAS  PubMed  PubMed Central  Article  Google Scholar 

110.
Korf, I. Gene finding in novel genomes. BMC Bioinform. 5, 59 (2004).
Article  Google Scholar 

111.
Schattner, P., Brooks, A. N. & Lowe, T. M. The tRNAscan-SE, snoscan and snoGPS web servers for the detection of tRNAs and snoRNAs. Nucleic Acids Res. 33, W686–W689 (2005).
CAS  PubMed  PubMed Central  Article  Google Scholar 

112.
Palmer, J. M. Funannotate: a fungal genome annotation and comparative genomics pipeline. https://github.com/nextgenusfs/funannotate (2016).

113.
Hoff, K. J., Lange, S., Lomsadze, A., Borodovsky, M. & Stanke, M. BRAKER1: Unsupervised RNA-Seq-based genome annotation with GeneMark-ET and AUGUSTUS. Bioinformatics 32, 767–769 (2016).
CAS  PubMed  Article  Google Scholar 

114.
Lomsadze, A., Burns, P. D. & Borodovsky, M. Integration of mapped RNA-Seq reads into automatic training of eukaryotic gene finding algorithm. Nucleic Acids Res. 42, e119 (2014).
PubMed  PubMed Central  Article  CAS  Google Scholar 

115.
Pertea, M. et al. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33, 290–295 (2015).
CAS  PubMed  PubMed Central  Article  Google Scholar 

116.
Trapnell, C. et al. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

117.
Haas, B. J. et al. Automated eukaryotic gene structure annotation using EVidenceModeler and the program to assemble spliced alignments. Genome Biol. 9, R7 (2008).
PubMed  PubMed Central  Article  CAS  Google Scholar 

118.
Nawrocki, E. P., Kolbe, D. L. & Eddy, S. R. Infernal 1.0: inference of RNA alignments. Bioinformatics 25, 1335–1337 (2009).
CAS  PubMed  PubMed Central  Article  Google Scholar 

119.
Griffiths-Jones, S., Bateman, A., Marshall, M., Khanna, A. & Eddy, S. R. Rfam: an RNA family database. Nucleic Acids Res. 31, 439–441 (2003).
CAS  PubMed  PubMed Central  Article  Google Scholar 

120.
Wucher,V. et al. FEELnc: A tool for Long non-coding RNAs annotation and its application to the dog transcriptome. bioRxiv https://doi.org/10.1101/064436 (2016).

121.
Thiel, T., Michalek, W., Varshney, R. K. & Graner, A. Exploiting EST databases for the development and characterization of gene-derived SSR-markers in barley (Hordeum vulgare L.). Theor. Appl. Genet. 106, 411–422 (2003).
CAS  PubMed  Article  Google Scholar 

122.
Rice, P., Longden, I. & Bleasby, A. EMBOSS: The European molecular biology open software suite. Trends. Genet. 16, 276–277 (2000).
CAS  PubMed  Article  Google Scholar 

123.
Jurka, J. W. RepBase. https://www.girinst.org/server/RepBase (2016).

124.
Smit, A. F. A. & Hubley, R. RepeatModeler Open-1.0. http://www.repeatmasker.org (2014).

125.
Price, A. L., Jones, N. C. & Pevzner, P. A. D. novo identification of repeat families in large genomes. Bioinformatics 21, i351–i358 (2005).
CAS  PubMed  Article  Google Scholar 

126.
Benson, G. Tandem repeats finder: A program to analyze DNA sequences. Nucleic Acids Res. 27, 573–580 (1999).
CAS  PubMed  PubMed Central  Article  Google Scholar 

127.
Jones, P. et al. InterProScan 5: genome-scale protein function classification. Bioinformatics 30, 1236–1240 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

128.
Huerta-Cepas, J. et al. eggNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences. Nucleic Acids Res. 44, D286–D293 (2016).
CAS  PubMed  Article  Google Scholar 

129.
Rawlings, N. D., Barrett, A. J. & Finn, R. Twenty years of the MEROPS database of proteolytic enzymes, their substrates and inhibitors. Nucleic Acids Res. 44, D343–D350 (2016).
CAS  PubMed  Article  Google Scholar 

130.
Yin, Y. et al. dbCAN: A web resource for automated carbohydrate-active enzyme annotation. Nucleic Acids Res. 40, W445–W451 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar 

131.
Petersen, T. N., Brunak, S., von Heijne, G. & Nielsen, H. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nat. Methods 8, 785–786 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

132.
Okonechnikov, K., Conesa, A. & García-Alcalde, F. Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics 32, 292–294 (2016).
CAS  PubMed  Google Scholar 

133.
Sterne-Weiler, T., Weatheritt, R. J., Best, A. J., Ha, K. C. H. & Blencowe, B. J. Efficient and accurate quantitative profiling of alternative splicing patterns of any complexity on a laptop. Mol. Cell 72, 187–200 (2018).
CAS  PubMed  Article  Google Scholar 

134.
Alexa, A., Rahnenführer, J. & Lengauer, T. Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics 22, 1600–1607 (2006).
CAS  PubMed  PubMed Central  Article  Google Scholar 

135.
Li, Y., Xiang, J. & Duan, C. Insulin-like growth factor-binding protein-3 plays an important role in regulating pharyngeal skeleton and inner ear formation and differentiation. J. Biol. Chem. 280, 3613–3620 (2005).
CAS  PubMed  Article  Google Scholar 

136.
Lin, J. M. et al. Actions of fibroblast growth factor-8 in bone cells in vitro. Am. J. Physiol. Endocrinol. Metab. 297, E142–E150 (2009).
CAS  PubMed  Article  Google Scholar 

137.
Nichols, J. T., Pan, L., Moens, C. B. & Kimmel, C. B. barx1 represses joints and promotes cartilage in the craniofacial skeleton. Development 140, 2765–2775 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

138.
Bush, J. O., Lan, Y. & Jiang, R. The cleft lip and palate defects in Dancer mutant mice result from gain of function of the Tbx10 gene. Proc. Natl. Acad. Sci. U. S. A. 101, 7022–7027 (2004).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

139.
Vieira, A. R. et al. Medical sequencing of candidate genes for nonsyndromic cleft lip and palate. PLoS Genet. 1, e64 (2005).
PubMed  PubMed Central  Article  CAS  Google Scholar 

140.
Papaioannou, V. E. The T-box gene family: Emerging roles in development, stem cells and cancer. Development 141, 3819–3833 (2014).
CAS  PubMed  PubMed Central  Article  Google Scholar 

141.
Kang, Y. J., Stevenson, A. K., Yau, P. M. & Kollmar, R. Sparc protein is required for normal growth of zebrafish otoliths. J. Assoc. Res. Otolaryngol. 9, 436–451 (2008).
PubMed  PubMed Central  Article  Google Scholar 

142.
Rosset, E. M. & Bradshaw, A. D. SPARC/osteonectin in mineralized tissue. Matrix Biol. 52–54, 78–87 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

143.
Zarelli, V. E. & Dawid, I. B. Inhibition of neural crest formation by Kctd15 involves regulation of transcription factor AP-2. Proc. Natl. Acad. Sci. U. S. A. 110, 2870–2875 (2013).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

144.
Zhang, Z., Huynh, T. & Baldini, A. Mesodermal expression of Tbx1 is necessary and sufficient for pharyngeal arch and cardiac outflow tract development. Development 133, 3587–3595 (2006).
CAS  PubMed  PubMed Central  Article  Google Scholar 

145.
Yutzey, K. E. DiGeorge syndrome, Tbx1, and retinoic acid signaling come full circle. Circ. Res. 106, 630–632 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

146.
Ghassibe-Sabbagh, M. et al. FAF1, a gene that is disrupted in cleft palate and has conserved function in Zebrafish. Am. J. Hum. Genet. 88, 150–161 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

147.
Wilm, T. P. & Solnica-Krezel, L. Essential roles of a zebrafish prdm1/blimp1 homolog in embryo patterning and organogenesis. Development 132, 393–404 (2005).
CAS  PubMed  Article  Google Scholar 

148.
Wang, L., Rajan, H., Pitman, J. L., McKeown, M. & Tsai, C. C. Histone deacetylase-associating Atrophin proteins are nuclear receptor corepressors. Genes Dev. 20, 525–530 (2006).
PubMed  PubMed Central  Article  CAS  Google Scholar 

149.
Plaster, N., Sonntag, C., Schilling, T. F. & Hammerschmidt, M. REREa/Atrophin-2 interacts with histone deacetylase and Fgf8 signaling to regulate multiple processes of zebrafish development. Dev. Dyn. 236, 1891–1904 (2007).
CAS  PubMed  Article  Google Scholar 

150.
Jordan, V. K. et al. Genotype–phenotype correlations in individuals with pathogenic RERE variants. Hum. Mutat. 39, 666–675 (2018).
CAS  PubMed  PubMed Central  Article  Google Scholar 

151.
Diepeveen, E. T., Kim, F. D. & Salzburger, W. Sequence analyses of the distal-less homeobox gene family in East African cichlid fishes reveal signatures of positive selection. BMC Evol. Biol. 13, 153 (2013).
PubMed  PubMed Central  Article  Google Scholar 

152.
Stock, D. W. et al. The evolution of the vertebrate Dlx gene family. Proc. Natl. Acad. Sci. USA 93, 10858–10863 (1996).
ADS  CAS  PubMed  Article  Google Scholar 

153.
Mark, M., Ghyselinck, N. B. & Chambon, P. Function of retinoic acid receptors during embryonic development. Nucl. Recept. Signal. 7, e002 (2009).
PubMed  PubMed Central  Article  CAS  Google Scholar 

154.
Linville, A., Radtke, K., Waxman, J. S., Yelon, D. & Schilling, T. F. Combinatorial roles for zebrafish retinoic acid receptors in the hindbrain, limbs and pharyngeal arches. Dev. Biol. 325, 60–70 (2009).
CAS  PubMed  Article  Google Scholar 

155.
Swartz, M. E., Sheehan-Rooney, K., Dixon, M. J. & Eberhart, J. K. Examination of a palatogenic gene program in Zebrafish. Dev. Dyn. 240, 2204–2220 (2011).
CAS  PubMed  PubMed Central  Article  Google Scholar 

156.
Iwata, J. et al. Transforming growth factor-beta regulates basal transcriptional regulatory machinery to control cell proliferation and differentiation in cranial neural crest-derived osteoprogenitor cells. J. Biol. Chem. 285, 4975–4982 (2010).
CAS  PubMed  Article  Google Scholar 

157.
Prochazkova, M., Prochazka, J., Marangoni, P. & Klein, O. D. Bones, Glands, Ears and More: The Multiple Roles of FGF10 in Craniofacial Development. Front Genet. 9, 542 (2018).
CAS  PubMed  PubMed Central  Article  Google Scholar 

158.
Du, J. et al. Different expression patterns of Gli1-3 in mouse embryonic maxillofacial development. Acta Histochem. 114, 620–625 (2012).
CAS  PubMed  Article  Google Scholar  More

Differentiation within A. pinguis cryptic species A
Genetic studies using combined DNA sequences from five chloroplasts (rbcL-a, matK, rpoC1, trnL-F, trnH-pabA) and 1 (ITS) nuclear genomes (4598 bp) showed some differentiation within the A. pinguis cryptic species A into three distinct groups (lineages) A1, A2 and A3 (Fig. 3). Most of investigated plants belonging to the lineage A1 originated from the Pieniny Mts. (PNN), however two samples A1 were collected at the Beskidy Mts. (BS) and one at the Tatry Mts. (T). All plants identified to A2 and A3 came from the Beskidy Mts. and the Tatry Mts., respectively. Maximum parsimony analyses of combined plastid loci and the nuclear ITS locus produced trees showing that the lineage A3 is genetically the most distinct, while A1 and A2 reveal more similarity. In the K2P mode, the percentage of variation in the sequences between lineages A1 and A2 equals to 0.20%, while for A1 and A3 it raised five times, i.e. 1.0%. The same occurred for the lineages A2 and A3, 1.0%.
Figure 3

Phylogram resulting from maximum likelihood (ML) analysis based on combined data of all sequences and showing genetic similarity and differentiation between lineages A1, A2 and A3 of A. pinguis cryptic species A. Bootstrap values are given at branches. A. maxima was used as an outgroup for tree rooting.

Full size image

Total and water soluble (active) alkaline elements
The total content of alkaline elements (Table 2) shows that the content of calcium (Ca) prevails over magnesium (Mg), potassium (K) and sodium (Na) at any investigated site, i.e. Pieniny Mts. (PNN), Beskidy Mts. (BS) and Tatry Mts. (T). Soil samples collected under the lineage A1 covered the whole three geographical distributions, where the site BS exhibited the highest Ca concentrations (31 459.6 mg kg−1) followed by PNN (16 178.8 mg kg−1) and finally T with 7 398.7 mg kg−1. Interestingly, the lineage A2 occurred only at the site BS characterised by high Ca content (27 303.4 mg kg−1), whereas A3 at the Tatry Mts. (T) where the level 22 227.6 mg kg−1 was recorded. These data imply that the lineage A1 may have developed site-specific adaptation mechanisms to various concentrations of calcium. In the case of A2 and A3, the observed Ca concentrations amounted to 27,303.4 and 22,227.6 mg kg−1, respectively and should be described as high.
Table 2 Total content of alkaline elements (Ca, Mg, K, Na) in the growth media of A. pinguis cryptic species A genetic lineage A1, A2, A3 at Pieniny, Beskidy and Tatry Mts.
Full size table

Variations in magnesium (Mg) concentrations for the lineage A1 followed another pattern differing from that observed in the case of calcium. Its contents varied accordingly: T  > BS  > PNN, with the highest levels recorded for A3 and A1 at the Tatry MTs. (T), respectively. It should be mentioned that both Ca and Mg are in most cases responsible (Ca much more) for geochemical reactions controlling the pH of the growth media. The role of potassium (K) as well as sodium (Na) is generally less pronounced in these reactions, but also their contents, which were very low appeared as the proof.
The evaluation of site-specific occurrence of the lineages A1, A2 and A3 should not be performed on the basis of total content solely of alkaline elements, since this fraction is mostly informative on the current status of Ca, Mg, K and Na. Therefore, we have tested the soil samples for recovering the concentrations expressed as active fractions (Table 3) potentially involved in the growth process of these lineages. The levels (percentage share into the total content) of active Ca are significantly low and varied as follows: PNN (3.27%)  > T (0.89%)  > BS (0.73%) for A1, but raised to 1.34% (BS) in the case of A2. The lineage A3 has recorded a concentration of 0.71%, slightly comparable to A1, but at the same site (T).
Table 3 Content of active forms of alkaline elements (Ca, Mg, K, Na) in the growth media of A. pinguis cryptic species A genetic lineage A1, A2, A3 at Pieniny, Beskidy and Tatry Mts.
Full size table

Should these Ca concentrations reflect any trend in site-specific behavior of Aneura pinguis cryptic species A. three lineages? Preliminary observations may be indicative of the calciphilous character of A1, specifically for the PNN site, followed by A2 in the case of BS. Lineages identified at the relatively lower share of active Ca, that is below 1.00% may fall into the acidophilous range. The percentage share of active Mg into its total concentrations followed similar distribution patterns like active Ca, with A1 recording 2.53% at the PNN site. Magnesium and calcium are divalent elements, which significantly control the alkalinity of soil environment.
In the case of the current study, the occurrence of this lineage (i.e. A1) at this site is not a random process. By applying the same criteria like for active Ca, it appeared that A1, A2 and A3 at the Beskidy as well as Tatry Mts. met the rule of active Mg  T (0.87).
Lineage A2 index: BS = 4.34.
Lineage A3 index: T = 1.31.
The respective pH values changed quite accordingly to the indices as shown below:
Lineage A1 site pH: BS (8.05)  > PNN (7.50)  > T (6.45).
Lineage A2 site pH: BS = 7.85.
Lineage A3 site pH: T = 7.08.
These ranges imply that genetic lineages A1 and A2 are by essence both calciphilous biotypes and may occur on sites rich in Ca, mostly alkaline as confirmed by the PNN and BS sites. On the other hand, some biotypes of the lineage A1 may be easily adapting also to low Ca concentrations, indicative of acidophilous features, as in the case of A3. Both (A1 and A3) occur at the Tatry MTs.
A detailed distribution of indices as well as respective pH is illustrated by the Figs. 4, 5 and 6, specifically for the genetic lineages A1, A2 and A3, respectively. The mean index values for the PNN site is 3.24 which discriminates the data into two groups: 60%  3.24. In the case of BS, the mean value amounted to 2.70, but for only two sampling sites. Therefore, the mean values of the singular site specific index shows a clear pattern, which strengthens the preferential adaptation of A1 in prevalence to alkalinity as follows: PPN (3.24)  > BS (2.70)  > T (0.87).
Figure 4

Singular site specific index of active forms of alkaline elements (Ca, Mg, K, Na) and pH in the growth media of A. pinguis cryptic species A genetic lineage A1 at Pieniny, Beskidy and Tatry Mts.

Full size image

Figure 5

Singular site specific index of active forms of alkaline elements (Ca, Mg, K, Na) and pH in the growth media of A. pinguis cryptic species A genetic lineage A2 at Beskidy Mts.

Full size image

Figure 6

Singular site specific index of active forms of alkaline elements (Ca, Mg, K, Na) and pH in the growth media of A. pinguis cryptic species A genetic lineage A3 at Tatry Mts.

Full size image

The genetic lineage A2 outlines a great variability in terms of the site specific index, which was slightly high (4.62) only for the sampling site BS 3–28. It should be mentioned that the mean value at this site raised up to 4.34, hence being 56% lower than the highest and next 49% higher than the lowest index. Curiously, the respective pH values did not vary significantly (7.83–7.89), which implies that A2 is decidedly calciphilous.
Indices reported in the Fig. 6 fluctuated widely from 0.66 to 2.38 with a mean of 1.31. Only two values were higher but the remaining, i.e. about 67% placed below. Such high share reveals that the genetic lineage A3 is basically acidophilus. This is decidedly outlined by significantly low values of indices as a consequence of low concentrations of active Ca. More

1.
Hiederer R, Köchy M. Global soil organic carbon estimates and the harmonized world soil database. EUR. 2011;79:25225.
Google Scholar 
2.
Scharlemann JPW, Tanner EVJ, Hiederer R, Kapos V. Global soil carbon: understanding and managing the largest terrestrial carbon pool. Carbon Manag. 2014;5:81–91.
CAS  Article  Google Scholar 

3.
Wieder WR, Bonan GB, Allison SD. Global soil carbon projections are improved by modelling microbial processes. Nat Clim Chang. 2013;3:909–12.
CAS  Article  Google Scholar 

4.
Schimel JP, Weintraub MN. The implications of exoenzyme activity on microbial carbon and nitrogen limitation in soil: a theoretical model. Soil Biol Biochem. 2003;35:549–63.
CAS  Article  Google Scholar 

5.
Huang Y, Guenet B, Ciais P, Janssens IA, Soong JL, Wang Y, et al. ORCHIMIC (v1.0), a microbe-mediated model for soil organic matter decomposition. Geosci Model Dev. 2018;11:2111–38.
CAS  Article  Google Scholar 

6.
Georgiou K, Abramoff RZ, Harte J, Riley WJ, Torn MS. Microbial community-level regulation explains soil carbon responses to long-term litter manipulations. Nat Commun. 2017;8:1223.
PubMed  PubMed Central  Article  CAS  Google Scholar 

7.
Kelleher BP, Simpson AJ. Humic substances in soils: are they really chemically distinct? Environ Sci Technol. 2006;40:4605–11.
CAS  PubMed  Article  Google Scholar 

8.
Wang C, Wang X, Pei G, Xia Z, Peng B, Sun L, et al. Stabilization of microbial residues in soil organic matter after two years of decomposition. Soil Biol Biochem. 2020;141:107687.
CAS  Article  Google Scholar 

9.
Cotrufo MF, Wallenstein M, Boot C, Denef K, Paul E. The microbial efficiency-matrix stabilization (MEMS) framework integrates plant litter decomposition with soil organic matter stabilization: do labile plant inputs form stable soil organic matter? Glob Change Biol. 2013;19:988–95.
Article  Google Scholar 

10.
Zhu X, Jackson RD, DeLucia EH, Tiedje JM, Liang C. The soil microbial carbon pump: from conceptual insights to empirical assessments. Glob Change Biol. 2020;26:6032–9.
Article  Google Scholar 

11.
Miltner A, Bombach P, Schmidt-Brücken B, Kästner M. SOM genesis: microbial biomass as a significant source. Biogeochemistry. 2012;111:41–55.
CAS  Article  Google Scholar 

12.
Torn MS, Trumbore SE, Chadwick OA, Vitousek PM, Hendricks DM. Mineral control of soil organic carbon storage and turnover. Nature. 1997;389:170–3.
CAS  Article  Google Scholar 

13.
Dwivedi D, Riley WJ, Torn MS, Spycher N, Maggi F, Tang JY. Mineral properties, microbes, transport, and plant-input profiles control vertical distribution and age of soil carbon stocks. Soil Biol Biochem. 2017;107:244–59.
CAS  Article  Google Scholar 

14.
Mikutta R, Kleber M, Torn MS, Jahn R. Stabilization of soil organic matter: association with minerals or chemical recalcitrance? Biogeochemistry. 2006;77:25–56.
CAS  Article  Google Scholar 

15.
Liang C, Balser TC. Microbial production of recalcitrant organic matter in global soils: Implications for productivity and climate policy. Nat Rev Microbiol. 2011;9:75–75.
CAS  PubMed  Article  Google Scholar 

16.
Khan KS, Mack R, Castillo X, Kaiser M, Joergensen RG. Microbial biomass, fungal and bacterial residues, and their relationships to the soil organic matter C/N/P/S ratios. Geoderma. 2016;271:115–23.
CAS  Article  Google Scholar 

17.
Liang C, Amelung W, Lehmann J, Kästner M. Quantitative assessment of microbial necromass contribution to soil organic matter. Glob Chang Biol. 2019;25:3578–90.
PubMed  Article  Google Scholar 

18.
Kögel-Knabner I. The macromolecular organic composition of plant and microbial residues as inputs to soil organic matter: fourteen years on. Soil Biol Biochem. 2017;105:A3–8.
Article  CAS  Google Scholar 

19.
Todd-Brown KEO, Randerson JT, Post WM, Hoffman FM, Tarnocai C, Schuur EAG, et al. Causes of variation in soil carbon simulations from CMIP5 Earth System Models and comparison with observations. Biogeosciences. 2013;10:1717–36.
Article  Google Scholar 

20.
Parton WJ, Schimel DS, Cole CV, Ojima DS. Analysis of factors controlling soil organic matter levels in great plains grasslands. Soil Sci Soc Am J. 1987;51:1173–9.
CAS  Article  Google Scholar 

21.
Wang G, Post WM, Mayes MA. Development of microbial‐enzyme‐mediated decomposition model parameters through steady‐state and dynamic analyses. Ecol Appl. 2013;23:255–72.
PubMed  Article  Google Scholar 

22.
Wang G, Mayes MA, Gu L, Schadt CW. Representation of dormant and active microbial dynamics for ecosystem modeling. PLoS ONE. 2014;9:e89252.
PubMed  PubMed Central  Article  CAS  Google Scholar 

23.
Wang G, Jagadamma S, Mayes MA, Schadt CW, Steinweg JM, Gu L, et al. Microbial dormancy improves development and experimental validation of ecosystem model. ISME J. 2015;9:226–37.
CAS  PubMed  Article  Google Scholar 

24.
German D, Marcelo K, Stone M, Allison S. The Michaelis–Menten kinetics of soil extracellular enzymes in response to temperature: a cross-latitudinal study. Glob Change Biol. 2012;18:1468–79.
Article  Google Scholar 

25.
Allison SD, Wallenstein MD, Bradford MA. Soil-carbon response to warming dependent on microbial physiology. Nat Geosci. 2010;3:336–40.
CAS  Article  Google Scholar 

26.
Li J, Wang G, Allison SD, Mayes MA, Luo Y. Soil carbon sensitivity to temperature and carbon use efficiency compared across microbial-ecosystem models of varying complexity. Biogeochemistry. 2014;119:67–84.
Article  Google Scholar 

27.
Wieder WR, Grandy AS, Kallenbach CM, Bonan GB. Integrating microbial physiology and physio-chemical principles in soils with the MIcrobial-MIneral Carbon Stabilization (MIMICS) model. Biogeosciences. 2014;11:3899–917.
Article  CAS  Google Scholar 

28.
Tang J, Riley WJ. Weaker soil carbon–climate feedbacks resulting from microbial and abiotic interactions. Nat Clim Chang. 2015;5:56–60.
CAS  Article  Google Scholar 

29.
Sulman BN, Moore JA, Abramoff R, Averill C, Kivlin S, Georgiou K, et al. Multiple models and experiments underscore large uncertainty in soil carbon dynamics. Biogeochemistry. 2018;141:109–23.
CAS  Article  Google Scholar 

30.
Sulman BN, Phillips RP, Oishi AC, Shevliakova E, Pacala SW. Microbe-driven turnover offsets mineral-mediated storage of soil carbon under elevated CO2. Nat Clim Change. 2014;4:1099–102.
CAS  Article  Google Scholar 

31.
Lawrence C, Neff J, Schimel J. Does adding microbial mechanisms of decomposition improve soil organic matter models? A comparison of four models using data from a pulsed rewetting experiment. Soil Biol Biochem. 2009;41:1923–34.
CAS  Article  Google Scholar 

32.
Wang X, Wang C, Cotrufo MF, Sun L, Jiang P, Liu Z, et al. Elevated temperature increases the accumulation of microbial necromass nitrogen in soil via increasing microbial turnover. Glob Change Biol. 2020;26:5277–89.
Article  Google Scholar 

33.
Throckmorton HM, Bird JA, Dane L, Firestone MK, Horwath WR. The source of microbial C has little impact on soil organic matter stabilisation in forest ecosystems. Ecol Lett. 2012;15:1257–65.
PubMed  Article  Google Scholar 

34.
Kindler R, Miltner A, Richnow H-H, Kästner M. Fate of gram-negative bacterial biomass in soil—mineralization and contribution to SOM. Soil Biol Biochem. 2006;38:2860–70.
CAS  Article  Google Scholar 

35.
Schweigert M, Herrmann S, Miltner A, Fester T, Kästner M. Fate of ectomycorrhizal fungal biomass in a soil bioreactor system and its contribution to soil organic matter formation. Soil Biol Biochem. 2015;88:120–7.
CAS  Article  Google Scholar 

36.
Derrien D, Amelung W. Computing the mean residence time of soil carbon fractions using stable isotopes: impacts of the model framework. Eur J Soil Sci. 2011;62:237–52.
Article  Google Scholar 

37.
Dormand JR, Prince PJ. A family of embedded Runge-Kutta formulae. J Comput Appl Math. 1980;6:19–26.
Article  Google Scholar 

38.
Shampine LF, Reichelt MW. The MATLAB ODE suite. Siam J Sci Comput. 1997;18:1–22.
Article  Google Scholar 

39.
Coleman TF, Li Y. On the convergence of reflective newton methods for large-scale nonlinear minimization subject to bounds. Math Program. 1994;67:189–224.
Article  Google Scholar 

40.
Coleman TF, Li Y. An interior trust region approach for nonlinear minimization subject to bounds. SIAM J Optim. 1996;6:418–45.
Article  Google Scholar 

41.
Moré JJ. The Levenberg–Marquardt algorithm: implementation and theory. In: Watson GA (ed). Numerical Analysis. Springer: Berlin, Heidelberg, 1978, p. 105–16.

42.
Leave-one-out cross-validation. In: Sammut C, Webb GI, editors. Encyclopedia of machine learning. Boston, MA: Springer USA; 2010. p. 600–1.

43.
Wang C, Qu L, Yang L, Liu D, Morrissey E, Miao R, et al. Large-scale importance of microbial carbon use efficiency and necromass to soil organic carbon. Glob Chang Biol. 2021.

44.
Farrell M, Prendergast-Miller M, Jones DL, Hill PW, Condron LM. Soil microbial organic nitrogen uptake is regulated by carbon availability. Soil Biol Biochem. 2014;77:261–7.
CAS  Article  Google Scholar 

45.
Hagerty SB, Allison SD, Schimel JP. Evaluating soil microbial carbon use efficiency explicitly as a function of cellular processes: implications for measurements and models. Biogeochemistry. 2018;140:269–83.
CAS  Article  Google Scholar 

46.
Qiao Y, Wang J, Liang G, Du Z, Zhou J, Zhu C, et al. Global variation of soil microbial carbon-use efficiency in relation to growth temperature and substrate supply. Sci Rep. 2019;9:5621.
PubMed  PubMed Central  Article  CAS  Google Scholar 

47.
Krinner G, Viovy N, de Noblet-Ducoudré N, Ogée J, Polcher J, Friedlingstein P, et al. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. Glob Biogeochem Cycles. 2005;19:GB1015.
Article  CAS  Google Scholar 

48.
Wang G, Post WM, Mayes MA, Frerichs JT, Sindhu J. Parameter estimation for models of ligninolytic and cellulolytic enzyme kinetics. Soil Biol Biochem. 2012;48:28–38.
Article  CAS  Google Scholar 

49.
Davidson EA, Janssens IA. Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. Nature. 2006;440:165–73.
CAS  PubMed  Article  Google Scholar 

50.
Fick SE, Hijmans RJ. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int J Climatol. 2017;37:4302–15.
Article  Google Scholar 

51.
Guevara M, Taufer M, Vargas R. Gap-free global annual soil moisture: 15 km grids for 1991–2018. Earth Syst Sci Data. 2020;2020:1–65.
Google Scholar 

52.
Kalnay E, Kanamitsu M, Kistler R, Collins W, Deaven D, Gandin L, et al. The NCEP/NCAR 40-year reanalysis project. Bull Am Meteorol Soc. 1996;77:437–72.
Article  Google Scholar 

53.
Batjes NH. Harmonized soil property values for broad-scale modelling (WISE30sec) with estimates of global soil carbon stocks. Geoderma. 2016;269:61–8.
CAS  Article  Google Scholar 

54.
Hengl T, Mendes de Jesus J, Heuvelink GBM, Ruiperez Gonzalez M, Kilibarda M, Blagotić A, et al. SoilGrids250m: global gridded soil information based on machine learning. PLoS ONE. 2017;12:e0169748.
PubMed  PubMed Central  Article  CAS  Google Scholar 

55.
Olson DM, Dinerstein E. The Global 200: a representation approach to conserving the earth’s most biologically valuable ecoregions. Conserv Biol. 1998;12:502–15.
Article  Google Scholar 

56.
Kögel-Knabner I. The macromolecular organic composition of plant and microbial residues as inputs to soil organic matter. Soil Biol Biochem. 2002;34:139–62.
Article  Google Scholar 

57.
Fernandez CW, Koide RT. Initial melanin and nitrogen concentrations control the decomposition of ectomycorrhizal fungal litter. Soil Biol Biochem. 2014;77:150–7.
CAS  Article  Google Scholar 

58.
Hemkemeyer M, Dohrmann AB, Christensen BT, Tebbe CC. Bacterial preferences for specific soil particle size fractions revealed by community analyses. Front Microbiol. 2018;9:149.
PubMed  PubMed Central  Article  Google Scholar 

59.
Mills A. Keeping in touch: microbial life on soil particle surfaces. Adv Agron. 2003;78:1–43.
Article  Google Scholar 

60.
Kindler R, Miltner A, Thullner M, Richnow H-H, Kästner M. Fate of bacterial biomass derived fatty acids in soil and their contribution to soil organic matter. Org Geochem. 2009;40:29–37.
CAS  Article  Google Scholar 

61.
Huang Y, Liang C, Duan X, Chen H, Li D. Variation of microbial residue contribution to soil organic carbon sequestration following land use change in a subtropical karst region. Geoderma. 2019;353:340–6.
CAS  Article  Google Scholar 

62.
Ahrens B, Braakhekke MC, Guggenberger G, Schrumpf M, Reichstein M. Contribution of sorption, DOC transport and microbial interactions to the 14C age of a soil organic carbon profile: insights from a calibrated process model. Soil Biol Biochem. 2015;88:390–402.
CAS  Article  Google Scholar 

63.
Nguyen RT, Harvey HR. Preservation via macromolecular associations during Botryococcus braunii decay: proteins in the Pula Kerogen. Org Geochem. 2003;34:1391–403.
CAS  Article  Google Scholar 

64.
Kallenbach CM, Frey SD, Grandy AS. Direct evidence for microbial-derived soil organic matter formation and its ecophysiological controls. Nat Commun. 2016;7:13630.
CAS  PubMed  PubMed Central  Article  Google Scholar 

65.
Puget P, Angers DA, Chenu C. Nature of carbohydrates associated with water-stable aggregates of two cultivated soils. Soil Biol Biochem. 1998;31:55–63.
Article  Google Scholar 

66.
Schmidt MWI, Torn MS, Abiven S, Dittmar T, Guggenberger G, Janssens IA, et al. Persistence of soil organic matter as an ecosystem property. Nature. 2011;478:49–56.
CAS  PubMed  Article  Google Scholar 

67.
Spence A, Simpson AJ, McNally DJ, Moran BW, McCaul MV, Hart K, et al. The degradation characteristics of microbial biomass in soil. Geochim Cosmochim Acta. 2011;75:2571–81.
CAS  Article  Google Scholar 

68.
Drigo B, Anderson IC, Kannangara GSK, Cairney JWG, Johnson D. Rapid incorporation of carbon from ectomycorrhizal mycelial necromass into soil fungal communities. Soil Biol Biochem. 2012;49:4–10.
CAS  Article  Google Scholar 

69.
Wang G, Chen S. A review on parameterization and uncertainty in modeling greenhouse gas emissions from soil. Geoderma. 2012;170:206–16.
CAS  Article  Google Scholar 

70.
Blagodatskaya Е, Blagodatsky S, Khomyakov N, Myachina O, Kuzyakov Y. Temperature sensitivity and enzymatic mechanisms of soil organic matter decomposition along an altitudinal gradient on Mount Kilimanjaro. Sci Rep. 2016;6:22240.
CAS  Article  Google Scholar 

71.
German DP, Weintraub MN, Grandy AS, Lauber CL, Rinkes ZL, Allison SD. Optimization of hydrolytic and oxidative enzyme methods for ecosystem studies. Soil Biol Biochem. 2011;43:1387–97.
CAS  Article  Google Scholar 

72.
Wu J, Xiao H. Measuring the gross turnover time of soil microbial biomass C under incubation. Acta Pedol Sin. 2004;41:401–7.
CAS  Google Scholar 

73.
Cheng W. Rhizosphere priming effect: Its functional relationships with microbial turnover, evapotranspiration, and C–N budgets. Soil Biol Biochem. 2009;41:1795–801.
CAS  Article  Google Scholar 

74.
Luo Z, Tang Z, Guo X, Jiang J, Sun OJ. Non-monotonic and distinct temperature responses of respiration of soil microbial functional groups. Soil Biol Biochem. 2020;148:107902.
CAS  Article  Google Scholar 

75.
de Graaff M-A, Classen AT, Castro HF, Schadt CW. Labile soil carbon inputs mediate the soil microbial community composition and plant residue decomposition rates. New Phytol. 2010;188:1055–64.
PubMed  Article  CAS  Google Scholar 

76.
Paul EA. The nature and dynamics of soil organic matter: plant inputs, microbial transformations, and organic matter stabilization. Soil Biol Biochem. 2016;98:109–26.
CAS  Article  Google Scholar 

77.
Crowther TW, Sokol NW, Oldfield EE, Maynard DS, Thomas SM, Bradford MA. Environmental stress response limits microbial necromass contributions to soil organic carbon. Soil Biol Biochem. 2015;85:153–61.
CAS  Article  Google Scholar 

78.
Ding X, Chen S, Zhang B, He H, Filley TR, Horwath WR. Warming yields distinct accumulation patterns of microbial residues in dry and wet alpine grasslands on the Qinghai-Tibetan Plateau. Biol Fertil Soils. 2020;56:881–92.
CAS  Article  Google Scholar 

79.
Mao D, Luo L, Wang Z, Zhang C, Ren C. Variations in net primary productivity and its relationships with warming climate in the permafrost zone of the Tibetan Plateau. J Geogr Sci. 2015;25:967–77.
Article  Google Scholar 

80.
Wu J, Feng Y, Zhang X, Wurst S, Tietjen B, Tarolli P, et al. Grazing exclusion by fencing non-linearly restored the degraded alpine grasslands on the Tibetan Plateau. Sci Rep. 2017;7:15202.
PubMed  PubMed Central  Article  CAS  Google Scholar 

81.
Li J, Wang G, Mayes MA, Allison SD, Frey SD, Shi Z, et al. Reduced carbon use efficiency and increased microbial turnover with soil warming. Glob Change Biol. 2019;25:900–10.
Article  Google Scholar 

82.
Chen G, Ma S, Tian D, Xiao W, Jiang L, Xing A, et al. Patterns and determinants of soil microbial residues from tropical to boreal forests. Soil Biol Biochem. 2020;151:108059.
CAS  Article  Google Scholar 

83.
Wang YP, Chen BC, Wieder WR, Leite M, Medlyn BE, Rasmussen M, et al. Oscillatory behavior of two nonlinear microbial models of soil carbon decomposition. Biogeosciences. 2014;11:1817–31.
CAS  Article  Google Scholar 

84.
Soares M, Rousk J. Microbial growth and carbon use efficiency in soil: links to fungal-bacterial dominance, SOC-quality and stoichiometry. Soil Biol Biochem. 2019;131:195–205.
CAS  Article  Google Scholar 

85.
Liang C, Cheng G, Wixon DL, Balser TC. An Absorbing Markov Chain approach to understanding the microbial role in soil carbon stabilization. Biogeochemistry. 2011;106:303–9.
Article  Google Scholar 

86.
Fan Z, Liang C. Significance of microbial asynchronous anabolism to soil carbon dynamics driven by litter inputs. Sci Rep. 2015;5:9575.
CAS  PubMed  PubMed Central  Article  Google Scholar  More

1.
Montoya, J. M., Pimm, S. L. & Solé, R. V. Ecological networks and their fragility. Nature 442, 259–264 (2006).
CAS  Article  Google Scholar 
2.
Faust, K. & Raes, J. Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538–550 (2012).
CAS  Article  Google Scholar 

3.
Pržulj, N. & Malod-Dognin, N. Network analytics in the age of big data. Science 353, 123–124 (2016).
Article  Google Scholar 

4.
Berry, D. & Widder, S. Deciphering microbial interactions and detecting keystone species with co-occurrence networks. Front. Microbiol. 5, 219 (2014).
Article  Google Scholar 

5.
Okuyama, T. & Holland, J. N. Network structural properties mediate the stability of mutualistic communities. Ecol. Lett. 11, 208–216 (2008).
Article  Google Scholar 

6.
Landi, P., Minoarivelo, H. O., Brännström, Å., Hui, C. & Dieckmann, U. Complexity and stability of ecological networks: a review of the theory. Popul. Ecol. 60, 319–345 (2018).
Article  Google Scholar 

7.
Hillebrand, H. et al. Decomposing multiple dimensions of stability in global change experiments. Ecol. Lett. 21, 21–30 (2018).
Article  Google Scholar 

8.
Toju, H. et al. Species-rich networks and eco-evolutionary synthesis at the metacommunity level. Nat. Ecol. Evol. 1, 0024 (2017).
Article  Google Scholar 

9.
Montesinos-Navarro, A., Hiraldo, F., Tella, J. L. & Blanco, G. Network structure embracing mutualism–antagonism continuums increases community robustness. Nat. Ecol. Evol. 1, 1661–1669 (2017).
Article  Google Scholar 

10.
Ullah, H., Nagelkerken, I., Goldenberg, S. U. & Fordham, D. A. Climate change could drive marine food web collapse through altered trophic flows and cyanobacterial proliferation. PLoS Biol. 16, e2003446 (2018).
Article  CAS  Google Scholar 

11.
Dunne, J. A., Williams, R. J. & Martinez, N. D. Food-web structure and network theory: the role of connectance and size. Proc. Natl Acad. Sci. USA 99, 12917–12922 (2002).
CAS  Article  Google Scholar 

12.
Thébault, E. & Fontaine, C. Stability of ecological communities and the architecture of mutualistic and trophic networks. Science 329, 853–856 (2010).
Article  CAS  Google Scholar 

13.
García-Palacios, P., Gross, N., Gaitán, J. & Maestre, F. T. Climate mediates the biodiversity–ecosystem stability relationship globally. Proc. Natl Acad. Sci. USA 115, 8400–8405 (2018).
Article  CAS  Google Scholar 

14.
IPCC Climate Change 2013: The Physical Science Basis (eds Stocker, T. F. et al.) (Cambridge Univ. Press, 2013).

15.
Xue, K. et al. Tundra soil carbon is vulnerable to rapid microbial decomposition under climate warming. Nat. Clim. Change 6, 595–600 (2016).
CAS  Article  Google Scholar 

16.
Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. & West, G. B. Toward a metabolic theory of ecology. Ecology 85, 1771–1789 (2004).
Article  Google Scholar 

17.
Guo, X. et al. Climate warming leads to divergent succession of grassland microbial communities. Nat. Clim. Change 8, 813–818 (2018).
Article  Google Scholar 

18.
Xu, X., Sherry, R. A., Niu, S., Li, D. & Luo, Y. Net primary productivity and rain-use efficiency as affected by warming, altered precipitation, and clipping in a mixed-grass prairie. Glob. Change Biol. 19, 2753–2764 (2013).
Article  Google Scholar 

19.
Guo, X. et al. Climate warming accelerates temporal scaling of grassland soil microbial biodiversity. Nat. Ecol. Evol. 3, 612–619 (2019).
Article  Google Scholar 

20.
Zhou, J. et al. Functional molecular ecological networks. mBio 1, e00169–10 (2010).
Google Scholar 

21.
Barabási, A.-L. & Oltvai, Z. N. Network biology: understanding the cell’s functional organization. Nat. Rev. Genet. 5, 101–113 (2004).
Article  CAS  Google Scholar 

22.
D’Amen, M., Mod, H. K., Gotelli, N. J. & Guisan, A. Disentangling biotic interactions, environmental filters, and dispersal limitation as drivers of species co-occurrence. Ecography 41, 1233–1244 (2018).
Article  Google Scholar 

23.
Barner, A. K., Coblentz, K. E., Hacker, S. D. & Menge, B. A. Fundamental contradictions among observational and experimental estimates of non-trophic species interactions. Ecology 99, 557–566 (2018).
Article  Google Scholar 

24.
Goberna, M. et al. Incorporating phylogenetic metrics to microbial co-occurrence networks based on amplicon sequences to discern community assembly processes. Mol. Ecol. Resour. 19, 1552–1564 (2019).
Article  Google Scholar 

25.
Carr, A., Diener, C., Baliga, N. S. & Gibbons, S. M. Use and abuse of correlation analyses in microbial ecology. ISME J. 13, 2647–2655 (2019).
Article  Google Scholar 

26.
Morton, J. T. et al. Establishing microbial composition measurement standards with reference frames. Nat. Commun. 10, 2719 (2019).
Article  CAS  Google Scholar 

27.
Fuhrman, J. A. Microbial community structure and its functional implications. Nature 459, 193–199 (2009).
CAS  Article  Google Scholar 

28.
Herren, C. M. & McMahon, K. D. Cohesion: a method for quantifying the connectivity of microbial communities. ISME J. 11, 2426–2438 (2017).
Article  Google Scholar 

29.
Zhou, J., Deng, Y., Luo, F., He, Z. & Yang, Y. Phylogenetic molecular ecological network of soil microbial communities in response to elevated CO2. mBio 2, e00122–11 (2011).
Article  Google Scholar 

30.
Banerjee, S., Schlaeppi, K. & van der Heijden, M. G. A. Keystone taxa as drivers of microbiome structure and functioning. Nat. Rev. Microbiol. 16, 567–576 (2018).
CAS  Article  Google Scholar 

31.
Zelikova, T. J. et al. Long-term exposure to elevated CO2 enhances plant community stability by suppressing dominant plant species in a mixed-grass prairie. Proc. Natl Acad. Sci. USA 111, 15456–15461 (2014).
CAS  Article  Google Scholar 

32.
Douglas, G. M. et al. PICRUSt2 for prediction of metagenome functions. Nat. Biotechnol. 38, 685–688 (2020).
CAS  Article  Google Scholar 

33.
MacArthur, R. Fluctuations of animal populations and a measure of community stability. Ecology 36, 533–536 (1955).
Article  Google Scholar 

34.
May, R. M. Stability and Complexity in Model Ecosystems (Princeton Univ. Press, 2019).

35.
Guo, X. et al. Gene-informed decomposition model predicts lower soil carbon loss due to persistent microbial adaptation to warming. Nat. Commun. 11, 4897 (2020).
CAS  Article  Google Scholar 

36.
Melillo, J. M. et al. Long-term pattern and magnitude of soil carbon feedback to the climate system in a warming world. Science 358, 101–105 (2017).
CAS  Article  Google Scholar 

37.
Zhou, J. et al. Microbial mediation of carbon-cycle feedbacks to climate warming. Nat. Clim. Change 2, 106–110 (2012).
CAS  Article  Google Scholar 

38.
Galiana, N. et al. The spatial scaling of species interaction networks. Nat. Ecol. Evol. 2, 782–790 (2018).
Article  Google Scholar 

39.
Bastolla, U. et al. The architecture of mutualistic networks minimizes competition and increases biodiversity. Nature 458, 1018–1020 (2009).
CAS  Article  Google Scholar 

40.
Cardinale, B. J. et al. Biodiversity loss and its impact on humanity. Nature 486, 59–67 (2012).
CAS  Article  Google Scholar 

41.
Li, D., Zhou, X., Wu, L., Zhou, J. & Luo, Y. Contrasting responses of heterotrophic and autotrophic respiration to experimental warming in a winter annual-dominated prairie. Glob. Change Biol. 19, 3553–3564 (2013).
Google Scholar 

42.
Treves, D. S., Xia, B., Zhou, J. & Tiedje, J. M. A two-species test of the hypothesis that spatial isolation influences microbial diversity in soil. Microb. Ecol. 45, 20–28 (2003).
CAS  Article  Google Scholar 

43.
Zhou, J., Xia, B., Huang, H., Palumbo, A. V. & Tiedje, J. M. Microbial diversity and heterogeneity in sandy subsurface soils. Appl. Environ. Microbiol. 70, 1723–1734 (2004).
CAS  Article  Google Scholar 

44.
Zhou, J. et al. Spatial and resource factors influencing high microbial diversity in soil. Appl. Environ. Microbiol. 68, 326–334 (2002).
CAS  Article  Google Scholar 

45.
O’Brien, S. L. et al. Spatial scale drives patterns in soil bacterial diversity. Environ. Microbiol. 18, 2039–2051 (2016).
Article  Google Scholar 

46.
Penton, C. R., Gupta, V. V. S. R., Yu, J. & Tiedje, J. M. Size matters: assessing optimum soil sample size for fungal and bacterial community structure analyses using high throughput sequencing of rRNA gene amplicons. Front. Microbiol. 7, 824 (2016).
Google Scholar 

47.
Zhou, J., Bruns, M. A. & Tiedje, J. M. DNA recovery from soils of diverse composition. Appl. Environ. Microbiol. 62, 316–322 (1996).
CAS  Article  Google Scholar 

48.
Hurt, R. A. et al. Simultaneous recovery of RNA and DNA from soils and sediments. Appl. Environ. Microbiol. 67, 4495–4503 (2001).
CAS  Article  Google Scholar 

49.
Peiffer, J. A. et al. Diversity and heritability of the maize rhizosphere microbiome under field conditions. Proc. Natl Acad. Sci. USA 110, 6548–6553 (2013).
CAS  Article  Google Scholar 

50.
Wu, L. et al. Phasing amplicon sequencing on Illumina Miseq for robust environmental microbial community analysis. BMC Microbiol. 15, 125 (2015).
Article  CAS  Google Scholar 

51.
Wen, C. et al. Evaluation of the reproducibility of amplicon sequencing with Illumina MiSeq platform. PLoS ONE 12, e0176716 (2017).
Article  CAS  Google Scholar 

52.
Zhou, J. et al. High-throughput metagenomic technologies for complex microbial community analysis: open and closed formats. mBio 6, e02288–14 (2015).
CAS  Article  Google Scholar 

53.
Zhou, J. et al. Reproducibility and quantitation of amplicon sequencing-based detection. ISME J. 5, 1303–1313 (2011).
CAS  Article  Google Scholar 

54.
Luo, F. et al. Constructing gene co-expression networks and predicting functions of unknown genes by random matrix theory. BMC Bioinformatics 8, 299 (2007).
Article  CAS  Google Scholar 

55.
Luo, F., Zhong, J., Yang, Y., Scheuermann, R. H. & Zhou, J. Application of random matrix theory to biological networks. Phys. Lett. A 357, 420–423 (2006).
CAS  Article  Google Scholar 

56.
Deng, Y. et al. Molecular ecological network analyses. BMC Bioinformatics 13, 113 (2012).
Article  Google Scholar 

57.
Shi, S. et al. The interconnected rhizosphere: high network complexity dominates rhizosphere assemblages. Ecol. Lett. 19, 926–936 (2016).
Article  Google Scholar 

58.
Mehta, M. L. Random Matrices 2nd edn (Elsevier, 2004).

59.
Plerou, V., Gopikrishnan, P., Rosenow, B., Amaral, L. A. N. & Stanley, H. E. Universal and non-universal properties of cross-correlations in financial time series. Phys. Rev. Lett. 83, 1471–1474 (1999).
CAS  Article  Google Scholar 

60.
Aitchison, J. The statistical analysis of compositional data. J. R. Stat. Soc. B 44, 139–160 (1982).
Google Scholar 

61.
Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J. Microbiome datasets are compositional: and this is not optional. Front. Microbiol. 8, 2224 (2017).
Article  Google Scholar 

62.
Pawlowsky-Glahn, V. & Egozcue, J. J. Compositional data and their analysis: an introduction. Geol. Soc. Spec. Publ. 264, 1–10 (2006).
CAS  Article  Google Scholar 

63.
Friedman, J. & Alm, E. J. Inferring correlation networks from genomic survey data. PLoS Comput. Biol. 8, e1002687 (2012).
CAS  Article  Google Scholar 

64.
Watts, S. C., Ritchie, S. C., Inouye, M. & Holt, K. E. FastSpar: rapid and scalable correlation estimation for compositional data. Bioinformatics 35, 1064–1066 (2019).
CAS  Article  Google Scholar 

65.
Weiss, S. et al. Correlation detection strategies in microbial data sets vary widely in sensitivity and precision. ISME J. 10, 1669–1681 (2016).
CAS  Article  Google Scholar 

66.
R: a language and environment for statistical computing (R Foundation for Statistical Computing, 2019).

67.
Goslee, S. C. & Urban, D. L. The ecodist package for dissimilarity-based analysis of ecological data. J. Stat. Softw. 22, 1–19 (2007).
Article  Google Scholar 

68.
Oksanen, J. et al. vegan: Community Ecology Package. Version 2.5-6 (2019).

69.
Lima-Mendez, G. et al. Determinants of community structure in the global plankton interactome. Science 348, 1262073 (2015).
Article  CAS  Google Scholar 

70.
Yuan, M.M. et al. Mengting-Maggie-Yuan/warming-network-complexity-stability: warming-network-complexity-stability-v1.0. Version 1.0 (Zenodo, 2021); https://doi.org/10.5281/zenodo.4383469

71.
He, Z. et al. GeoChip 3.0 as a high-throughput tool for analyzing microbial community composition, structure and functional activity. ISME J. 4, 1167–1179 (2010).
CAS  Article  Google Scholar 

72.
He, Z. et al. GeoChip: a comprehensive microarray for investigating biogeochemical, ecological and environmental processes. ISME J. 1, 67–77 (2007).
CAS  Article  Google Scholar 

73.
Ning, D., Deng, Y., Tiedje, J. M. & Zhou, J. A general framework for quantitatively assessing ecological stochasticity. Proc. Natl Acad. Sci. USA 116, 16892–16898 (2019).
CAS  Article  Google Scholar 

74.
Zhou, J. & Ning, D. Stochastic community assembly: does it matter in microbial ecology? Microbiol. Mol. Biol. Rev. 81, e00002–e00017 (2017).
Article  Google Scholar 

75.
Csárdi, G. & Nepusz, T. The igraph software package for complex network research. InterJ. Complex Syst. 1695, 1–9 (2006).
Google Scholar 

76.
Maslov, S. & Sneppen, K. Specificity and stability in topology of protein networks. Science 296, 910–913 (2002).
CAS  Article  Google Scholar 

77.
Almeida‐Neto, M., Guimarães, P., Guimarães, P. R., Loyola, R. D. & Ulrich, W. A consistent metric for nestedness analysis in ecological systems: reconciling concept and measurement. Oikos 117, 1227–1239 (2008).
Article  Google Scholar 

78.
Guimerà, R. & Nunes Amaral, L. A. Functional cartography of complex metabolic networks. Nature 433, 895–900 (2005).
Article  CAS  Google Scholar 

79.
Olesen, J. M., Bascompte, J., Dupont, Y. L. & Jordano, P. The modularity of pollination networks. Proc. Natl Acad. Sci. USA 104, 19891–19896 (2007).
CAS  Article  Google Scholar 

80.
Banerjee, S., Schlaeppi, K. & van der Heijden, M. G. A. Reply to ‘can we predict microbial keystones?’. Nat. Rev. Microbiol. 17, 194 (2019).
CAS  Article  Google Scholar 

81.
Röttjers, L. & Faust, K. Can we predict keystones? Nat. Rev. Microbiol. 17, 193 (2019).
Article  CAS  Google Scholar 

82.
Langfelder, P. & Horvath, S. Eigengene networks for studying the relationships between co-expression modules. BMC Syst. Biol. 1, 54 (2007).
Article  CAS  Google Scholar 

83.
Hautier, Y. et al. Eutrophication weakens stabilizing effects of diversity in natural grasslands. Nature 508, 521–525 (2014).
CAS  Article  Google Scholar 

84.
Hui, C., McGeoch, M. A., Harrison, A. E. S. & Bronstein, E. J. L. Zeta diversity as a concept and metric that unifies incidence-based biodiversity patterns. Am. Nat. 184, 684–694 (2014).
Article  Google Scholar 

85.
Shi, Z. et al. Functional gene array-based ultrasensitive and quantitative detection of microbial populations in complex communities. mSystems 4, e00296–19 (2019).
Google Scholar 

86.
Sun, S., Jones, R. B. & Fodor, A. A. Inference-based accuracy of metagenome prediction tools varies across sample types and functional categories. Microbiome 8, 46 (2020).
Article  Google Scholar  More