Modified from hybridViability.Rmd (Oct 2017)

Purpose

Test for trade-offs between male and female function in F2 hybrids of two sympatric moth-pollinated plant species, Schiedea kaalae and S. hookeri.

Data Import

setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
pol <- read.table("f2pollenEH.csv", header=T, sep="\t", colClasses=c(date="Date"))
pol <- pol[!is.na(pol$V1) & pol$dadfullcross!="9999",]

pol$v <- rowSums(pol[,grepl("V", colnames(pol))])/10 # average viable counts
pol$i <- rowSums(pol[,grepl("I", colnames(pol))])/10 # average inviable counts
pol$vpf <- (pol$v * 444.4/20) * 10      # number viable per flower (5 flrs)
pol$ipf <- (pol$i * 444.4/20) * 10      # number inviable per flower (5 flrs)
pol$vp <- pol$vpf / (pol$vpf + pol$ipf) # proportion viable
pol$tpf <- pol$vpf + pol$ipf
pol$t <- 10*(pol$v + pol$i)

polm2 <- pol[grepl("x",pol$momfullcross),]
polm2.split <- cbind(polm2$crossid, as.data.frame(matrix(unlist(matrix(lapply(matrix(unlist(strsplit(as.character(polm2$momfullcross), " x ", fixed=T)), nrow=2), function(x) { spl <- unlist(strsplit(unlist(strsplit(x, " ", fixed=T)), "-", fixed=T)); c(spl[1],gsub("-NA","",paste(spl[2],spl[3],spl[4], sep="-"))) }), ncol=2)), ncol=4, byrow=T)))
colnames(polm2.split) <- c("crossid","mommompop","mommomid","momdadpop","momdadid")

pold2 <- pol[grepl("x",pol$dadfullcross),]
pold2.split <- cbind(pold2$crossid, as.data.frame(matrix(unlist(matrix(lapply(matrix(unlist(strsplit(as.character(pold2$dadfullcross), " x ", fixed=T)), nrow=2), function(x) { spl <- unlist(strsplit(unlist(strsplit(x, " ", fixed=T)), "-", fixed=T)); c(spl[1],gsub("-NA","",paste(spl[2],spl[3],spl[4], sep="-"))) }), ncol=2)), ncol=4, byrow=T)))
colnames(pold2.split) <- c("crossid","dadmompop","dadmomid","daddadpop","daddadid")

polm1 <- pol[!grepl("x",pol$momfullcross),]
polm1.split <- cbind(polm1$crossid, as.data.frame(matrix(unlist(matrix(lapply(as.character(polm1$momfullcross), function(x) { spl <- unlist(strsplit(unlist(strsplit(x, " ", fixed=T)), "-", fixed=T)); c(spl[1],gsub("-NA","",paste(spl[2],spl[3],spl[4], sep="-"))) }), ncol=1)), ncol=2, byrow=T)))
colnames(polm1.split) <- c("crossid","mompop","momid")

pold1 <- pol[!grepl("x",pol$dadfullcross),]
pold1.split <- cbind(pold1$crossid, as.data.frame(matrix(unlist(matrix(lapply(as.character(pold1$dadfullcross), function(x) { spl <- unlist(strsplit(unlist(strsplit(x, " ", fixed=T)), "-", fixed=T)); c(spl[1],gsub("-NA","",paste(spl[2],spl[3],spl[4], sep="-"))) }), ncol=1)), ncol=2, byrow=T)))
colnames(pold1.split) <- c("crossid","dadpop","dadid")


pol <- merge(pol, polm2.split, by="crossid", all.x=T)
pol <- merge(pol, pold2.split, by="crossid", all.x=T)
pol <- merge(pol, polm1.split, by="crossid", all.x=T)
pol <- merge(pol, pold1.split, by="crossid", all.x=T)

popcols <- c("mompop","dadpop","mommompop","momdadpop","dadmompop","daddadpop")
pol[popcols] <- lapply(pol[popcols], function(x) { sapply(x, mapvalues, from = c("794","866","899","879","892","904","3587"), to = c("WK","WK","WK","879WKG","892WKG","904WPG","3587WP"), warn_missing=F) })

kaalpops <- c("904WPG","3587WP", "892WKG")
pop2sp <- function(x) { factor(ifelse(x %in% kaalpops, "kaal", ifelse(is.na(x),NA,"hook")))}
pol$mommomspecies <- pop2sp(pol$mommompop)
pol$momdadspecies <- pop2sp(pol$momdadpop)
pol$dadmomspecies <- pop2sp(pol$dadmompop)
pol$daddadspecies <- pop2sp(pol$daddadpop)
pol$momcross <- factor(toupper(paste0(substr(pol$mommomspecies,0,1), substr(pol$momdadspecies,0,1))))
pol$dadcross <- factor(toupper(paste0(substr(pol$dadmomspecies,0,1), substr(pol$daddadspecies,0,1))))

pol$momcross <- factor(ifelse(pol$momcross=="NANA", toupper(substr(as.character(pop2sp(pol$mompop)),0,1)), as.character(pol$momcross)))
pol$dadcross <- factor(ifelse(pol$dadcross=="NANA", toupper(substr(as.character(pop2sp(pol$dadpop)),0,1)), as.character(pol$dadcross)))
# is.na(pol$dadcross) <- is.na(pol$momcross) <- pol$crosstype=="control"

pol$dadcross <- factor(pol$dadcross, c("H",  "HK", "KH", "K"))
pol$momcross  <- factor(pol$momcross, c("H",  "HK", "KH", "K"))
pol$crosstype <- with(pol, factor(paste(momcross, dadcross, sep="x")))
pol$species <- factor(ifelse(pol$mompop %in% kaalpops, "kaal", "hook"))
pol$dadsp <-   factor(ifelse(pol$dadpop %in% kaalpops, "kaal", "hook"))
pol <- within(pol, smompop <- as.factor(paste(species,mompop,sep="")))
pol <- within(pol, mompid <- as.factor(paste(mompop,momid,sep=".")))
pol <- within(pol, dadpid <- as.factor(paste(dadpop,dadid,sep=".")))

pol$mompop <- as.character(pol$mompop)
pol$mompop[is.na(pol$mompop)] <- as.character(pol$momcross[is.na(pol$mompop)])
pol$mompop <- as.factor(pol$mompop)
pol$dadpop <- as.character(pol$dadpop)
pol$dadpop[is.na(pol$dadpop)] <- as.character(pol$dadcross[is.na(pol$dadpop)])
pol$dadpop <- as.factor(pol$dadpop)
setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
sds <- read.table("f2seedsEH.csv", header=T, sep="\t", colClasses=c(date_emasculation="Date",   date_pollination="Date"))

msds <- sds %>% 
  group_by(crossid,pollen_donor) %>% 
  dplyr::summarize(viable_seeds = mean(viable_seeds), inviable_seeds = mean(inviable_seeds)) %>%
  mutate(pviable_seeds = viable_seeds/(viable_seeds+inviable_seeds))

sdspol <- merge(msds,pol,by="crossid")

str(sdspol)
   'data.frame':    43 obs. of  59 variables:
    $ crossid       : Factor w/ 45 levels "101-1","111-1",..: 1 2 3 4 5 6 7 8 8 10 ...
    $ pollen_donor  : Factor w/ 2 levels "H","K": 1 2 1 2 1 1 2 2 1 2 ...
    $ viable_seeds  : num  12.29 1.2 1.67 1.75 7.33 ...
    $ inviable_seeds: num  0.571 1.8 1.333 0 1.667 ...
    $ pviable_seeds : num  0.956 0.4 0.556 1 0.815 ...
    $ date          : Date, format: "2018-07-23" "2018-08-22" ...
    $ momfullcross  : Factor w/ 25 levels "3587-14","3587-7 x 879 10-1 ID 1",..: 14 15 16 16 17 18 18 19 19 20 ...
    $ dadfullcross  : Factor w/ 29 levels "3587-7 x 879 G-2 ID 4",..: 4 9 15 20 11 21 21 19 19 18 ...
    $ V1            : int  32 26 19 32 37 5 14 38 38 15 ...
    $ I1            : int  9 7 1 10 19 25 14 18 18 12 ...
    $ V2            : int  75 40 25 38 27 3 18 37 37 13 ...
    $ I2            : int  12 5 5 17 19 25 15 26 26 10 ...
    $ V3            : int  67 37 22 31 25 8 10 44 44 21 ...
    $ I3            : int  7 9 0 8 11 18 13 29 29 17 ...
    $ V4            : int  46 32 32 38 81 7 14 44 44 13 ...
    $ I4            : int  11 6 4 6 53 26 10 22 22 16 ...
    $ V5            : int  66 29 14 39 29 3 4 30 30 15 ...
    $ I5            : int  7 3 9 6 20 16 6 19 19 11 ...
    $ V6            : int  74 33 15 33 21 9 11 46 46 17 ...
    $ I6            : int  8 8 2 7 11 21 9 21 21 17 ...
    $ V7            : int  62 26 12 38 20 6 12 38 38 18 ...
    $ I7            : int  8 9 3 6 10 13 10 27 27 17 ...
    $ V8            : int  52 31 26 30 30 4 11 37 37 16 ...
    $ I8            : int  11 10 5 10 20 14 7 10 10 13 ...
    $ V9            : int  61 30 23 32 28 8 12 49 49 16 ...
    $ I9            : int  10 5 6 6 12 20 6 17 17 16 ...
    $ V10           : int  53 29 22 48 24 8 9 35 35 19 ...
    $ I10           : int  6 7 2 9 24 21 7 16 16 18 ...
    $ v             : num  58.8 31.3 21 35.9 32.2 6.1 11.5 39.8 39.8 16.3 ...
    $ i             : num  8.9 6.9 3.7 8.5 19.9 19.9 9.7 20.5 20.5 14.7 ...
    $ vpf           : num  13065 6955 4666 7977 7155 ...
    $ ipf           : num  1978 1533 822 1889 4422 ...
    $ vp            : num  0.869 0.819 0.85 0.809 0.618 ...
    $ tpf           : num  15043 8488 5488 9866 11577 ...
    $ t             : num  677 382 247 444 521 260 212 603 603 310 ...
    $ mommompop     : Factor w/ 4 levels "3587WP","879WKG",..: 2 NA 3 3 3 3 3 3 3 NA ...
    $ mommomid      : Factor w/ 8 levels "1","10-1","2-2",..: 8 NA 1 1 4 4 4 5 5 NA ...
    $ momdadpop     : Factor w/ 3 levels "879WKG","892WKG",..: 3 NA 1 1 1 1 1 1 1 NA ...
    $ momdadid      : Factor w/ 10 levels "1-ID-1","10-1-ID",..: 7 NA 4 4 2 4 4 9 9 NA ...
    $ dadmompop     : Factor w/ 4 levels "3587WP","879WKG",..: NA 2 NA NA NA 3 3 3 3 2 ...
    $ dadmomid      : Factor w/ 10 levels "1","10","10-1",..: NA 3 NA NA NA 2 2 1 1 10 ...
    $ daddadpop     : Factor w/ 3 levels "879WKG","892WKG",..: NA 3 NA NA NA 1 1 1 1 3 ...
    $ daddadid      : Factor w/ 15 levels "10-1-ID","10-ID-7",..: NA 12 NA NA NA 1 1 3 3 11 ...
    $ mompop        : Factor w/ 6 levels "3587WP","879WKG",..: 5 3 6 6 6 6 6 6 6 4 ...
    $ momid         : Factor w/ 10 levels "1","10-1","10-1-C1",..: NA 1 NA NA NA NA NA NA NA 6 ...
    $ dadpop        : Factor w/ 7 levels "3587WP","879WKG",..: 7 5 2 3 2 6 6 6 6 5 ...
    $ dadid         : Factor w/ 10 levels "10-1","10-4/3/13",..: 4 NA 10 2 6 NA NA NA NA NA ...
    $ mommomspecies : Factor w/ 2 levels "hook","kaal": 1 NA 2 2 2 2 2 2 2 NA ...
    $ momdadspecies : Factor w/ 2 levels "hook","kaal": 2 NA 1 1 1 1 1 1 1 NA ...
    $ dadmomspecies : Factor w/ 2 levels "hook","kaal": NA 1 NA NA NA 2 2 2 2 1 ...
    $ daddadspecies : Factor w/ 2 levels "hook","kaal": NA 2 NA NA NA 1 1 1 1 2 ...
    $ momcross      : Factor w/ 4 levels "H","HK","KH",..: 2 4 3 3 3 3 3 3 3 4 ...
    $ dadcross      : Factor w/ 4 levels "H","HK","KH",..: 1 2 1 4 1 3 3 3 3 2 ...
    $ crosstype     : Factor w/ 10 levels "HKxH","HKxHK",..: 1 9 6 7 6 8 8 8 8 9 ...
    $ species       : Factor w/ 2 levels "hook","kaal": 1 2 1 1 1 1 1 1 1 2 ...
    $ dadsp         : Factor w/ 2 levels "hook","kaal": 1 1 1 2 1 1 1 1 1 1 ...
    $ smompop       : Factor w/ 5 levels "hook879WKG","hookNA",..: 2 4 2 2 2 2 2 2 2 5 ...
    $ mompid        : Factor w/ 11 levels "3587WP.14","879WKG.10-1",..: 11 8 11 11 11 11 11 11 11 9 ...
    $ dadpid        : Factor w/ 11 levels "3587WP.A","879WKG.10-1",..: 11 10 5 6 4 10 10 10 10 10 ...

Inventory

Take out F2s before analysis

table(sdspol$crosstype)
   
    HKxH HKxHK  HKxK  HxHK  HxKH  KHxH  KHxK KHxKH  KxHK  KxKH 
       3     9     2     5     3     3     4     8     3     3
sdspol <- sdspol %>% filter(!(crosstype %in% c("HKxHK","KHxKH")))

Analyses

Viable seeds vs. viable pollen

(sdspol.plot <- ggplot(sdspol, aes(x=vpf, y=viable_seeds, color=crosstype, linetype=pollen_donor, shape=pollen_donor)) + 
   geom_point(size=3) +
   #geom_smooth(method = "lm", se=F, aes(group=crosstype)) + 
   # geom_quantile(aes(group=1))+
   xlab("Viable pollen grains per flower") + 
   ylab("Mean number of viable seeds") +
   scale_color_brewer("Backcross Type", type="qual", palette="Paired") +
   #scale_shape_manual("Maternal Plant", values=c(15, 18, 16, 17))+
   scale_shape_manual("Pollen donor", values=c(15, 19))+
   theme_minimal()+
   theme(legend.key.size = unit(1.7,"lines")))

#ggsave("viableseeds_viablepollen.png", sdspol.plot, type="cairo", width=6, height=6)

Variables: viable pollen grains per flower, mean number of viable seeds.

  1. Test for an overall relationship at the crosstype level, taking the mean values for each crosstype (low N)
sdspol.mean <- sdspol %>% group_by(crosstype) %>% summarize_at(c("viable_seeds", "vpf"), mean)
sdspol.mean.mod <- lm(viable_seeds ~ vpf, data=sdspol.mean)
summary(aov(sdspol.mean.mod))
               Df Sum Sq Mean Sq F value Pr(>F)
   vpf          1  0.148  0.1476   0.057  0.819
   Residuals    6 15.459  2.5765
  1. Test for an overall relationship at the plant level, ignoring crosstype
sdspol.overall.mod <- lm(viable_seeds ~ vpf, data=sdspol)
summary(aov(sdspol.overall.mod))
               Df Sum Sq Mean Sq F value Pr(>F)  
   vpf          1  45.67   45.67   4.578 0.0428 *
   Residuals   24 239.43    9.98                 
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. Test for a relationship within each cross type
sdspol.mod <- lm(viable_seeds ~ vpf * crosstype, data=sdspol)
summary(aov(sdspol.mod))
                 Df Sum Sq Mean Sq F value  Pr(>F)   
   vpf            1  45.67   45.67  12.514 0.00538 **
   crosstype      7 101.11   14.44   3.957 0.02485 * 
   vpf:crosstype  7 101.83   14.55   3.986 0.02429 * 
   Residuals     10  36.50    3.65                   
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(sdspol.mod)
   
   Call:
   lm(formula = viable_seeds ~ vpf * crosstype, data = sdspol)
   
   Residuals:
        Min       1Q   Median       3Q      Max 
   -3.01790 -0.88314  0.08254  0.89785  2.02652 
   
   Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
   (Intercept)       -1.458e+01  4.083e+00  -3.572 0.005076 ** 
   vpf                2.029e-03  4.162e-04   4.875 0.000647 ***
   crosstypeHKxK      1.506e+01  4.638e+00   3.247 0.008760 ** 
   crosstypeHxHK      1.588e+01  5.373e+00   2.956 0.014375 *  
   crosstypeHxKH      8.514e+00  4.753e+00   1.791 0.103539    
   crosstypeKHxH      1.243e+01  5.706e+00   2.178 0.054385 .  
   crosstypeKHxK      1.487e+01  4.976e+00   2.988 0.013622 *  
   crosstypeKxHK      1.996e+01  6.037e+00   3.307 0.007921 ** 
   crosstypeKxKH      2.228e+01  5.434e+00   4.100 0.002146 ** 
   vpf:crosstypeHKxK -9.577e-04  8.926e-04  -1.073 0.308498    
   vpf:crosstypeHxHK -2.142e-03  6.142e-04  -3.488 0.005841 ** 
   vpf:crosstypeHxKH -3.249e-04  5.578e-04  -0.583 0.573112    
   vpf:crosstypeKHxH -7.853e-04  8.528e-04  -0.921 0.378769    
   vpf:crosstypeKHxK -1.876e-03  6.830e-04  -2.746 0.020613 *  
   vpf:crosstypeKxHK -2.477e-03  9.109e-04  -2.719 0.021592 *  
   vpf:crosstypeKxKH -2.888e-03  8.356e-04  -3.457 0.006159 ** 
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
   
   Residual standard error: 1.91 on 10 degrees of freedom
   Multiple R-squared:  0.872,  Adjusted R-squared:   0.68 
   F-statistic: 4.541 on 15 and 10 DF,  p-value: 0.01014