grosbeak ancestry associations

0.1 Investigate genotype/phenotype associations

Code
library(ggplot2)
library(ggExtra)
#load dataset from publicly available repository
samps<-read.csv("https://raw.githubusercontent.com/DevonDeRaad/grosbeak.rad.hybridtransect/refs/heads/main/data/grosbeak.sample.info.csv")

#plot overall phenotype against ancestry for males
#isolate samples with "male total phenotype score" (113 out of 138)
male.tot<-samps[!is.na(samps$male.total),]
#plot
plot(male.tot$mac.mel.q,male.tot$male.total)

Code
plot(male.tot$mac.mel.q,male.tot$tarsus)

Code
plot(male.tot$mac.mel.q,male.tot$wing)

Code
plot(male.tot$mac.mel.q,male.tot$nape)

Code
plot(male.tot$mac.mel.q,male.tot$back)

Code
plot(male.tot$mac.mel.q,male.tot$rump)

Code
plot(male.tot$mac.mel.q,male.tot$breast.underwing.coverts)

Code
plot(male.tot$mac.mel.q,male.tot$flanks)

0.2 Check to make sure that date collected is not strongly influencing testis size

Code
#subset to only males
male<-samps[samps$sex == "male",]
#with testis measurements
male<-samps[!is.na(samps$Rtestiswidth),]

#plot
plot(male$day.year, male$total.testis.area)
#calculate stats
cor(male$day.year,male$total.testis.area)
[1] 0.1060319
Code
lm<-lm(total.testis.area ~ day.year, data = male)
summary(lm)

Call:
lm(formula = total.testis.area ~ day.year, data = male)

Residuals:
     Min       1Q   Median       3Q      Max 
-120.490  -30.779   -3.303   35.510  131.594 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  68.9564   100.4290   0.687    0.494
day.year      0.7290     0.6489   1.123    0.264

Residual standard error: 49.11 on 111 degrees of freedom
Multiple R-squared:  0.01124,   Adjusted R-squared:  0.002335 
F-statistic: 1.262 on 1 and 111 DF,  p-value: 0.2637
Code
#plot with best fit linear regression line laid over
p<-ggplot(male, aes(x=day.year, y=total.testis.area)) +
      geom_point(cex = 3, alpha = .5)+
      guides(shape = guide_legend(override.aes = list(size = 5), order=2, label.theme= element_text(face="italic")))+
      xlab(paste("Day of the year (0 - 365)"))+
      ylab(paste("Testis area"))+
      labs(title = "") +
      scale_color_manual(values = c("#ef3b2c","#fff319"),na.value = "grey") +
      #ylim(c(0,12)) +
      #scale_y_continuous(breaks=seq(0, 1, 1/4)) +
      geom_abline(slope = coef(lm)[["day.year"]], intercept = coef(lm)[["(Intercept)"]]) +
      theme_classic()+
      theme(legend.position="none")
      
ggMarginal(p, type = "histogram") 

Code
#save
pp<-ggMarginal(p, type = "histogram") 
#ggsave("~/Desktop/grosbeak.rad/testis.date.plot.pdf", pp, width = 6,height = 4,units = "in")

0.3 Now test for ancestry / testis area association

Code
#calculate 'minor parent ancestry'
x<-c()
for (i in 1:nrow(male)){
  x[i]<-min(male$mac.lud.q[i],male$mac.mel.q[i])
}
hist(x, breaks=20)

Code
male$minor.ancestry<-x

#plot the association between 'minor parent ancestry' (essentially admixture proportion) and testis size
plot(male$minor.ancestry, male$total.testis.area)

Code
#try binning for hybrid vs not
mean(male$total.testis.area[male$minor.ancestry <.02])
[1] 180.4512
Code
mean(male$total.testis.area[male$minor.ancestry >.02])
[1] 184.871
Code
#try different cutoff
mean(male$total.testis.area[male$minor.ancestry <.2])
[1] 181.1869
Code
mean(male$total.testis.area[male$minor.ancestry >.2])
[1] 190.1667
Code
#try separating by class and species
mean(male$total.testis.area[male$mac.lud.q < .05])
[1] 171.5167
Code
mean(male$total.testis.area[male$mac.lud.q > .05 & male$mac.lud.q < .5])
[1] 186.9231
Code
mean(male$total.testis.area[male$mac.lud.q > .5 & male$mac.lud.q < .95])
[1] 172.8
Code
mean(male$total.testis.area[male$mac.lud.q > .95])
[1] 198.3714
Code
#make a column assigning samples as pure parental or hybrid
male$class<-"x"
male$class[male$mac.lud.q < .05]<-"melanocephalus"
male$class[male$mac.lud.q > .05 & male$mac.lud.q < .5]<-"melanocephalus hybrid"
male$class[male$mac.lud.q > .5 & male$mac.lud.q < .95]<-"ludovicianus hybrid"
male$class[male$mac.lud.q > .95]<-"ludovicianus"

#plot the boxplots
ggplot(male, aes(class, total.testis.area)) + 
  geom_boxplot() + 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Code
#boxplot(total.testis.area ~ class, data = male)
#no apparent differences based on ancestry class

#plot individual testis measurements
plot(male$minor.ancestry, male$Ltestislength)

Code
plot(male$minor.ancestry, male$Ltestiswidth)

Code
plot(male$minor.ancestry, male$Rtestislength)

Code
plot(male$minor.ancestry, male$Rtestiswidth)

Code
#calculate stats
cor(male$minor.ancestry,male$total.testis.area)
[1] -0.02463896
Code
lm<-lm(total.testis.area ~ minor.ancestry, data = male)
summary(lm)

Call:
lm(formula = total.testis.area ~ minor.ancestry, data = male)

Residuals:
     Min       1Q   Median       3Q      Max 
-122.146  -30.146   -3.945   33.854  132.854 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     182.146      5.002   36.41   <2e-16 ***
minor.ancestry  -13.964     53.776   -0.26    0.796    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 49.38 on 111 degrees of freedom
Multiple R-squared:  0.0006071, Adjusted R-squared:  -0.008396 
F-statistic: 0.06743 on 1 and 111 DF,  p-value: 0.7956

0.4 Plot final figure

Code
#plot with best fit linear regression line laid over
p<-ggplot(male, aes(x=minor.ancestry, y=total.testis.area)) +
      geom_point(cex = 3, alpha = .5)+
      guides(shape = guide_legend(override.aes = list(size = 5), order=2, label.theme= element_text(face="italic")))+
      xlab(paste("Minor parent ancestry proportion"))+
      ylab(paste("Testis area"))+
      labs(title = "") +
      scale_color_manual(values = c("#ef3b2c","#fff319"),na.value = "grey") +
      #ylim(c(0,12)) +
      #scale_y_continuous(breaks=seq(0, 1, 1/4)) +
      geom_abline(slope = coef(lm)[["minor.ancestry"]], intercept = coef(lm)[["(Intercept)"]]) +
      theme_classic()+
      theme(legend.position="none")
      
ggMarginal(p, type = "histogram")

Code
#save
pp<-ggMarginal(p, type = "histogram") 
#ggsave("~/Desktop/grosbeak.rad/testis.ancestry.plot.pdf", pp, width = 4,height = 4,units = "in")