##### Load in required packages ####
library(car)
library(lme4)
# for data visualisation
library(ggplot2)
library(ggpubr)
library(ggsci)

##### Load in and check data ####
RpData <- read.csv("Aphid_Fitness_Results_Rpadi.csv")
SaData <- read.csv("Aphid_Fitness_Results_Savenae.csv")

FitnessData <- subset(RpData, Dead %in% c ("Alive")) # Only include reps that are alive
FitnessDataSa <- subset(SaData, Dead %in% c ("Alive"))# Only include reps that are alive

# Manipulate data - make sure factors are coded as factors
FitnessDataSa$Population <- as.factor(FitnessDataSa$Population)
FitnessDataSa$Block <- as.factor(FitnessDataSa$Block)
FitnessDataSa$Morph <- as.factor(FitnessDataSa$Morph)
FitnessDataSa$Endosymbiont <- as.factor(FitnessDataSa$Endosymbiont)
FitnessDataSa$Genotype <- as.factor(FitnessDataSa$Genotype)

#
FitnessData$Population <- as.factor(FitnessData$Population)
FitnessData$Block <- as.factor(FitnessData$Block)
FitnessData$Morph <- as.factor(FitnessData$Morph)
FitnessData$Endosymbiont <- as.factor(FitnessData$Endosymbiont)
FitnessData$Genotype <- as.factor(FitnessData$Genotype)

# Calclate the rm metric
FitnessData$rm <- 0.74*((log2(FitnessData$Offspring_PRT))/FitnessData$Pre_Reproductive_Time)
FitnessDataSa$rm <- 0.74*((log2(FitnessDataSa$Offspring_PRT))/FitnessDataSa$Pre_Reproductive_Time)

####### Data pre-processing #####
## SA data
FitnessDataSa$Morph <- car::recode(FitnessDataSa$Morph, "'Winged' = 'Alate'; 'Unwinged' = 'Apterous'") # recode
FitnessDataSa$Morph <- factor(FitnessDataSa$Morph, levels = c("Apterous", "Alate"))

## RP data
FitnessData$Morph <- car::recode(FitnessData$Morph, "'Winged' = 'Alate'; 'Unwinged' = 'Apterous'")
FitnessData$Morph <- factor(FitnessData$Morph, levels = c("Apterous", "Alate"))

# Reorder data
summary(FitnessDataSa$Population)
FitnessDataSa$Population = factor(FitnessDataSa$Population, levels = c("SA-1", "SA-2", "SA-3", "SA-4", "SA-5",
                                                                       "SA-6", "SA-8", "SA-9", "SA-10", "SA-11", "SA-12", 
                                                                       "SA-13", "SA-14", "SA-15", "SA-16", "SA-17", "SA-18", 
                                                                       "SA-19", "SA-20", "SA-21", "SA-22", "SA-23", "SA-24", 
                                                                       "SA-25", "SA-26"))

FitnessData$Population = factor(FitnessData$Population, levels = c("RP-1", "RP-2", "RP-3", "RP-4", "RP-5",
                                                                   "RP-6", "RP-7", "RP-12"))


### DATA ANALYSIS #####
### SA analysis

## Analysis of winged moprh
MSWing <- glm(Winged_Binary ~ 
                Genotype + Genotype : Endosymbiont + Block, 
              family = "binomial"(link = "logit"), data = FitnessDataSa)
plot(MSWing)
summary(MSWing)
Anova(MSWing)
AIC(MSWing)

## Analysis of development time
# Lmer 
MS1.a <- lmer(Pre_Reproductive_Time ~ 
                (Genotype + Genotype : Endosymbiont) + Morph + Genotype : Morph  + 
                (1|Block), data = FitnessDataSa)
summary(MS1.a)
Anova(MS1.a)
plot(MS1.a)

## Analysis of rm
# Lmer
MS2.a <- lmer(rm ~ 
                (Genotype + Genotype : Endosymbiont) + Morph + Genotype : Morph  + 
                (1|Block), data = FitnessDataSa)
summary(MS2.a)
Anova(MS2.a)
plot(MS2.a)

#####
### RP analysis
## Analysis of winged moprh
MRWing <- glm(Winged_Binary ~ 
                Genotype + Genotype : Endosymbiont + Block, 
              family = "binomial"(link = "logit"), data = FitnessData)
plot(MRWing)
summary(MRWing)
Anova(MRWing)
AIC(MRWing)

## Analysis of development time
# Lmer
MR1.a <- lmer(Pre_Reproductive_Time ~ 
                (Genotype + Genotype : Endosymbiont) + Morph + Genotype : Morph +
                (1|Block), data = FitnessData)
summary(MR1.a)
Anova(MR1.a)
plot(MR1.a)

## Analysis of rm
# Lmer
MR2.a <- lmer(rm ~ (Genotype + Genotype : Endosymbiont) + Morph + Genotype : Morph +
                (1|Block), 
              data = FitnessData)
summary(MR2.a)
Anova(MR2.a)
plot(MR2.a)
AIC(MR2.a)

##### Plots #####
# Change levels in endosymbiont factor
FitnessDataSa$Endosymbiont <- factor(FitnessDataSa$Endosymbiont, levels = c("None", "F. symbiotica", "R. insecticola", "R. insecticola and F. symbiotica"))
#
FitnessData$Endosymbiont <- factor(FitnessData$Endosymbiont, levels = c("None", "F.symbiotica", "H.defensa"))

# Plot themes
PlotTheme <- theme_classic(base_size = 11)
Jitter_Theme <- geom_jitter(size = 0.3, alpha = 0.15)
Boxplot_Theme <- geom_boxplot(outlier.alpha = 0.00)

### Plot for S avenae - Pre reproductive development time
SaDev <- ggplot(FitnessDataSa, 
                aes(x = Genotype, y = Pre_Reproductive_Time, fill = Endosymbiont))+
  Boxplot_Theme+Jitter_Theme+PlotTheme+
  xlab("S. avenae genotype")+
  ylab(element_blank())+facet_wrap(~Morph)+
  theme(legend.position = "bottom", legend.text = element_text(size = 8))+
  scale_fill_manual(values = c("#8dd3c7","#ffffb3", "#bebada", "#fb8072"))
SaDev

### Plot for R padi - Pre reproductive development time
RpDev <- ggplot(FitnessData, 
                aes(x = Genotype, y = Pre_Reproductive_Time, fill = Endosymbiont))+
  Boxplot_Theme+Jitter_Theme+PlotTheme+
  xlab("R. padi genotype")+
  ylab(element_blank())+facet_wrap(~Morph)+
  theme(legend.position = "bottom", legend.text = element_text(size = 8))+
  scale_fill_manual(values = c("#8dd3c7","#ffffb3", "#febadf"))
RpDev

# Combine plots and export
AphidDev <- ggarrange(SaDev, RpDev, ncol = 1, nrow = 2, labels = "AUTO", font.label = list(size = 7))
AphidDev
AphidDev2 <- annotate_figure(AphidDev, 
                             left = text_grob("Development time (pre-reproductive period)", rot = 90, size = 10))
AphidDev2
ggsave("Fig_2_Final.png", AphidDev2, units = "mm", width = 160, height = 130)

### Plot for S avenae - rm
SaFec <- ggplot(FitnessDataSa, 
                aes(x = Genotype, y = rm, fill = Endosymbiont))+
  Boxplot_Theme+Jitter_Theme+PlotTheme+
  xlab("S. avenae genotype")+
  ylab(element_blank())+facet_wrap(~Morph)+
  theme(legend.position = "bottom", legend.text = element_text(size = 8))+
  scale_fill_manual(values = c("#8dd3c7","#ffffb3", "#bebada", "#fb8072"))
SaFec

### Plot for R padi - rm
RpFec <- ggplot(FitnessData, 
                aes(x = Genotype, y = rm, fill = Endosymbiont))+
  Boxplot_Theme+Jitter_Theme+PlotTheme+
  xlab("R. padi genotype")+
  ylab(element_blank())+facet_wrap(~Morph)+
  theme(legend.position = "bottom", legend.text = element_text(size = 8))+
  scale_fill_manual(values = c("#8dd3c7","#ffffb3", "#febadf"))
RpFec

# Combine plots and export
AphidFec <- ggarrange(SaFec, RpFec, ncol = 1, nrow = 2, labels = "AUTO", font.label = list(size = 7))
AphidFec
AphidFec2 <- annotate_figure(AphidFec, 
                             left = text_grob("Intrinsic rate of fitness (rm)", rot = 90, size = 10))
AphidFec2
ggsave("Fig_3_Final.png", AphidFec2, units = "mm", width = 160, height = 130)

### Alate production plot - load in dataframe of apterous aphid proportion for each genotype x endosymbiont combination
ApterousDat <- read.csv("Prop_Winged.csv")
ApterousDat$Endosymbiont <- factor(ApterousDat$Endosymbiont, levels = c("None", "F. symbiotica", "R. insecticola", "R. insecticola and F. symbiotica", "H. defensa"))

# divide in S avenae and R padi
SaW <- subset(WingedDat, Species %in% c ("S_avenae"))
RpW <- subset(WingedDat, Species %in% c ("R_padi"))

### Plot for S. avenae - alate proportion
SaWing <- ggplot(SaW, 
                 aes(x = Genotype, y = Prop_Winged, fill = Endosymbiont))+
  Boxplot_Theme+Jitter_Theme+PlotTheme+
  xlab("S. avenae genotype")+
  ylab(element_blank())+
  theme(legend.position = "bottom", legend.text = element_text(size = 8))+
  scale_fill_manual(values = c("#8dd3c7","#ffffb3", "#bebada", "#fb8072"))+
  ylim(0.0,1.0)+
  annotate("text", label = c("F.s"), x = 1, y = 1.0, size = 3)+
  annotate("text", label = c("R.i"), x = 2, y = 1.0, size = 3)+
  annotate("text", label = c("R.i"), x = 2.85, y = 1.0, size = 3)+
  annotate("text", label = c("R.i + F.s"), x = 3.20, y = 1.0, size = 3)+
  annotate("text", label = c("None"), x = 3.75, y = 1.0, size = 3)+
  annotate("text", label = c("R.i"), x = 4, y = 1.0, size = 3)+
  annotate("text", label = c("R.i + F.s"), x = 4.35, y = 1.0, size = 3)
SaWing

### Plot for R padi - alate proportion
RpWing <- ggplot(RpW, 
                 aes(x = Genotype, y = Prop_Winged, fill = Endosymbiont))+
  Boxplot_Theme+Jitter_Theme+PlotTheme+
  xlab("R. padi genotype")+
  ylab(element_blank())+
  theme(legend.position = "bottom", legend.text = element_text(size = 8))+
  scale_fill_manual(values = c("#8dd3c7","#ffffb3", "#febadf"))+
  ylim(0.0,1.0)+
  annotate("text", label = c("None"), x = 0.8, y = 1.0, size = 3)+
  annotate("text", label = c("F.s"), x = 1.2, y = 1.0, size = 3)+
  annotate("text", label = c("None"), x = 1.8, y = 1.0, size = 3)+
  annotate("text", label = c("H.d"), x = 2.2, y = 1.0, size = 3)+
  annotate("text", label = c("None"), x = 3, y = 1.00, size = 3)
RpWing

# Combine plots and export
AphidWing <- ggarrange(SaWing, RpWing, ncol = 1, nrow = 2, labels = "AUTO", font.label = list(size = 7))
AphidWing
AphidWing2 <- annotate_figure(AphidWing, 
                              left = text_grob("Proportion of apterous adults", rot = 90, size = 10))
AphidWing2
ggsave("Fig_1_Final.png", AphidWing2, units = "mm",width = 160, height = 130)