Note that I have attached a .docx version of the review in case the tables and other elements are not displayed correctly here.

Overall, the authors have made a significant effort to allow the reproducibility of their analyses. The methods reported correspond the R script and models provided. The models are running properly (although I did not run them entirely, as the process is time-consuming), and the authors have uploaded .rds files containing the results of the models without fully executing them. Additionally, the readme file attached to the data and code accurately reports the dataset's structure and content.
However, I noted some issues that need to be addressed before the manuscript is published. More specifically, the raw data are unavailable; however, some outliers have been removed, and certain conversions have been applied to the current dataset. For full reproducibility, it would be beneficial to upload the raw data along with the procedures used to remove the outliers and the R code for the necessary transformation of tarsus length. Alternatively, it is essential to indicate where readers can find the raw data (e.g., SPI-Birds network with the URL, if possible), state which outliers have been removed (if applicable), and specify the conversions performed on this dataset within the Readme file.
Also, I’ve noticed some inconsistencies between the values reported in the manuscript and those I found by inspecting the data (see below for details). In addition, the figures are not directly reproducible (i.e., Figures 2 and 3) with the provided code (but see R code at the end of this review to programmatically reproduce Table 1 and the tables included in Figure 2).
Finally, I suggest adding a Licence file to the repository and citing the licence of the data and code within the Readme files (for help in choosing a licence, see https://choosealicense.com/).

Detailed comments:
- In the dataset ThompsonEA_UrbanVp_NstTarsus.csv, the number of days between hatching and measurements ranges from 13 to 17 days, rather than the 12 to 17 days reported in L209 and L361.
- The raw data and subsequent conversion of tarsus length referenced in L211 are not included in the dataset.
- The range of adult tarsus lengths is between 12.70 and 25.49 mm, which differs from the reported range of 11 to 28 mm in L318.
- There are 12 cities included in the analyses of nestling tarsus lengths, rather than the 13 cities stated in L354, as both Harjavalta and Barcelona are missing.
- The range of nestling tarsus lengths is from 10.2 to 25.9 mm, instead of the 10.2 to 26.2 mm reported in L356. Additionally, the range is 12.11 to 21.00 mm, rather than the 10 to 21.8 mm stated in L357.
- In Figure 2:
Table A appears to report inconsistent values. For example, the number of years sampled for the urban habitat in Budapest is 3, calculated as 2021 - 2019 + 1. However, for Antwerp, the reported number is 24, while using the same calculation, we should have 2018 - 1994 + 1, which equals 25 (See the table reproduced below and corresponding R code at the end of this file).
Overall (after having merged the three datasets), the range of years for the Urban area in Antwerp appears to be from 1999 to 2022, while the table claims it is from 1997 to 2022.
For Barcelona, the data indicate that the urban habitat was sampled from 1992 to 2018, whereas I found the range to be 1994-2018. Also, there are gaps in the sampling sequences between 2005 and 2009, as well as from 2014 to 2018. Therefore, I recommend reporting the actual number of years sampled instead of providing a range. For example, it would be more accurate to state the sampling period as "17: 1994 – 2018" for Barcelona.

City forest forest urban urban
1 ANT 25: 1994 – 2018 8 24: 1999 – 2022 1
3 BAR 23: 1998 – 2021 1 17: 1994 – 2018 2
5 BUD 2: 2020 – 2021 3 3: 2019 – 2021 3
7 GLA 4: 2017 – 2020 3 6: 2014 – 2020 2
9 GOT 3: 2017 – 2019 2 3: 2017 – 2019 2
11 HAR 25: 1991 – 2018 12 21: 1991 – 2018 7
13 HEL 3: 2018 – 2020 4 3: 2018 – 2020 2
15 MAL 7: 2013 – 2020 5 7: 2013 – 2020 4
17 MON 28: 1991 – 2018 1 8: 2007 – 2018 8
19 MUN 10: 2010 – 2019 12 2: 2014 – 2015 1
21 PAR 12: 2010 – 2021 5 9: 2012 – 2021 9
23 STR 8: 2014 – 2022 1 8: 2014 – 2022 8
25 VES 8: 2013 – 2020 2 8: 2013 – 2020 3
27 WAR 5: 2016 – 2020 2 5: 2016 – 2020 5

For the table C in Figure 2, I found an inconsistency between the dataset and the reported data for nestling tarsus measurements in blue tits. The correct number of individuals seems 2816 instead of 2905. This discrepancy may be related to the lower number of subpopulations I have identified, which is 25 instead of 26 in the forest. Additionally, it appears that the number of subpopulations has been inverted between forest and urban habitats for the Lay data in both great tit and blue tit (See the table reproduced below and corresponding R code at the end of this file).
StudySystems Individuals UrbanSubpops ForestSubpops
Adult tarsus - Great tit 13 13554 50 49
Adult tarsus - Blue tit 8 8414 30 28
Nestling tarsus - Great tit 12 7505 48 47
Nestling tarsus - Blue tit 8 2816 28 25
Lay date - Great tit 13 6378 53 58
Lay date - Blue tit 9 4419 36 37

• The code to reproduce Figure 3 and Table 1 is not included. Although the results presented in Table 1 can be retrieved in the model outputs, it would be helpful (and could have ease my work) to provide the code to reproduce both Table 1 and Figure 3, even in a basic, unstyled format such as what I did for Figure 2 and the chunks of table 1 attached to the review (see brmsRes function provided at the end of the review).
  For instance, you can display the results of the model for Adult Tarsus measure in Great tit using brmsRes(summary(Comp_tarsAdGT_Hab), "Adult tarsus GT").

By inspecting model output (using .rds file but not running the whole model) I found inconsistencies in Table 1: the sample size for adult tarsus model in blue tit is 8411 instead of 8414. Also, for the same model, some values do not match those reported in the manuscript (see Table 1 and the .html tables attached to the review especially AT_BT.html). While this discrepancy does not alter the overall conclusion, it may be attributed to differences in sample size. I suggest double-checking this.

########################### R Code ##############################

############################

merge the three datasets

############################

import the data

dtars<-read.csv("ThompsonEA_UrbanVp_AdultTarsus.csv")
dnestTars<-read.csv("ThompsonEA_UrbanVp_NstTarsus.csv")
dld <-read.csv("ThompsonEA_UrbanVp_layDate.csv")

Merge tarsus and nestling datasets

merged_1_2 <- merge(dtars, dnestTars,
by = c("City", "BreedingSeason", "ID_clu", "Habitat", "ID", "Species",
"mean_propurb_1000", "mean_propurb_100",
"mean_shdi_1000", "mean_shdi_100", "mean_lat",
"area", "years_clu"),
all = TRUE)

Merge the result with lay date dataset

AllData <- merge(merged_1_2, dld,
by = c("City", "BreedingSeason", "ID_clu", "Habitat", "ID", "Species",
"mean_propurb_1000", "mean_propurb_100",
"mean_shdi_1000", "mean_shdi_100", "mean_lat",
"area", "years_clu"),
all = TRUE)

############################

reproduce the table A from Fig 2

############################

Retrieve year range and number of years per City and Habitat

yearSummary <- aggregate(BreedingSeason ~ City + Habitat, data = AllData, FUN = function(x) {
rangeY <- range(x)
numY <- length(unique(x))
paste0(numY, ": ", rangeY[1], " – ", rangeY[2])
})

Count number of subpopulations per City and Habitat

clusterSummary <- aggregate(ID_clu ~ City + Habitat, data = AllData, FUN = function(x) length(unique(x)))
names(clusterSummary)[3] <- "N subpopulation"

Merge summary and subpopulation count

finalTab <- merge(yearSummary, clusterSummary, by = c("City", "Habitat"))

reshape to wide format to match display

wideTab <- reshape(finalTab,
idvar = "City",
timevar = "Habitat",
direction = "wide")
names(wideTab) <- gsub("BreedingSeason\.", "", names(wideTab))
names(wideTab) <- gsub("N subpopulation\.", "", names(wideTab))

Display the result

wideTab

############################

reproduce the table C from Fig 2

############################

a quick function to summarize the datasets

summarize <- function(df, species) {
dfSp <- subset(df, Species == species)
nCities <- length(unique(dfSp$City))
nInd <- nrow(dfSp)

nURBclu <- length(unique(dfSp$ID?klzzwxh:0097?clu[dfSp$Habitat == "urban"]))
nFORclu <- length(unique(dfSp$ID?klzzwxh:0098?clu[dfSp$Habitat == "forest"]))

return(c(StudySystems = nCities,
Individuals = nInd,
UrbanSubpops = nURBclu,
ForestSubpops = nFORclu))
}

Species

gT <- "PARMAJ"
bT <- "CYACAE"

Make the summary

summaryTab <- rbind(
summarize(dtars, gT),
summarize(dtars, bT),
summarize(dnestTars, gT),
summarize(dnestTars, bT),
summarize(dld, gT),
summarize(dld, bT)
)

Change row and column names

rownames(summaryTab) <- c(
"Adult tarsus - Great tit", "Adult tarsus - Blue tit",
"Nestling tarsus - Great tit", "Nestling tarsus - Blue tit",
"Lay date - Great tit", "Lay date - Blue tit"
)

Display the result

summaryTab

-------------------- a function to display brms output ---------

brmsRes <- function(model, traitLab = "Trait Summary") {
if (!requireNamespace("knitr", quietly = TRUE)) install.packages("knitr")
if (!requireNamespace("htmltools", quietly = TRUE)) install.packages("htmltools")

# a quick function to format CI
CI <- function(est, l95, u95, digits = 3) {
paste0(round(est, digits), " [", round(l95, digits), ", ", round(u95, digits), "]")
}

# ----------------------
# Fixed effects (mean)
# ----------------------
fixef <- model$fixed meanEffects <- fixef[!grepl("^sigma?klzzwxh:0123?", rownames(fixef)), ] rownames(meanEffects) <- gsub("^(Intercept)$", "Intercept (βm0)", rownames(meanEffects))
rownames(meanEffects) <- gsub("^Habitaturban$", "Habitat (urban)", rownames(meanEffects)) rownames(meanEffects) <- gsub("^meanlat.scaled$", "Latitude", rownames(meanEffects))
rownames(meanEffects) <- gsub("^SexM$", "Sex (male)", rownames(meanEffects))

meanFormatted <- data.frame(
Term = rownames(meanEffects),
Estimate_CI = mapply(CI, meanEffects$Estimate, meanEffects$l-95% CI, meanEffects$u-95% CI)
)

# ----------------------
# Random effects for the mean part
# ----------------------
rand <- model$random

randMeanFormatted <- data.frame(
Term = c("Year", "System", "Cluster: Forest", "Cluster: Urban"),
Estimate_CI = c(
CI(rand$BreedingSeason$Estimate, rand$BreedingSeason$l-95% CI, rand$BreedingSeason$u-95% CI),
CI(rand$City$Estimate[1], rand$City$l-95% CI[1], rand$City$u-95% CI[1]),
CI(rand$ID?klzzwxh:0150?clu$Estimate[1], rand$ID?klzzwxh:0151?clu$l-95% CI[1], rand$ID?klzzwxh:0152?clu$u-95% CI[1]),
CI(rand$ID?klzzwxh:0153?clu$Estimate[3], rand$ID?klzzwxh:0154?clu$l-95% CI[3], rand$ID?klzzwxh:0155?clu$u-95% CI[3])
)
)

# ----------------------
# Fixed Effects (Dispersion)
# ----------------------
dispEffects <- fixef[grepl("^sigma_", rownames(fixef)), , drop = FALSE]
rownames(dispEffects) <- gsub("^sigma_Intercept$", "Intercept (βv0,exp)", rownames(dispEffects)) rownames(dispEffects) <- gsub("^sigma?klzzwxh:0178?Habitaturban$", "Habitat (urban)", rownames(dispEffects))
rownames(dispEffects) <- gsub("^sigma_meanshdi_1000.scaled$", "Heterogeneity (1000m)", rownames(dispEffects)) rownames(dispEffects) <- gsub("^sigma?klzzwxh:0181?meanlat.scaled$", "Latitude", rownames(dispEffects))
rownames(dispEffects) <- gsub("^sigma_area.scaled$", "Cluster area", rownames(dispEffects)) rownames(dispEffects) <- gsub("^sigma?klzzwxh:0183?years?klzzwxh:0184?clu.scaled$", "Cluster years", rownames(dispEffects))
rownames(dispEffects) <- gsub("^sigma_SexM$", "A) Sex (male) C) Age (1)", rownames(dispEffects))

dispFormatted <- data.frame(
Term = rownames(dispEffects),
Estimate_CI = mapply(CI, dispEffects$Estimate, dispEffects$l-95% CI, dispEffects$u-95% CI)
)

# ----------------------
# Random effects for the dispersion part
# ----------------------
randDispFormatted <- data.frame(
Term = c("System (Intercept)", "System (rmean, dispersion)", "Forest (intercept)", "Forest (rmean, dispersion)", "Urban (intercept)", "Urban (rmean, dispersion)"),
Estimate_CI = c(
CI(rand$City$Estimate[2], rand$City$l-95% CI[2], rand$City$u-95% CI[2]),
CI(rand$City$Estimate[3], rand$City$l-95% CI[3], rand$City$u-95% CI[3]),
CI(rand$ID?klzzwxh:0216?clu$Estimate[2], rand$ID?klzzwxh:0217?clu$l-95% CI[2], rand$ID?klzzwxh:0218?clu$u-95% CI[2]),
CI(rand$ID?klzzwxh:0219?clu$Estimate[5], rand$ID?klzzwxh:0220?clu$l-95% CI[5], rand$ID?klzzwxh:0221?clu$u-95% CI[5]),
CI(rand$ID?klzzwxh:0222?clu$Estimate[4], rand$ID?klzzwxh:0223?clu$l-95% CI[4], rand$ID?klzzwxh:0224?clu$u-95% CI[4]),
CI(rand$ID?klzzwxh:0225?clu$Estimate[6], rand$ID?klzzwxh:0226?clu$l-95% CI[6], rand$ID?klzzwxh:0227?clu$u-95% CI[6])
)
)

# ---------------
# Build final HTML table
# ---------------
fullDat <- rbind(
data.frame(Term = "Fixed Effects (Mean part)", Estimate_CI = ""),
meanFormatted,
data.frame(Term = "Random Effects (Mean part)", Estimate_CI = ""),
randMeanFormatted,
data.frame(Term = "Fixed Effects (Dispersion part)", Estimate_CI = ""),
dispFormatted,
data.frame(Term = "Random Effects (Dispersion part)", Estimate_CI = ""),
randDispFormatted
)

# Print kable
htmltab <- knitr::kable(
fullDat,
caption = paste0(traitLab, " | ", "n = ", model$nobs),
col.names = c("Parameter", "Estimate [95% CI]"),
align = "l",
format = "html"
)
htmltools::browsable(htmltools::HTML(htmltab))
}

Thank you to the authors for carefully addressing my comments and suggestions. Based on their response, the improvements made to the code and data documentation now substantially enhance reproducibility.
One small note: the link to the Dryad repository currently appears to be broken—this may simply be due to the ongoing review process. It would be good to ensure that the link is active and accessible upon publication.
I'm happy to have contributed to this process and appreciate the authors' thoughtful responses.