OSCR

Evolutionary stasis in synapsid encephalization during the end-permian mass extinction.

Code ↔ Paper

6 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 6 matches
  1. [1] § Phylogenetic comparative analyses ↔ Tree_calibration_Benoitetal_2026.Rmd, lines 76–120 · score 0.84 · consensus.edges, multi2di, extinction rates, calibrated trees, durationfreqcont, cal3
  2. [2] § Phylogenetic comparative analyses ↔ Script_Benoitetal_2026.Rmd, lines 1002–1041 · score 0.75 · RRphylo, terminal taxa, Biarmosuchia, Dinocephalia, Gorgonopsia, Therocephalia
  3. [3] § Phylogenetic comparative analyses ↔ Script_Benoitetal_2026.Rmd, lines 213–260 · score 0.67 · Akaike Information Criterion, finite sample, evolutionary models, Pagel, body mass, Brownian
  4. [4] § Phylogenetic comparative analyses ↔ Script_Benoitetal_2026.Rmd, lines 85–108 · score 0.63 · lambda model, best fit, Phylopars, imputation, imputed, BIC
  5. [5] § Phylogenetic comparative analyses ↔ Script_Benoitetal_2026.Rmd, lines 262–307 · score 0.62 · extremely low, Phylogenetic signal, weight, body mass, coefficient, AICc
  6. [6] § Phylogenetic comparative analyses ↔ Script_Benoitetal_2026.Rmd, lines 213–260 · score 0.62 · variance covariance, evolutionary models, PGLS, body mass, coefficient, AICc

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R Markdown · 1,096 lines · 50 KB · no license · 5 matches

  1. ---
  2. title: "Script Benoit et al. 2026"
  3. author: "Lucas Legendre"
  4. date: "`r Sys.Date()`"
  5. output:
  6. html_document:
  7. df_print: paged
  8. toc_float:
  9. collapsed: false
  10. smooth_scroll: false
  11. toc: true
  12. ---
  14. ```{r setup, include=FALSE}
  15. knitr::opts_chunk$set(echo = TRUE)
  16. ```
  18. Compiled under R version 4.5.3 (2026-03-11)
  20. <b>WARNING</b>: edit the working directory to your preferred folder.
  22. This document details all analyses performed in R for the study:
  23. Benoit J., Legendre L.J., Midzuk A., Fernandez V., Araújo R., Browning C. Abdala F., Botha J., Angielczyk K. 2026. Evolutionary stasis in synapsid encephalization during the End-Permian mass extinction. *Scientific Reports*.
  25. For more information regarding the study, datasets, and analyses, please refer to the Main Text and Supplementary Information of this paper. If you have any additional questions, feel free to email me at [lucasjlegendre\@gmail.com](mailto:[email hidden]){.email}.
  27. ### Loading packages
  29. ```{r, message = FALSE}
  30. library(ape)
  31. library(nlme)
  32. library(phytools)
  33. library(phylolm)
  34. library(geiger)
  35. library(windex)
  36. library(AICcmodavg)
  37. library(evobiR)
  38. library(tidyverse)
  39. library(RRphylo)
  40. library(Rphylopars)
  41. library(hrbrthemes)
  42. library(RColorBrewer)
  43. library(strap)
  44. library(slouch)
  45. library(scales)
  46. ```
  48. ## Seed and cores
  49. ```{r}
  50. set.seed(101)
  51. cc <- max(1, floor(0.5 * parallel::detectCores()))
  52. ```
  54. ## Data and tree
  56. - Data
  57. ```{r}
  58. maindata<-read.csv("maindataNew.csv", header=TRUE, sep=";")
  59. ```
  61. - Tree
  62. ```{r}
  63. # Calibrated tree (see separate tree calibration script)
  64. ttree<-read.nexus("CalibratedTreeNew.nex")
  65. plotTree(ttree, fsize=0.5, lwd=1, ftype="i")
  66. ttree$root.time<-263
  68. setdiff(ttree$tip.label, maindata$Taxon)
  69. setdiff(maindata$Taxon, ttree$tip.label)
  70. # Perfect match between tip labels and taxa in the dataset
  72. # Reorder data
  73. maindata2<-ReorderData(ttree, maindata, taxa.names=1)
  74. rownames(maindata2)<-maindata2$Taxon
  75. ```
  77. ## Preliminary analyses – evolutionary models, phenograms, ancestral reconstructions
  79. - Log-conversion
  80. ```{r}
  81. logmain<-maindata2[,c("Taxon", "Order", "Body_mass", "EQ_noOB", "EQ_OB")]
  82. logmain[,3:5]<-log1p(logmain[,3:5]) # log1p to get only positive values
  83. ```
  85. - Predict missing data in the dataset, using `Rphylopars` (Goolsby et al., 2017)
  86. ```{r}
  87. # Data
  88. parsdata<-logmain[,c(1,4,5)]
  89. colnames(parsdata)[1]<-"species"
  91. # Imputation
  92. modelpars<-c("BM", "lambda", "EB", "star") # Evolutionary models
  93. PPE<-list(); BICdat<-data.frame()
  94. for (i in 1:length(modelpars)) {
  95. PPE[[i]]<-phylopars(trait_data = parsdata, tree = ttree, model=modelpars[i])
  96. BICdat[i,1]<-BIC(PPE[[i]]) # BIC to compare models
  97. }
  98. rownames(BICdat)<-modelpars; colnames(BICdat)<-"BIC"
  99. BICdat # The lambda model is the best fit
  101. # Extract imputed values and add them to the main dataset
  102. imp <- PPE[[which.min(BICdat$BIC)]]$anc_recon[1:nrow(logmain),
  103. c("EQ_noOB","EQ_OB")]
  104. for(tr in c("EQ_noOB","EQ_OB")) {
  105. idx_na <- is.na(logmain[[tr]])
  106. logmain[[tr]][idx_na] <- imp[rownames(logmain)[idx_na], tr]
  107. }
  108. ```
  110. - Phylogenetic signal (Pagel's lambda) in `phytools` (Revell, 2024)
  111. ```{r}
  112. # Phylogenetic signal (Pagel's lambda)
  113. var=list(); phy=list(); phydat=data.frame()
  114. for (i in 3:5) {
  115. var<-logmain[,i]; names(var)<-rownames(logmain)
  116. var2<-var[!is.na(var)]
  117. treevar<-drop.tip(ttree,setdiff(ttree$tip.label,names(var2)))
  118. phy[[i]]<-phylosig(treevar, var2, method="lambda", test=TRUE)
  119. phydat[i,1]<-colnames(logmain)[i]; phydat[i,2]<-phy[[i]]$lambda; phydat[i,3]<-phy[[i]]$P
  120. }
  121. colnames(phydat)<-c("Variable", "Lambda", "P")
  122. phydat[c(3:5),]
  123. ```
  124. Signal is very high for body mass (λ > 0.99), high for EQ with and without olfactory bulbs (λ = 0.73 and λ = 0.70, respectively), and highly significant (p < 0.001) for all three traits.
  126. - Test for best evolutionary model, using `geiger` (Pennell et al., 2014)
  127. ```{r, warning = FALSE}
  128. models=c("BM", "OU", "EB", "rate_trend","lambda", "white")
  129. # (For more info, check '?fitContinuous')
  131. var=list(); fit=list(); mod=list()
  132. for (i in 3:5) {
  133. var<-logmain[,i]; names(var)<-rownames(logmain)
  134. var2<-var[!is.na(var)]
  135. treevar<-multi2di(drop.tip(ttree,setdiff(ttree$tip.label,names(var2))))
  136. for (m in 1:length(models)) {
  137. fit[[m]]=fitContinuous(treevar, var2, model=models[m], ncores=2)
  138. }
  139. mod[[i]]<-modSel.geiger(fit[[1]],fit[[2]],fit[[3]],fit[[4]],fit[[5]],fit[[6]])
  140. }
  141. mod
  142. ```
  143. Best model for body mass is an Early Burst model (EB). A Lambda model fits both types of EQ best – both traits follow a Brownian motion model; we can reconstruct ancestral states using such a model and map them on the tree.
  145. - Phenogram for each trait
  146. ```{r}
  147. for (i in 3:5) {
  148. trait<-logmain[,i]; names(trait)<-rownames(logmain)
  149. trait<-na.omit(trait)
  150. treena<-drop.tip(ttree, setdiff(ttree$tip.label, names(trait)))
  151. phenogram(treena, trait, ,ftype="i", spread.labels=TRUE,
  152. spread.cost=c(1,0),fsize=0.7,color=palette()[4],
  153. xlab="Time (Ma)",ylab=paste("log",colnames(logmain)[i]),cex.axis=0.8,
  154. axes=FALSE,las=1, lwd=1)
  155. axis(1,at=round(seq(0,max(nodeHeights(treena)),length.out=6),1),
  156. label=round(seq(max(nodeHeights(treena)),0,length.out=6),1),
  157. cex.axis=0.8)
  158. axis(2,las=1,cex.axis=0.8)
  159. grid()
  160. title(paste(colnames(logmain[i])))
  161. }
  162. ```
  164. Some clear shifts in all traits, but they seem to occur in Mammaliaformes and less inclusive clades.
  166. - Ancestral state reconstructions
  167. ```{r}
  168. # Custom color palette
  169. pal<-c("yellow1","orange1","orangered3","purple4")
  171. # Compile the best tree with corresponding value of lambda for both types of EQ
  172. lambdaval1<-phydat[4,2]; lambdaval2<-phydat[5,2]
  173. treeEQnoOB<-ttree; treeEQOB<-ttree
  174. treeEQnoOB<-rescale(treeEQnoOB, "lambda", lambdaval1)
  175. treeEQOB<-rescale(treeEQOB, "lambda", lambdaval2)
  176. treelist<-c(ttree, treeEQnoOB, treeEQOB)
  178. dataplot=list(); fit=list(); obj=list()
  179. for (i in c(3:5)) {
  180. dataplot[[i]]<-logmain[,i]; names(dataplot[[i]])<-rownames(logmain)
  181. dataplot[[i]]<-na.omit(dataplot[[i]])
  182. treeplot<-treelist[[i-2]]
  183. fit[[i]]<-fastAnc(treeplot, dataplot[[i]], vars=TRUE, CI=TRUE)
  184. obj[[i]]<-setMap(contMap(treeplot, dataplot[[i]], plot=FALSE),
  185. colors=pal)
  186. plot(obj[[i]], fsize=0.5, ftype="i")
  187. title(paste('Ancestral state reconstruction for', colnames(logmain)[i]))
  188. }
  189. ```
  191. As expected from phenograms, we do not see clear shifts in any major therapsid clade outside Mammaliaformes. The pattern in EQ seems to be the exact reverse of that of body mass – as lineages get smaller, their EQ increases. We can test that further for specific lineages and see whether there are any specific shifts in the tree.
  193. ### How well correlated to body mass is the EQ?
  195. We can test it using Phylogenetic Generalized Least Squares (PGLS) regressions.
  197. - Import function `fitEvolPar` to estimate the best fit of the main parameter (alpha, lambda, or g) in the corresponding evolutionary model (OU, lambda, or EB, respectively)
  199. (To download it from GitHub: <https://github.com/LucasLegendre/fitEvolPar>)
  201. ```{r}
  202. sys.source("fitEvolPar.R", envir = knitr::knit_global())
  203. ```
  205. ```{r}
  206. # Compile the VCV matrix (necessary for a non-ultrametric tree, which is always going to be the case with fossils in the sample)
  207. Wt<-diag(vcv.phylo(ttree))
  209. # Data frame for 'fitEvolPar'
  210. EQnoOBd<-logmain[,3:4]
  211. ```
  213. #### Compile individual regressions
  215. We also test for best evolutionary model in the variance-covariance (VCV) matrix in the error term, using the Akaike Information Criterion corrected for finite sample sizes (AICc).
  217. - For EQ with no OB
  218. ```{r, warning=FALSE}
  219. # PGLS models
  220. BM<-gls(EQ_noOB~Body_mass, data=logmain,
  221. correlation=corBrownian(phy=ttree, form=~1),
  222. weights=varFixed(~Wt), method="ML")
  223. OU<-gls(EQ_noOB~Body_mass, data=logmain,
  224. correlation=corMartins(0.1, phy=ttree, form=~1),
  225. weights=varFixed(~Wt), method="ML")
  226. Lambda<-gls(EQ_noOB~Body_mass, data=logmain,
  227. correlation=corPagel(1, phy=ttree, form=~1),
  228. weights=varFixed(~Wt), method="ML")
  229. EB<-gls(EQ_noOB~Body_mass, data=logmain,
  230. correlation=corBlomberg(fitEvolPar(EQnoOBd, ttree,"EB"),
  231. phy=ttree, fixed = TRUE, form=~1),
  232. weights=varFixed(~Wt), method="ML")
  233. OLS<-gls(EQ_noOB~Body_mass, data=logmain, method="ML")
  235. # AICc
  236. Cand.models = list()
  237. Cand.models[[1]] = BM
  238. Cand.models[[2]] = OU
  239. Cand.models[[3]] = Lambda
  240. Cand.models[[4]] = EB
  241. Cand.models[[5]] = OLS
  243. Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
  244. aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
  246. # Best model
  247. best1<-lm(EQ_noOB~Body_mass, data=logmain)
  248. summary(best1)
  250. # Plot
  251. ggplot(logmain, aes(Body_mass, EQ_noOB, color=Order)) +
  252. geom_point(size=5) +
  253. # geom_text(aes(label=Taxon), hjust=-0.1, vjust=0.4) +
  254. xlab("ln body mass (g)") +
  255. ylab("ln EQ (no OB)") +
  256. geom_abline(intercept=best1$coefficients[1], slope=best1$coefficients[2],
  257. colour="royalblue", linewidth=1.3) +
  258. scale_color_brewer(palette="Dark2") +
  259. theme_ipsum(axis_title_size=15, base_size = 20, axis_text_size = 12)
  260. ```
  262. Best model is a white model (i.e. no phylogenetic signal). The model is significant (p < 0.01), but the pseudo R-squared is extremely low (R^2 = 0.05). Looking at the plot, there does indeed not seem to be any clear effect of body mass on the EQ.
  264. - For EQ with OB
  265. ```{r, warning=FALSE}
  266. # PGLS models
  267. BM<-gls(EQ_OB~Body_mass, data=logmain,
  268. correlation=corBrownian(phy=ttree, form=~1),
  269. weights=varFixed(~Wt), method="ML")
  270. OU<-gls(EQ_OB~Body_mass, data=logmain,
  271. correlation=corMartins(0.1, phy=ttree, form=~1),
  272. weights=varFixed(~Wt), method="ML")
  273. Lambda<-gls(EQ_OB~Body_mass, data=logmain,
  274. correlation=corPagel(1, phy=ttree, form=~1),
  275. weights=varFixed(~Wt), method="ML")
  276. EB<-gls(EQ_OB~Body_mass, data=logmain,
  277. correlation=corBlomberg(fitEvolPar(EQnoOBd, ttree,"EB"),
  278. phy=ttree, fixed = TRUE, form=~1),
  279. weights=varFixed(~Wt), method="ML")
  280. OLS<-gls(EQ_OB~Body_mass, data=logmain, method="ML")
  282. # AICc
  283. Cand.models = list()
  284. Cand.models[[1]] = BM
  285. Cand.models[[2]] = OU
  286. Cand.models[[3]] = Lambda
  287. Cand.models[[4]] = EB
  288. Cand.models[[5]] = OLS
  290. Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
  291. aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
  293. # Best model
  294. best2<-lm(EQ_OB~Body_mass, data=logmain)
  295. summary(best2)
  297. # Plot
  298. ggplot(logmain, aes(Body_mass, EQ_OB, color=Order)) +
  299. geom_point(size=5) +
  300. # geom_text(aes(label=Taxon), hjust=-0.1, vjust=0.4) +
  301. xlab("ln body mass (g)") +
  302. ylab("ln EQ (OB)") +
  303. geom_abline(intercept=best2$coefficients[1], slope=best2$coefficients[2],
  304. colour="royalblue", linewidth=1.3) +
  305. scale_color_brewer(palette="Dark2") +
  306. theme_ipsum(axis_title_size=15, base_size = 20, axis_text_size = 12)
  307. ```
  309. Similar result as without OB: the regression is significant (p < 0.01), but R-squared is very low (R^2 = 0.07), and there is clearly no signal in the plot.
  311. ## Evolutionary rates and testing for shifts, using `RRphylo`
  313. More info on the method (phylogenetic ridge regressions) and how to use it in the very detailed help section and tutorials of this package; see also Kratsch & McHardy (2014) and Castiglione et al. (2018, 2019).
  315. Since we did not recover body mass as a good predictor of EQ, we will not use it here as a covariate to look at the evolution of EQ and test for shifts in values and/or evolutionary rates.
  317. ### For EQ with no OB
  319. - Prepare the data and tree
  320. ```{r}
  321. # Dataset
  322. EQ1<-logmain[,4]; names(EQ1)<-rownames(logmain)
  324. # Proportion of clusters
  325. cc <- 2/parallel::detectCores()
  327. # Compile RRphylo
  328. RREQ1<-RRphylo(tree=ttree, y=EQ1, clus=cc)
  330. # To visualize node numbers in the tree
  331. nodetree<-ttree; nodetree$edge.length<-rep(1, 250) # Version of the tree where all branches have a length of 1 (easier to visualize)
  332. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  333. labelnodes(text=1:ttree$Nnode,node=1:ttree$Nnode+Ntip(ttree),
  334. interactive=FALSE,circle.exp=0.4,cex=0.5)
  335. ```
  337. #### Let's test for significant shifts in evolutionary rates in the tree, using `search.shift`
  339. - In the whole tree (as a first step, to see if any main trends can be identified)
  340. ```{r}
  341. shiftsWholeR<-search.shift(RREQ1, status.type= "clade")
  342. shiftsWholeR$all.clades
  343. as.numeric(rownames(shiftsWholeR$all.clades))-length(ttree$tip.label) # actual node numbers as visualized in plotTree
  345. # Significant shifts
  346. posshift<-as.numeric(names(which(shiftsWholeR$all.clades[,2]>0.975)))
  347. posshift # Significant positive shifts
  348. negshift<-as.numeric(names(which(shiftsWholeR$all.clades[,2]<0.025)))
  349. negshift # Significant negative shifts
  351. # Plot them on the tree
  352. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  353. labelnodes(text=1:ttree$Nnode,node=1:ttree$Nnode+Ntip(ttree),
  354. interactive=FALSE,circle.exp=0.4,cex=0.5)
  355. nodelabels(text="+", node=posshift, adj=c(1.5, 1), frame="none", col="red2")
  356. nodelabels(text="-", node=negshift, adj=c(1.5, 1), frame="none", col="steelblue2")
  357. ```
  359. **Two major series of shifts, both in Eucynodontia**: higher rates in Cynognathia (suggesting a selective pressure on EQ in that clade), and lower rates in Probainognathia (suggesting relaxed constraints, potentially associated with ancestral stabilizing selection). Eucynodontia and the immediately more inclusive clade (with *Platycraniellus*) also show a shift towards a lower rate, suggesting **the trend in Probainognathia is ancestral to eucynodonts**, whereas **that in Cynognathia is a later, distinct shift**.
  361. - Let's make a phylomorphospace for EQ with no OB and body mass, with colors for these two clades:
  362. ```{r}
  363. # Data
  364. EQnoOBdat<-logmain[,3:4]
  366. # Basic phylomorphospace, if you need it – remove the 'label="off"' part if you want taxa names, but it gets very crowded!
  367. par(mar = c(5.1, 4.1, 4.1, 2.1))
  368. phylomorphospace(ttree,EQnoOBdat,bty="l",label="off",xlab="log (Body mass) (g)",
  369. ylab="log (EQ) (no OB)",node.size=c(0,1))
  370. title(main="Phylomorphospace plot",font.main=3)
  372. # Basic tree with labelled clades
  373. plotTree(ttree,fsize=0.5,ftype="i")
  374. nodelabels(frame="circ",bg="white",cex=0.3)
  375. cladelabels(ttree,c("Probainognathia","Cynognathia"),c(144,179))
  377. # Define colors for the tree
  378. painted<-paintSubTree(ttree, 144, "Probainognathia")
  379. painted<-paintSubTree(painted, 179, "Cynognathia")
  381. # Plot it on the tree
  382. plot(painted, fsize=0.5, ftype="i")
  383. cladelabels(ttree, c("Probainognathia","Cynognathia"), c(144,179), offset=0.5)
  385. # Plot it on the phylomorphospace
  386. par(mar = c(5.1, 4.1, 4.1, 2.1))
  387. phylomorphospace(painted,EQnoOBdat,bty="l",label="off",xlab="log (Body mass) (g)",
  388. ylab="log (EQ) (no OB)",node.size=c(0,1.2),node.by.map=TRUE)
  389. title(main="Phylomorphospace plot",font.main=3)
  390. legend(x="topleft",legend=c("Cynognathia","Probainognathia"),
  391. pch=21,pt.cex=1.5,pt.bg=palette()[2:4],bty="n")
  393. # Plot it on the phenogram
  394. EQnoOB<-logmain[,4]; names(EQnoOB)<-rownames(logmain)
  395. par(mar = c(5.1, 4.1, 4.1, 2.1))
  396. phenogram(painted, EQnoOB, ,ftype="off", spread.labels=TRUE,
  397. spread.cost=c(1,0),fsize=0.7,
  398. xlab="Time (Ma)",ylab="log (EQ) (no OB)",cex.axis=0.8,
  399. axes=FALSE,las=1, lwd=1)
  400. axis(1,at=round(seq(0,max(nodeHeights(painted)),length.out=6),1),
  401. label=round(seq(max(nodeHeights(painted)),0,length.out=6),1),
  402. cex.axis=0.8)
  403. axis(2,las=1,cex.axis=0.8)
  404. grid()
  405. title(paste("Phenogram for EQ (no OB)"))
  406. legend(x="bottomright",legend=c("Cynognathia","Probainognathia"),
  407. pch=21,pt.cex=1.5,pt.bg=palette()[2:4],bty="n")
  408. ```
  410. As expected, Cynognathia show very high variation in a short amount of time, while the earliest nodes within Probainognathia seem to show limited variation over an extended time frame compared to more inclusive clades in the tree.
  412. Let's go a bit further... Same test, but for specific clades in the tree!
  414. We will be testing for shifts in the following clades:
  415. Biarmosuchia, Dinocephalia, Neotherapsida, Dicynodontia, Theriodontia, Gorgonopsia, Eutheriodontia, Therocephalia, Cynodontia, Eucynodontia, Cynognathia, Probainognathia, Mammaliamorpha, Mammaliaformes, Mammalia.
  417. - Specific shifts in those clades
  418. ```{r, warning=FALSE}
  419. # Node numbers
  420. nodes<-c(124,120,3,93,4,90,5,83,6,17,53,18,26,28,32)
  421. cornodes<-nodes+length(ttree$tip.label)
  423. # Look for significant shifts
  424. shiftsR<-search.shift(RREQ1, status.type= "clade", node=cornodes)
  425. shiftsR$single.clades
  426. which(shiftsR$single.clades[,2]>0.975); which(shiftsR$single.clades[,2]<0.025)
  428. posshift<-as.numeric(names(which(shiftsR$single.clades[,2]>0.975)))
  429. posshift # Significant positive shifts
  430. negshift<-as.numeric(names(which(shiftsR$single.clades[,2]<0.025)))
  431. negshift # Significant negative shifts
  433. # Plot them on the tree
  434. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  435. labelnodes(text=1:ttree$Nnode,node=1:ttree$Nnode+Ntip(ttree),
  436. interactive=FALSE,circle.exp=0.4,cex=0.5)
  437. #nodelabels(text="+", node=posshift, adj=c(1.5, 1), frame="none", col="red2")
  438. nodelabels(text="-", node=negshift, adj=c(1.5, 1), frame="none", col="steelblue2")
  439. ```
  441. Only the ancestral negative shift in Cynodontia, followed by similar shifts in Probainognathia, Mammaliamorpha, Mammaliaformes, and Mammalia are recovered.
  443. - We can plot the shifts on the tree
  444. ```{r}
  445. RRplotR<-plotRR(RREQ1, y=EQnoOB, multivariate = "rates")
  447. ## Recompile individual shifts for the plot
  448. posshift<-as.numeric(names(which(shiftsWholeR$all.clades[,2]>0.975)))
  449. negshift<-as.numeric(names(which(shiftsWholeR$all.clades[,2]<0.025)))
  450. EQ1shift<-search.shift(RREQ1, status.type= "clade", node=c(posshift, negshift))
  451. EQ1shift$single.clades<-as.data.frame(EQ1shift$single.clades)
  452. which(EQ1shift$single.clades[,2]>0.975); which(EQ1shift$single.clades[,2]<0.025)
  454. ## Values of EQ
  455. RRplotR$plotRRphen(variable=1, tree.args=list(no.margin=TRUE),
  456. colorbar.args=list(x=170, y=90, title.pos="bottom"))
  458. ## Rates and rate shifts on the same plot
  459. RRplotR$plotRRrates(variable=1, tree.args=list(no.margin=TRUE),
  460. colorbar.args=list(x=170, y=90, title.pos="bottom"))
  461. addShift(EQ1shift, symbols.args=list(lwd=2,fg=c(pos="red2",neg="steelblue2"),
  462. bg=scales::alpha(c(rep("red2",length(posshift)),
  463. rep("steelblue2",length(negshift))),0.3)))
  465. ## Only rate shifts
  466. plotShiftR<-plotShift(RREQ1, EQ1shift)
  467. plotShiftR$plotClades(tree.args=list(no.margin=TRUE),
  468. symbols.args=list(lwd=2,fg=c(pos="red2",neg="steelblue2"),
  469. bg=scales::alpha(c(rep("red2",length(posshift)),
  470. rep("steelblue2",length(negshift))),0.3)))
  471. ```
  473. - Finally, let's test for the significance of those shifts with jacknife/multiple permutation tests, using `overfitRR`. We use a threshold of 0.75 (3/4 of simulated trees recovered with the same shift) to assess a shift as significant.
  474. ```{r}
  475. # Generate alternative topologies
  476. resampleEQ1<-resampleTree(RREQ1$tree, node=c(posshift, negshift), nsim=1000)
  478. # Recompile regressions for the new trees
  479. overfitEQ1<-overfitRR(RREQ1, y=EQnoOB, phylo.list=resampleEQ1, clus=cc)
  481. # Test for shift significance with new topologies
  482. overshift1<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=148)
  483. overshift2<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=c(186,143))
  484. overshift3<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=c(150,144))
  485. overshift4<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=c(151,145))
  486. overshift5<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=c(152,142))
  487. overshift6<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=c(153,147))
  488. overshift7<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=c(154,185))
  489. overshift8<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=c(155,149))
  490. overshift9<-overfitSS(RR=RREQ1, oveRR=overfitEQ1, node=c(156,146))
  492. # Display the results
  493. overshift1$shift.results$clade
  494. overshift2$shift.results$clade
  495. overshift3$shift.results$clade
  496. overshift4$shift.results$clade
  497. overshift5$shift.results$clade
  498. overshift6$shift.results$clade
  499. overshift7$shift.results$clade
  500. overshift8$shift.results$clade
  501. overshift9$shift.results$clade
  503. ```
  504. Nodes 16 through 30 (i.e. Probainognathia through Mammaliaformes) all have significant negative shifts, whereas all other nodes (i.e. clades within Cynognathia) do not.
  505. => **Only the rate shifts in Probainognathia (negative shifts, i.e. stabilizing selection of EQ) are recovered as significant.**
  507. Now that we have identified shifts in evolutionary rates, what about phenotypic trends (i.e. changes in the actual values of EQ through time)?
  509. #### Testing for significant phenotypic trends for EQ, using `search.trend`
  511. See the vignette and help section of the function, as well as Castiglione et al. (2019). We are testing for significant deviations from a Brownian Motion (BM) model in the value of our trait of interest (EQ).
  513. - For the whole tree
  514. ```{r}
  515. ST<-search.trend(RR=RREQ1, y=EQnoOB, nsim=1000)
  517. # Display the results
  518. rateres<-cbind(ST$rate.regression[1,1:4], NA); colnames(rateres)[5]<-"dev"
  519. phenres<-c(ST$phenotypic.regression[1,1:3], NA, ST$phenotypic.regression[1,4])
  520. names(phenres)[c(4,5)]<-c("spread", "dev")
  521. STres<-rbind(rateres, phenres)
  522. rownames(STres)<-c("rescaled absolute rate regression", "phenotypic regression")
  523. STres
  524. ```
  525. The absolute rate regression is not significant (p.random < 0.95), meaning the rates for the whole tree are not significantly different than expected under BM (makes sense, since we only identified localized shifts in the tree earlier). Accordingly, the 'spread' metric is very close to 1 (the expected value under BM).
  527. However, **the phenotypic regression is highly significant** (p.random = 1). The "dev" metric quantifies the deviation of the phenotypic mean of our trait from the root value, expressed in standard deviations of the distribution. Here, dev = 1.85, meaning that the value of EQ is almost twice as high as the standard deviation expected from BM. *There is a strong trend in our phenotype!* Since the slope is positive, this means we have a significant increase in EQ somewhere in the tree.
  529. - Let's test the same nodes we tested earlier for shifts
  530. ```{r}
  531. STclade1<-search.trend(RR=RREQ1, y=EQnoOB, node=c(250,246,129), nsim=1000)
  532. STclade2<-search.trend(RR=RREQ1, y=EQnoOB, node=c(219,130,216), nsim=1000)
  533. STclade3<-search.trend(RR=RREQ1, y=EQnoOB, node=c(131,209,179), nsim=1000)
  534. STclade4<-search.trend(RR=RREQ1, y=EQnoOB, node=c(143,132,144), nsim=1000)
  535. STclade5<-search.trend(RR=RREQ1, y=EQnoOB, node=c(152,154,158), nsim=1000)
  537. # Plot trends in EQ value on the tree
  538. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  539. labelnodes(text=1:RREQ1$tree$Nnode,node=1:RREQ1$tree$Nnode+Ntip(RREQ1$tree),
  540. interactive=FALSE,circle.exp=0.4,cex=0.5)
  541. nodelabels(text="EQ+", node=c(129, 130, 131, 132, 143, 144, 152), adj=c(1.1, 1), frame="none", col="red2")
  543. # To get the values of EQ at nodes with significant shifts
  544. values<-as.data.frame(cbind(as.numeric(rownames(RREQ1$aces)), expm1(RREQ1$aces)))
  545. Ivalues<-values %>% filter (V1 %in% c(129:132, 143, 144, 152))
  546. rownames(Ivalues)<-c("Neotherapsida", "Theriodontia", "Eutheriodontia",
  547. "Cynodontia", "Eucynodontia", "Probainognathia",
  548. "Mammaliamorpha")
  549. colnames(Ivalues)<-c("Node number", "EQ")
  550. Ivalues
  551. ```
  552. For absolute values of EQ, we have **two trends of significant increases**:
  553. - One *early* in therapsid evolution: Neotherapsida, Theriodontia, Eutheriodontia, and Cynodontia
  554. - One *closer to the mammalian lineage*: Eucynodontia, Probainognathia, and Mammaliamorpha
  556. - However, we still need to test for the significance of these trends using `overfitRR`. As for shifts, we use a threshold of p = 0.75 for significance.
  557. ```{r}
  558. # Generate alternative topologies
  559. resampleEQ1trend<-resampleTree(RREQ1$tree, nsim=100)
  561. # Recompile regressions for the new trees
  562. overfitEQ1trend<-overfitRR(RREQ1, y=EQnoOB, phylo.list=resampleEQ1trend, clus=cc)
  564. # Test for trend significance with new topologies
  565. overtrend1<-overfitST(RR=RREQ1, y=EQnoOB, oveRR=overfitEQ1trend, node=c(129,143), clus=cc)
  566. overtrend2<-overfitST(RR=RREQ1, y=EQnoOB, oveRR=overfitEQ1trend, node=c(130,144), clus=cc)
  567. overtrend3<-overfitST(RR=RREQ1, y=EQnoOB, oveRR=overfitEQ1trend, node=c(131,152), clus=cc)
  568. overtrend4<-overfitST(RR=RREQ1, y=EQnoOB, oveRR=overfitEQ1trend, node=c(132), clus=cc)
  570. # Display results
  571. ## For the whole tree
  572. overtrend1$trend.results$tree$phenotype
  573. ## For specific nodes
  574. overtrend1$trend.results$node$phenotype
  575. overtrend2$trend.results$node$phenotype
  576. overtrend3$trend.results$node$phenotype
  577. overtrend4$trend.results$node$phenotype
  578. ```
  579. The trend of EQ increase is **highly significant** (p > 0.95) for the whole tree and for all tested nodes.
  581. - Let's do a nicer plot of the phenotypic mean on the tree so we can visualize the trend better
  582. ```{r}
  583. # Extract ancestral states for the phenotypic mean (EQ)
  584. phenoplot<-ST$trend.data$phenotypeVStime
  585. phenonames<-rownames(phenoplot); phenonameshort<-phenonames[which(nchar(rownames(phenoplot))==3)]
  586. phenoanc<-phenoplot[which(nchar(rownames(phenoplot))==3),1]; names(phenoanc)<-phenonameshort
  588. # Plot them on the calibrated tree
  589. phenoMapCal<-contMap(RREQ1$tree, EQnoOB, method="user", anc.states=phenoanc, plot=FALSE)
  590. ## Previous palette
  591. plot(setMap(phenoMapCal,colors=colorRampPalette(pal)(10)),fsize=0.5,lwd=4,cex=c(0.5,0.3))
  592. ## Rainbow palette
  593. plot(setMap(phenoMapCal,colors=rev(brewer.pal(10, "Spectral"))),fsize=0.5,lwd=4,cex=c(0.5,0.3))
  595. # Plot them on the tree with all branch lengths equal (easier to visualize)
  596. phenoMapEqual<-contMap(nodetree, EQnoOB, method="user", anc.states=phenoanc, plot=FALSE)
  597. ## Previous palette
  598. plot(setMap(phenoMapEqual,colors=colorRampPalette(pal)(10)),fsize=0.5,lwd=4,cex=c(0.5,0.3))
  599. ## Rainbow palette
  600. plot(setMap(phenoMapEqual,colors=rev(brewer.pal(10, "Spectral"))),fsize=0.5,lwd=4,cex=c(0.5,0.3))
  601. ```
  603. - Visualize the significant phenotypic shifts on the time-calibrated tree with added geological time scale, using `strap` (Bell & Lloyd, 2015)
  604. ```{r}
  605. # Simple version
  606. geoscalePhylo(ttree, units=c("Period", "Epoch"), boxes="Epoch", cex.tip=0.4,
  607. cex.ts=0.8, tick.scale="Age", x.lim=c(270, 55))
  608. nodelabels(text="EQ+", node=c(129, 130, 131, 132, 143, 144, 152),
  609. adj=c(1.1, 1), frame="none", col="red2")
  611. # Version with stratigraphic ranges of individual taxa in the tree
  612. ages<-maindata[,5:6]; rownames(ages)<-maindata$Taxon
  613. geoscalePhylo(ttree, ages, units=c("Period", "Epoch"), boxes="Epoch", cex.tip=0.4,
  614. cex.ts=0.8, tick.scale="Age", x.lim=c(270, 55))
  615. nodelabels(text="EQ+", node=c(129, 130, 131, 132, 143, 144, 152),
  616. adj=c(1.1, 1), frame="none", col="red2")
  617. ```
  619. - Same with significant shifts in evolutionary rate
  620. ```{r}
  621. # Simple version
  622. geoscalePhylo(ttree, units=c("Period", "Epoch"), boxes="Epoch", cex.tip=0.4,
  623. cex.ts=0.8, tick.scale="Age", x.lim=c(270, 55))
  624. nodelabels(text="-", node=c(144:152, 155), cex=1.5,
  625. adj=c(1.1, 1), frame="none", col="steelblue2")
  627. # Version with stratigraphic ranges of individual taxa in the tree
  628. ages<-maindata[,5:6]; rownames(ages)<-maindata$Taxon
  629. geoscalePhylo(ttree, ages, units=c("Period", "Epoch"), boxes="Epoch", cex.tip=0.4,
  630. cex.ts=0.8, tick.scale="Age", x.lim=c(270, 55))
  631. nodelabels(text="-", node=c(144:152, 155), cex=1.5,
  632. adj=c(1.1, 1), frame="none", col="steelblue2")
  633. ```
  635. ### For EQ with OB
  637. - Prepare the data and tree
  638. ```{r}
  639. # Dataset
  640. EQ2<-logmain[,5]; names(EQ2)<-rownames(logmain)
  642. # Proportion of clusters
  643. cc <- 2/parallel::detectCores()
  645. # Compile RRphylo
  646. RREQ2<-RRphylo(tree=ttree, y=EQ2, clus=cc)
  648. # To visualize node numbers in the tree
  649. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  650. labelnodes(text=1:RREQ2$tree$Nnode,node=1:RREQ2$tree$Nnode+Ntip(RREQ2$tree),
  651. interactive=FALSE,circle.exp=0.4,cex=0.5)
  652. ```
  654. #### Let's test for significant shifts in evolutionary rates in the tree, using `search.shift`
  656. - In the whole tree (as a first step, to see if any main trends can be identified)
  657. ```{r}
  658. shiftsWholeR<-search.shift(RREQ2, status.type= "clade")
  659. shiftsWholeR$all.clades
  660. as.numeric(rownames(shiftsWholeR$all.clades))-length(ttree$tip.label) # actual node numbers as visualized in plotTree
  662. # Significant shifts
  663. posshift<-as.numeric(names(which(shiftsWholeR$all.clades[,2]>0.975)))
  664. posshift # Significant positive shifts
  665. negshift<-as.numeric(names(which(shiftsWholeR$all.clades[,2]<0.025)))
  666. negshift # Significant negative shifts
  668. # Plot them on the tree
  669. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  670. labelnodes(text=1:RREQ2$tree$Nnode,node=1:RREQ2$tree$Nnode+Ntip(RREQ2$tree),
  671. interactive=FALSE,circle.exp=0.4,cex=0.5)
  672. nodelabels(text="+", node=posshift, adj=c(1.5, 1), frame="none", col="red2")
  673. nodelabels(text="-", node=negshift, adj=c(1.5, 1), frame="none", col="steelblue2")
  674. ```
  676. **Two major series of shifts, both in Eucynodontia**: higher rates in Cynognathia (suggesting a selective pressure on EQ in that clade), and lower rates in Probainognathia (suggesting relaxed constraints, potentially associated with ancestral stabilizing selection).
  678. Results are very similar to those obtained with EQ with no OB.
  680. - Let's make a phylomorphospace for EQ with OB and body mass, with colors for these two clades:
  681. ```{r}
  682. # Data
  683. EQOBdat<-logmain[,c(3,5)]
  685. # Basic phylomorphospace, if you need it – remove the 'label="off"' part if you want taxa names, but it gets very crowded!
  686. par(mar = c(5.1, 4.1, 4.1, 2.1))
  687. phylomorphospace(ttree,EQOBdat,bty="l",label="off",xlab="log (Body mass) (g)",
  688. ylab="log (EQ) (with OB)",node.size=c(0,1))
  689. title(main="Phylomorphospace plot",font.main=3)
  691. # Basic tree with labelled clades
  692. plotTree(ttree,fsize=0.5,ftype="i")
  693. nodelabels(frame="circ",bg="white",cex=0.3)
  694. cladelabels(ttree,c("Probainognathia","Cynognathia"),c(144,179))
  696. # Define colors for the tree
  697. painted<-paintSubTree(ttree, 144, "Cynognathia")
  698. painted<-paintSubTree(painted, 179, "Probainognathia")
  700. # Plot it on the tree
  701. plot(painted, fsize=0.5, ftype="i")
  702. cladelabels(ttree, c("Probainognathia","Cynognathia"), c(144,179), offset=0.5)
  704. # Plot it on the phylomorphospace
  705. par(mar = c(5.1, 4.1, 4.1, 2.1))
  706. phylomorphospace(painted,EQOBdat,bty="l",label="off",xlab="log (Body mass) (g)",
  707. ylab="log (EQ) (with OB)",node.size=c(0,1.2),node.by.map=TRUE)
  708. title(main="Phylomorphospace plot",font.main=3)
  709. legend(x="topleft",legend=c("Cynognathia","Probainognathia"),
  710. pch=21,pt.cex=1.5,pt.bg=palette()[2:4],bty="n")
  712. # Plot it on the phenogram
  713. EQOB<-logmain[,5]; names(EQOB)<-rownames(logmain)
  714. par(mar = c(5.1, 4.1, 4.1, 2.1))
  715. phenogram(painted, EQOB, ,ftype="off", spread.labels=TRUE,
  716. spread.cost=c(1,0),fsize=0.7,
  717. xlab="Time (Ma)",ylab="log (EQ) (with OB)",cex.axis=0.8,
  718. axes=FALSE,las=1, lwd=1)
  719. axis(1,at=round(seq(0,max(nodeHeights(painted)),length.out=6),1),
  720. label=round(seq(max(nodeHeights(painted)),0,length.out=6),1),
  721. cex.axis=0.8)
  722. axis(2,las=1,cex.axis=0.8)
  723. grid()
  724. title(paste("Phenogram for EQ (with OB)"))
  725. legend(x="bottomright",legend=c("Probainognathia","Cynognathia"),
  726. pch=21,pt.cex=1.5,pt.bg=palette()[2:4],bty="n")
  727. ```
  729. As expected, Cynognathia show very high variation in a short amount of time, while the earliest nodes within Probainognathia seem to show limited variation over an extended time frame compared to more inclusive clades in the tree.
  731. Let's go a bit further... Same test, but for specific clades in the tree!
  733. We will be testing for shifts in the following clades:
  734. Biarmosuchia, Dinocephalia, Neotherapsida, Dicynodontia, Theriodontia, Gorgonopsia, Eutheriodontia, Therocephalia, Cynodontia, Eucynodontia, Cynognathia, Probainognathia, Mammaliamorpha, Mammaliaformes, Mammalia.
  736. - Specific shifts in those clades
  737. ```{r, warning=FALSE}
  738. # Node numbers
  739. nodes<-c(124,120,3,93,4,90,5,83,6,17,53,18,26,28,32)
  740. cornodes<-nodes+length(ttree$tip.label)
  742. # Look for significant shifts
  743. shiftsR<-search.shift(RREQ2, status.type= "clade", node=cornodes)
  744. shiftsR$single.clades
  745. which(shiftsR$single.clades[,2]>0.975); which(shiftsR$single.clades[,2]<0.025)
  747. posshift<-as.numeric(names(which(shiftsR$single.clades[,2]>0.975)))
  748. posshift # Significant positive shifts
  749. negshift<-as.numeric(names(which(shiftsR$single.clades[,2]<0.025)))
  750. negshift # Significant negative shifts
  752. # Plot them on the tree
  753. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  754. labelnodes(text=1:RREQ2$tree$Nnode,node=1:RREQ2$tree$Nnode+Ntip(RREQ2$tree),
  755. interactive=FALSE,circle.exp=0.4,cex=0.5)
  756. #nodelabels(text="+", node=posshift, adj=c(1.5, 1), frame="none", col="red2")
  757. nodelabels(text="-", node=negshift, adj=c(1.5, 1), frame="none", col="steelblue2")
  758. ```
  760. Results are congruent with those of the initial 'auto-recognize' search, but only clades that had already been identified as presenting a significant negative shift **in Probainognathia** are recovered as such.
  762. - We can plot the shifts on the tree
  763. ```{r}
  764. RRplotR<-plotRR(RREQ2, y=EQOB, multivariate = "rates")
  766. ## Recompile individual shifts for the plot
  767. posshift<-as.numeric(names(which(shiftsWholeR$all.clades[,2]>0.975)))
  768. negshift<-as.numeric(names(which(shiftsWholeR$all.clades[,2]<0.025)))
  769. EQ2shift<-search.shift(RREQ2, status.type= "clade", node=c(posshift, negshift))
  770. EQ2shift$single.clades<-as.data.frame(EQ2shift$single.clades)
  771. which(EQ2shift$single.clades[,2]>0.975); which(EQ2shift$single.clades[,2]<0.025)
  773. ## Values of EQ
  774. RRplotR$plotRRphen(variable=1, tree.args=list(no.margin=TRUE),
  775. colorbar.args=list(x=170, y=90, title.pos="bottom"))
  777. ## Rates and rate shifts on the same plot
  778. RRplotR$plotRRrates(variable=1, tree.args=list(no.margin=TRUE),
  779. colorbar.args=list(x=170, y=90, title.pos="bottom"))
  780. addShift(EQ2shift, symbols.args=list(lwd=2,fg=c(pos="red2",neg="steelblue2"),
  781. bg=scales::alpha(c(rep("red2",length(posshift)),rep("steelblue2",length(negshift))),0.3)))
  783. ## Only rate shifts
  784. plotShiftR<-plotShift(RREQ2, EQ2shift)
  785. plotShiftR$plotClades(tree.args=list(no.margin=TRUE),
  786. symbols.args=list(lwd=2,fg=c(pos="red2",neg="steelblue2"),
  787. bg=scales::alpha(c(rep("red2",length(posshift)),rep("steelblue2",length(negshift))),0.3)))
  788. ```
  790. - Finally, let's test for the significance of those shifts with jacknife/multiple permutation tests, using `overfitRR`. We use a threshold of 0.75 (3/4 of simulated trees recovered with the same shift) to assess a shift as significant.
  791. ```{r}
  792. # Generate alternative topologies
  793. resampleEQ2<-resampleTree(RREQ2$tree, node=c(posshift, negshift), nsim=1000)
  795. # Recompile regressions for the new trees
  796. overfitEQ2<-overfitRR(RREQ2, y=EQOB, phylo.list=resampleEQ2, clus=cc)
  798. # Test for shift significance with new topologies
  799. overshift1<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(181,144))
  800. overshift2<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(184,145))
  801. overshift3<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(185,146))
  802. overshift4<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(186,147))
  803. overshift5<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(156,148))
  804. overshift6<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(155,149))
  805. overshift7<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(154,150))
  806. overshift8<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(153,151))
  807. overshift9<-overfitSS(RR=RREQ2, oveRR=overfitEQ2, node=c(152,183))
  809. # Display the results
  810. overshift1$shift.results$clade$single.clades
  811. overshift2$shift.results$clade$single.clades
  812. overshift3$shift.results$clade$single.clades
  813. overshift4$shift.results$clade$single.clades
  814. overshift5$shift.results$clade$single.clades
  815. overshift6$shift.results$clade$single.clades
  816. overshift7$shift.results$clade$single.clades
  817. overshift8$shift.results$clade$single.clades
  818. overshift9$shift.results$clade$single.clades
  819. ```
  820. Nodes 18 through 30 (i.e. all clades within Probainognathia) all have significant negative shifts, whereas nodes within Cynodontia (positive shifts) do not.
  821. => **Only the rate shifts in Probainognathia and right outside Eucynodontia (negative, i.e. stabilizing selection of EQ) are recovered as significant.**
  823. Now that we have identified shifts in evolutionary rates, what about phenotypic trends (i.e. changes in the actual values of EQ through time)?
  825. #### Testing for significant phenotypic trends for EQ, using `search.trend`
  827. See the vignette and help section of the function, as well as Castiglione et al. (2019). We are testing for significant deviations from a Brownian Motion (BM) model in the value of our trait of interest (EQ).
  829. - For the whole tree
  830. ```{r}
  831. ST<-search.trend(RR=RREQ2, y=EQOB, nsim=1000)
  833. # Display the results
  834. rateres<-cbind(ST$rate.regression[1,1:4], NA); colnames(rateres)[5]<-"dev"
  835. phenres<-c(ST$phenotypic.regression[1,1:3], NA, ST$phenotypic.regression[1,4])
  836. names(phenres)[c(4,5)]<-c("spread", "dev")
  837. STres<-rbind(rateres, phenres)
  838. rownames(STres)<-c("rescaled absolute rate regression", "phenotypic regression")
  839. STres
  840. ```
  841. The absolute rate regression is not significant (p.random < 0.95), meaning the rates for the whole tree are not significantly different than expected under BM (makes sense, since we only identified localized shifts in the tree earlier). Accordingly, the 'spread' metric is very close to 1 (the expected value under BM).
  843. However, **the phenotypic regression is highly significant** (p.random = 1). The "dev" metric quantifies the deviation of the phenotypic mean of our trait from the root value, expressed in standard deviations of the distribution. Here, dev = 2, meaning that the value of EQ is almost twice as high as the standard deviation expected from BM. *There is a strong trend in our phenotype!* Since the slope is positive, this means we have a significant increase in EQ somewhere in the tree.
  845. - Let's test the same nodes we tested earlier for shifts
  846. ```{r}
  847. STclade1<-search.trend(RR=RREQ2, y=EQOB, node=c(250,246,129), nsim=1000)
  848. STclade2<-search.trend(RR=RREQ2, y=EQOB, node=c(219,130,216), nsim=1000)
  849. STclade3<-search.trend(RR=RREQ2, y=EQOB, node=c(131,209,179), nsim=1000)
  850. STclade4<-search.trend(RR=RREQ2, y=EQOB, node=c(143,132,144), nsim=1000)
  851. STclade5<-search.trend(RR=RREQ2, y=EQOB, node=c(152,154,158), nsim=1000)
  853. # Plot trends in EQ value on the tree
  854. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  855. labelnodes(text=1:RREQ2$tree$Nnode,node=1:RREQ2$tree$Nnode+Ntip(RREQ2$tree),
  856. interactive=FALSE,circle.exp=0.4,cex=0.5)
  857. nodelabels(text="EQ+", node=c(129, 130, 131, 132, 143, 144, 152, 154),
  858. adj=c(1.1, 1), frame="none", col="red2")
  860. # To get the values of EQ at nodes with significant shifts
  861. values2<-as.data.frame(cbind(as.numeric(rownames(RREQ2$aces)), expm1(RREQ2$aces)))
  862. Ivalues2<-values2 %>% filter (V1 %in% c(129:132, 143, 144, 152, 154))
  863. rownames(Ivalues2)<-c("Neotherapsida", "Theriodontia", "Eutheriodontia",
  864. "Cynodontia", "Eucynodontia", "Probainognathia",
  865. "Mammaliamorpha", "Mammaliaformes")
  866. colnames(Ivalues2)<-c("Node number", "EQ")
  867. Ivalues2
  868. ```
  870. Significant shifts in rates are the same as for our previous analyses: a positive shift for Cynognathia and a negative shift for Probainognathia.
  872. For absolute values of EQ, we have **two trends of significant increases**:
  873. - One *early* in therapsid evolution: Neotherapsida, Theriodontia, Eutheriodontia, and Cynodontia
  874. - One *closer to the mammalian lineage*: Eucynodontia, Probainognathia, Mammaliamorpha, and Mammaliaformes
  876. - However, we still need to test for the significance of these trends using `overfitRR`. As for shifts, we use a threshold of p = 0.75 for significance.
  877. ```{r}
  878. # Generate alternative topologies
  879. resampleEQ2trend<-resampleTree(RREQ2$tree, nsim=100)
  881. # Recompile regressions for the new trees
  882. overfitEQ2trend<-overfitRR(RREQ2, y=EQOB, phylo.list=resampleEQ2trend, clus=cc)
  884. # Test for trend significance with new topologies
  885. overtrend1<-overfitST(RR=RREQ2, y=EQOB, oveRR=overfitEQ2trend, node=c(129,154), clus=cc)
  886. overtrend2<-overfitST(RR=RREQ2, y=EQOB, oveRR=overfitEQ2trend, node=c(130,144), clus=cc)
  887. overtrend3<-overfitST(RR=RREQ2, y=EQOB, oveRR=overfitEQ2trend, node=c(131,152), clus=cc)
  888. overtrend4<-overfitST(RR=RREQ2, y=EQOB, oveRR=overfitEQ2trend, node=c(132,143), clus=cc)
  890. # Display results
  891. ## For the whole tree
  892. overtrend1$trend.results$tree$phenotype
  893. ## For specific nodes
  894. overtrend1$trend.results$node$phenotype
  895. overtrend2$trend.results$node$phenotype
  896. overtrend3$trend.results$node$phenotype
  897. overtrend4$trend.results$node$phenotype
  898. ```
  899. The trend of EQ increase is **highly significant** (p > 0.95) for the whole tree and for all tested nodes except Mammaliaformes, which makes the results identical to those of analyses with no OB.
  901. - Let's do a nicer plot of the phenotypic mean on the tree so we can visualize the trend better
  902. ```{r}
  903. # Extract ancestral states for the phenotypic mean (EQ)
  904. phenoplot<-ST$trend.data$phenotypeVStime
  905. phenonames<-rownames(phenoplot); phenonameshort<-phenonames[which(nchar(rownames(phenoplot))==3)]
  906. phenoanc<-phenoplot[which(nchar(rownames(phenoplot))==3),1]; names(phenoanc)<-phenonameshort
  908. # Plot them on the calibrated tree
  909. phenoMapCal<-contMap(RREQ1$tree, EQOB, method="user", anc.states=phenoanc, plot=FALSE)
  910. ## Previous palette
  911. plot(setMap(phenoMapCal,colors=colorRampPalette(pal)(10)),fsize=0.5,lwd=4,cex=c(0.5,0.3))
  912. ## Rainbow palette
  913. plot(setMap(phenoMapCal,colors=rev(brewer.pal(10, "Spectral"))),fsize=0.5,lwd=4,cex=c(0.5,0.3))
  915. # Plot them on the tree with all branch lengths equal (easier to visualize)
  916. phenoMapEqual<-contMap(nodetree, EQOB, method="user", anc.states=phenoanc, plot=FALSE)
  917. ## Previous palette
  918. plot(setMap(phenoMapEqual,colors=colorRampPalette(pal)(10)),fsize=0.5,lwd=4,cex=c(0.5,0.3))
  919. ## Rainbow palette
  920. plot(setMap(phenoMapEqual,colors=rev(brewer.pal(10, "Spectral"))),fsize=0.5,lwd=4,cex=c(0.5,0.3))
  921. ```
  923. - Same with evolutionary rates
  924. ```{r}
  925. # Extract ancestral states for evolutionary rates
  926. rateplot<-ST$trend.data$rateVStime
  927. ratenames<-rownames(rateplot)
  928. ratenameshort<-ratenames[which(nchar(rownames(rateplot))==3)]
  929. rateanc<-rateplot[which(nchar(rownames(rateplot))==3),1]
  930. names(rateanc)<-ratenameshort
  932. # Plot them on the calibrated tree
  933. rateMapCal<-contMap(RREQ1$tree, EQOB, method="user", anc.states=rateanc, plot=FALSE)
  934. ## Previous palette
  935. plot(setMap(rateMapCal,colors=colorRampPalette(pal)(10)),
  936. fsize=0.5,lwd=4,cex=c(0.5,0.3))
  937. ## Rainbow palette
  938. plot(setMap(rateMapCal,colors=rev(brewer.pal(10,"Spectral"))),
  939. fsize=0.5,lwd=4,cex=c(0.5,0.3))
  941. # Plot them on the tree with all branch lengths equal (easier to visualize)
  942. rateMapEqual<-contMap(nodetree, EQOB, method="user", anc.states=rateanc, plot=FALSE)
  943. ## Previous palette
  944. plot(setMap(rateMapEqual,colors=colorRampPalette(pal)(10)),
  945. fsize=0.5,lwd=4,cex=c(0.5,0.3))
  946. ## Rainbow palette
  947. plot(setMap(rateMapEqual,colors=rev(brewer.pal(10,"Spectral"))),
  948. fsize=0.5,lwd=4,cex=c(0.5,0.3))
  949. ```
  952. - Visualize the phenotypic shifts on the time-calibrated tree with added geological time scale, using `strap` (Bell & Lloyd, 2015)
  953. ```{r}
  954. # Simple version
  955. geoscalePhylo(ttree, units=c("Period", "Epoch"), boxes="Epoch", cex.tip=0.4,
  956. cex.ts=0.8, tick.scale="Age", x.lim=c(270, 55))
  957. nodelabels(text="EQ+", node=c(129, 130, 131, 132, 143, 144, 152),
  958. adj=c(1.1, 1), frame="none", col="red2")
  960. # Version with stratigraphic ranges of individual taxa in the tree
  961. ages<-maindata[,5:6]; rownames(ages)<-maindata$Taxon
  962. geoscalePhylo(ttree, ages, units=c("Period", "Epoch"), boxes="Epoch", cex.tip=0.4,
  963. cex.ts=0.8, tick.scale="Age", x.lim=c(270, 55))
  964. nodelabels(text="EQ+", node=c(129, 130, 131, 132, 143, 144, 152),
  965. adj=c(1.1, 1), frame="none", col="red2")
  966. ```
  968. - Same with significant shifts in evolutionary rate
  969. ```{r}
  970. # Simple version
  971. geoscalePhylo(ttree, units=c("Period", "Epoch"), boxes="Epoch", cex.tip=0.4,
  972. cex.ts=0.8, tick.scale="Age", x.lim=c(270, 55))
  973. nodelabels(text="-", node=c(142, 144:146, 149:152, 154, 155), cex=1.5,
  974. adj=c(1.1, 1), frame="none", col="steelblue2")
  976. # Version with stratigraphic ranges of individual taxa in the tree
  977. ages<-maindata[,5:6]; rownames(ages)<-maindata$Taxon
  978. geoscalePhylo(ttree, ages, units=c("Period", "Epoch"), boxes="Epoch", cex.tip=0.4,
  979. cex.ts=0.8, tick.scale="Age", x.lim=c(270, 55))
  980. nodelabels(text="-", node=c(142, 144:146, 149:152, 154, 155), cex=1.5,
  981. adj=c(1.1, 1), frame="none", col="steelblue2")
  982. ```
  984. ## New analyses of evolutionary shifts and regimes in an OU framework
  986. ### Using `phylolm` (Ho & Ané, 2014)
  987. ```{r, warning=FALSE}
  988. ## For EQ with no OB
  989. result<-OUshifts(EQ1, ttree, method="mbic", nmax=ttree$Nnode)
  990. result$mean; result$pshift; result$shift
  991. par(mar=c(0,0,0,0)); plot.OUshifts(result, cex=0.5, show.data=FALSE)
  993. ## For EQ with OB
  994. result2<-OUshifts(EQ2, ttree, method="mbic", nmax=ttree$Nnode)
  995. result2$mean; result2$pshift; result2$shift
  996. par(mar=c(0,0,0,0)); plot.OUshifts(result2, cex=0.5, show.data=FALSE)
  998. ## Tree with edge labels to visualize numbers
  999. plotTree(nodetree, fsize=0.5, lwd=1, ftype="i")
  1000. edgelabels(frame="circle", cex=0.5, bg="white")
  1001. ```
  1002. For both EQ (with and without OB), the results are identical: four shifts recovered, only one of which is at an internal branch: branch 29, which leads to node 30: Mammalia + *Cifelliodon* + *Hadrocodium*. It is a positive shift.
  1004. The other three shifts are recovered at the following terminal taxa (one therocephalian and two non-eucynodont cynodonts):
  1005. - *Galesaurus planiceps* 1 (positive shift);
  1006. - *Thrinaxodon liorhinus* 3 (negative shift);
  1007. - *Tetracynodon darti* 1 (negative shift).
  1009. This seems to be congruent with our previous results obtained with `RRphylo`: the only significant shifts occur near the origin of mammals, with more inclusive clades not affected by the PTME and not showing any conspicuous patterns. We do not, however, recover any shifts near the base of the tree.
  1011. ### Using `slouch` (Kopperud et al., 2024)
  1013. - Preparing the nodes in the tree
  1014. ```{r}
  1015. # Defining regimes in the tree (different clades of interest – same as those used with `RRphylo`)
  1016. nodenames<-c(rep("Therapsida", 2), "Neotherapsida", "Theriodontia",
  1017. "Eutheriodontia", rep("Cynodontia", 11), "Eucynodontia",
  1018. rep("Probainognathia", 8), rep("Mammaliamorpha", 2),
  1019. rep("Mammaliaformes", 4), rep("Mammalia", 10),
  1020. rep("Mammaliamorpha", 5), rep("Probainognathia", 6),
  1021. rep("Cynognathia", 22), rep("Cynodontia", 8),
  1022. rep("Therocephalia",7), rep("Gorgonopsia", 3),
  1023. rep("Dicynodontia", 27), rep("Dinocephalia", 4),
  1024. rep("Biarmosuchia", 2))
  1026. ttree$node.label<-nodenames # Add them as node names in the tree
  1028. # Map them on the tree
  1029. tipnames<-c(rep("Biarmosuchia", 3), rep("Dinocephalia", 2),
  1030. rep("Dicynodontia", 18), "Gorgonopsia", rep("Therocephalia", 8),
  1031. rep("Cynodontia", 13), rep("Cynognathia", 13),
  1032. rep("Probainognathia", 8), rep("Mammaliamorpha", 4),
  1033. rep("Mammaliaformes", 4), rep("Mammalia", 11),
  1034. rep("Dicynodontia", 4), rep("Cynognathia", 7), rep("Dinocephalia", 3),
  1035. rep("Dicynodontia", 6), rep("Gorgonopsia", 3), rep("Cynodontia", 6),
  1036. rep("Cynognathia", 3), rep("Probainognathia", 6),
  1037. rep("Mammaliamorpha", 3))
  1039. regimes<-c(as.factor(tipnames), as.factor(nodenames))
  1040. levels(regimes)<-c("Therapsida", "Biarmosuchia", "Dinocephalia", "Neotherapsida", "Dicynodontia", "Theriodontia", "Gorgonopsia", "Eutheriodontia", "Therocephalia", "Cynodontia", "Eucynodontia", "Cynognathia", "Probainognathia", "Mammaliamorpha", "Mammaliaformes", "Mammalia")
  1041. ```
  1043. - Color palette for different regimes
  1044. ```{r}
  1045. # Custom palette with 16 colors
  1046. rainbowpal<-c(brewer.pal(9,"Set1"), brewer.pal(7,"Set2"))
  1047. show_col(rainbowpal) # Looks good
  1048. names(rainbowpal)<-levels(regimes)
  1050. edge_regimes<-factor(regimes[ttree$edge[,2]])
  1052. par(mar=c(1,1,1,1)+0.1)
  1053. plot(ttree,
  1054. edge.color = rainbowpal[edge_regimes],
  1055. edge.width = 3, cex = 0.6) # Groups are correctly separated by color
  1056. ```
  1058. - Build the model with those regimes
  1059. ```{r}
  1060. slouchmod <- slouch.fit(phy = ttree,
  1061. species = logmain$Taxon,
  1062. response = logmain$EQ_noOB,
  1063. fixed.fact = as.factor(tipnames))
  1064. summary(slouchmod)
  1066. # BM models for comparison
  1067. ## With no regimes (simplest)
  1068. BMmod <- brown.fit(phy = ttree,
  1069. species = logmain$Taxon,
  1070. response = logmain$EQ_noOB)
  1071. summary(BMmod)
  1073. ## With regimes (trend)
  1074. BMmodtrend <- brown.fit(phy = ttree,
  1075. species = logmain$Taxon,
  1076. response = logmain$EQ_noOB,
  1077. fixed.fact = as.factor(tipnames))
  1078. summary(BMmodtrend)
  1080. slouchmod$modfit$AICc; BMmod$modfit$AICc; BMmodtrend$modfit$AICc
  1081. # OU is a better fit than BM here
  1083. slouchmod$beta_primary$coefficients # optima for each node of interest
  1084. ```
  1086. ## References
  1087. - Bell, M.A., Lloyd, G.T., 2015. strap: an R package for plotting phylogenies against stratigraphy and assessing their stratigraphic congruence. *Palaeontology* 58, 379–389. <https://doi.org/10.1111/pala.12142>
  1088. - Castiglione, S., Serio, C., Mondanaro, A., di Febbraro, M., Profico, A., Girardi, G., Raia, P., 2019. Simultaneous detection of macroevolutionary patterns in phenotypic means and rate of change with and within phylogenetic trees including extinct species. *PLOS ONE* 14, e0210101. <https://doi.org/10.1371/journal.pone.0210101>
  1089. - Castiglione, S., Tesone, G., Piccolo, M., Melchionna, M., Mondanaro, A., Serio, C., di Febbraro, M., Raia, P., 2018. A new method for testing evolutionary rate variation and shifts in phenotypic evolution. *Methods in Ecology and Evolution* 9, 974–983. <https://doi.org/10.1111/2041-210X.12954>
  1090. - Ho, L.S.T., Ané, C., 2014. A linear-time algorithm for Gaussian and non-Gaussian trait evolution models. *Systematic Biology* 63, 397–408. <https://doi.org/10.1093/sysbio/syu005>
  1091. - Goolsby, E.W., Bruggeman, J., Ané, C., 2017. Rphylopars: fast multivariate phylogenetic comparative methods for missing data and within-species variation. *Methods in Ecology and Evolution* 8, 22–27. <https://doi.org/10.1111/2041-210X.12612>
  1092. - Kopperud, B.T., Pienaar, J., Voje, K.L., Orzack, S.H., Hansen, T.F., Grabowski, M., 2024. *slouch: Stochastic Linear Ornstein-Uhlenbeck Comparative Hypotheses.* R package, available at: <https://cran.r-project.org/web/packages/slouch/index.html>
  1093. - Kratsch, C., McHardy, A.C., 2014. RidgeRace: ridge regression for continuous ancestral character estimation on phylogenetic trees. *Bioinformatics* 30, i527–i533. <https://doi.org/10.1093/bioinformatics/btu477>
  1094. - Pennell, M.W., Eastman, J.M., Slater, G.J., Brown, J.W., Uyeda, J.C., FitzJohn, R.G., Alfaro, M.E., Harmon, L.J., 2014. geiger v2.0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. *Bioinformatics* 30, 2216–2218. <https://doi.org/10.1093/bioinformatics/btu181>
  1095. - Revell, L.J., 2024. phytools 2.0: an updated R ecosystem for phylogenetic comparative methods (and other things). *PeerJ* 12, e16505. <https://doi.org/10.7717/peerj.16505>

Script_Benoitetal_2026.Rmd at commit 521ae20, no license · at the source

Overview

Authors: Julien Benoit1, Lucas J. Legendre2, Ricardo Araújo3, Vincent Fernandez4, Adam Midzuk1, Claire Browning5,6, Fernando Abdala1,7, Jennifer Botha1,8, Kenneth D. Angielczyk1,9
ORCID iDs: Lucas J. Legendre
  1. Evolutionary Studies Institute, School of Geosciences, University of the Witwatersrand,Private Bag 3, Johannesburg, 2050 WITS South Africa
  2. Department of Earth and Planetary Sciences, The University of Texas at Austin,2305 Speedway Stop C1160, Austin, TX 78712 USA
  3. Centro de Recursos Naturais e Ambiente (CERENA), Instituto Superior Técnico, Universidade de Lisboa,Lisboa, Portugal
  4. European Synchrotron Radiation Facility,71 rue des Martyrs, Grenoble, France
  5. Iziko Museums of South Africa,Cape Town, 8000 South Africa
  6. Department of Geological Sciences, University Avenue South, University of Cape Town Rondebosch,13 University Avenue South, Cape Town, South Africa
  7. Unidad Ejecutora Lillo (CONICET-Fundación Miguel Lillo),Miguel Lillo 251, Tucumán San Miguel de Tucumán, Argentina
  8. GENUS: DSTI-NRF Centre of Excellence in Palaeosciences, University of the Witwatersrand,Johannesburg, South Africa
  9. Negaunee Integrative Research Center, Field Museum of Natural History,1400 South DuSable Lake Shore Drive, Chicago, IL 60605 USA
Institutions: University of the Witwatersrand (South Africa); The University of Texas at Austin (United States); University of Lisbon (Portugal); European Synchrotron Radiation Facility (France); Iziko Museums of South Africa (South Africa); University of Cape Town (South Africa); Fundación Miguel Lillo (Argentina); Field Museum of Natural History (United States)
Journal: Scientific reports, volume 16, issue 1, article 22211
Dates: received 27 February 2026; accepted 11 May 2026; published online 15 May 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s41598-026-53133-y · PMID 42141026 · PMCID PMC13369837 · OpenAlex W7161248481
Open access: gold, a free copy (OpenAlex)
Status: code verified
Methods: Statistics
Keywords: Biological crisis, Encephalization, Brain, Therapsida, Permian-Triassic extinction, Ecology, Evolution, Neuroscience
MeSH: Biological Evolution*, Brain*, Extinction, Biological*, Fossils*, Mammals*, Animals, Phylogeny (* major topic)
Topic: Protist diversity and phylogeny (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: DSTI-NRF African Origins Platform (AOP210218587003, AOP240418214774, AOP240326210961); Consejo Nacional de Investigaciones Científicas y Técnicas; GENUS: DSTI-NRF Centre of Excellence in Palaeosciences
Citations: not cited yet (Europe PMC); 119 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repository

Its files are read in the Code ↔ Paper reader above, with 6 matches between paragraphs and lines of code.

LucasLegendre/synapsid_EQ

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 521ae2067152c788ee0540a74c053e26702a176c, 22 April 2026
Languages: R (2)
Size: 7 files, 2 scripts
Software Heritage: not archived
Found in: the text, “Phylogenetic comparative analyses”
Holds: README, CITATION.cff, 2 notebooks
Not found: license file, environment file, tests, continuous integration, documentation
Tools: nlme (2 files), tidyverse (2 files)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
3 files
At the source: github.com/LucasLegendre/synapsid_EQ/

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 2 scripts, each with its path and the digest of its content;
  • 6 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Code and data availability statement

The paper has a code and data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • no repository, dataset or request procedure was recognized in it

Read it in the paper: doi.org/10.1038/s41598-026-53133-y.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 8 keywords, 7 MeSH terms, 3 funders, 113 references.

Cite

This paper

Benoit, J., Legendre, L. J., Araújo, R., Fernandez, V., Midzuk, A., Browning, C., Abdala, F., Botha, J., & Angielczyk, K. D. (2026). Evolutionary stasis in synapsid encephalization during the end-permian mass extinction. Scientific reports, 16(1), 22211. https://doi.org/10.1038/s41598-026-53133-y

BibTeX

@article{benoit2026evolutionary,
author = {Benoit, Julien and Legendre, Lucas J. and Araújo, Ricardo and Fernandez, Vincent and Midzuk, Adam and Browning, Claire and Abdala, Fernando and Botha, Jennifer and Angielczyk, Kenneth D.},
title = {{Evolutionary stasis in synapsid encephalization during the end-permian mass extinction}},
journal = {Scientific reports},
year = {2026},
month = may,
volume = {16},
number = {1},
pages = {22211},
publisher = {Nature Publishing Group},
issn = {2045-2322},
doi = {10.1038/s41598-026-53133-y},
url = {https://doi.org/10.1038/s41598-026-53133-y},
pmid = {42141026},
pmcid = {PMC13369837}
}

RIS

TY - JOUR
AU - Benoit, Julien
AU - Legendre, Lucas J.
AU - Araújo, Ricardo
AU - Fernandez, Vincent
AU - Midzuk, Adam
AU - Browning, Claire
AU - Abdala, Fernando
AU - Botha, Jennifer
AU - Angielczyk, Kenneth D.
TI - Evolutionary stasis in synapsid encephalization during the end-permian mass extinction
T2 - Scientific reports
J2 - Sci Rep
PY - 2026
DA - 2026/05/15
VL - 16
IS - 1
SP - 22211
SN - 2045-2322
PB - Nature Publishing Group
DO - 10.1038/s41598-026-53133-y
UR - https://doi.org/10.1038/s41598-026-53133-y
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41598-026-53133-y",
"type": "article-journal",
"title": "Evolutionary stasis in synapsid encephalization during the end-permian mass extinction",
"container-title": "Scientific reports",
"author": [
{
"family": "Benoit",
"given": "Julien"
},
{
"family": "Legendre",
"given": "Lucas J."
},
{
"family": "Araújo",
"given": "Ricardo"
},
{
"family": "Fernandez",
"given": "Vincent"
},
{
"family": "Midzuk",
"given": "Adam"
},
{
"family": "Browning",
"given": "Claire"
},
{
"family": "Abdala",
"given": "Fernando"
},
{
"family": "Botha",
"given": "Jennifer"
},
{
"family": "Angielczyk",
"given": "Kenneth D."
}
],
"container-title-short": "Sci Rep",
"volume": "16",
"issue": "1",
"page": "22211",
"DOI": "10.1038/s41598-026-53133-y",
"PMID": "42141026",
"PMCID": "PMC13369837",
"ISSN": "2045-2322",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41598-026-53133-y",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
15
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1371/journal.pbio.3003771 [code]
Bipedalism and brain expansion explain human handedness.
Journal: PLoS biology
In common: 3 references
[2] doi:10.3390/biology15161373
Elevational Variation in Organ Size and Brain-Viscera Energetic Trade-Offs in the Minshan Toad (&lt;i&gt;Bufo minshanicus&lt;/i&gt;).
Journal: Biology
In common: 3 references
[3] doi:10.1038/s41586-026-10809-9 [code]
Exceptional brain and ecological diversity in the earliest snakes.
Journal: Nature
In common: 3 references
[4] doi:10.1093/nar/gkag544 [code]
Single-base resolution atlas reveals moderate conservation and regulatory diversity of m6A modifications across mammals.
Journal: Nucleic acids research
In common: nlme, tidyverse, 1 reference
[5] doi:10.3897/bdj.14.e191439 [code]
Assessing patterns of extinction risk amongst mammal species in Nigeria: A comparative analysis of human impact.
Journal: Biodiversity data journal
In common: tidyverse, 2 references
[6] doi:10.1038/s41514-026-00439-w [code]
Effects of a three-month exercise programme on cognition, mood and neurogenesis: the NeuroFit randomised controlled trial.
Journal: npj aging
In common: nlme, tidyverse
[7] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: nlme, tidyverse
[8] doi:10.1126/sciadv.aec9291 [code]
Computational mechanisms of perception in autism revealed using games inspired by rodent operant tasks.
Journal: Science advances
In common: nlme, tidyverse
[9] doi:10.1093/braincomms/fcag343 [code]
Long-term brain volume trajectories and lifestyle associations in cognitively normal adults: the BRAIN-STRIDE study.
Journal: Brain communications
In common: nlme, tidyverse
[10] doi:10.1162/imag.a.1337 [code]
Data quality biases normative models derived from fetal brain MRI.
Journal: Imaging neuroscience (Cambridge, Mass.)
In common: nlme, tidyverse

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.