eduardofox2/Wu&al_D_citri_2018

## Wu&al_D_citri_2018
##Script prepared by Eduardo G. P. Fox & Wu Fengnian for a paper analysing infection by a pest insect##############
#Written on 2016 in P.R. China, Guangzhou

#Host Plant-Mediated Interactions Between ‘Candidatus Liberibacter asiaticus’ and Its Vector Diaphorina citri Kuwayama (Hemiptera: Liviidae)
#Published at https://doi.org/10.1093/jee/toy182

#Order of itens presented according with appearance on published paper.

#(1)Developmental period of different life stages of Diaphorina citri fed on healthy and “Candidatus Liberibacter asiaticus”-infected hosts

#Developmental period of Diaphorina citri on healthy Citrus reticulata cv. Mashuiju
Dev.CLas.neg.reticulata.Egg<-c(3.60, 3.40, 3.20,  4.10,  3.90,  3.50)
Dev.CLas.neg.reticulata.1st<-c(1.18, 1.57, 1.03,  1.80,  1.20,  2.10)
Dev.CLas.neg.reticulata.2nd<-c(1.50, 1.36, 1.57,  2.10,  1.48,  1.80)
Dev.CLas.neg.reticulata.3rd<-c(1.71, 1.99, 1.57,  1.84,  2.67,  2.35)
Dev.CLas.neg.reticulata.4th<-c(2.57, 2.46, 2.60,  2.68,  2.43,  2.85)
Dev.CLas.neg.reticulata.5th<-c(3.64, 4.02, 3.22,  4.26,  3.87,  3.25)
Dev.CLas.neg.reticulata.Total<-c(10.60, 11.38, 9.98, 12.68, 11.65, 12.35)


#Developmental period of Diaphorina citri on Clas-infected Citrus reticulata cv. Mashuiju
Dev.CLas.pos.reticulata.Egg<-c(4.00, 4.10, 3.90, 3.40, 4.10, 4.50)
Dev.CLas.pos.reticulata.1st<-c(1.11, 1.36, 1.00, 1.30, 1.20, 1.90)
Dev.CLas.pos.reticulata.2nd<-c(1.55, 1.06, 2.04, 1.80, 1.60, 2.30)
Dev.CLas.pos.reticulata.3rd<-c(1.55, 1.99, 1.03, 1.60, 2.00, 1.80)
Dev.CLas.pos.reticulata.4th<-c(1.77, 1.76, 1.80, 1.64, 1.86, 1.58)
Dev.CLas.pos.reticulata.5th<-c(1.77, 2.37, 1.37, 2.32, 3.13, 2.30)
Dev.CLas.pos.reticulata.Total<-c(7.75, 8.52, 7.22, 8.66, 9.79, 9.88)


#Developmental period of Diaphorina citri on healthy Citrus sunki Hort. ex Tanaka
Dev.CLas.neg.sunki.Egg<-c(4.20,  3.50,  5.20,  5.95,  5.91,  5.29)
Dev.CLas.neg.sunki.1st<-c(2.05,  1.90,  2.26,  2.55,  2.93,  2.00)
Dev.CLas.neg.sunki.2nd<-c(2.41,  2.24,  2.66,  1.04,  2.04,  2.86)
Dev.CLas.neg.sunki.3rd<-c(2.34,  2.57,  2.13,  2.96,  2.72,  3.43)
Dev.CLas.neg.sunki.4th<-c(2.87,  2.46,  3.27,  4.90,  4.17,  4.71)
Dev.CLas.neg.sunki.5th<-c(3.82,  4.13,  3.69,  3.52,  4.13,  4.14)
Dev.CLas.neg.sunki.Total<-c(13.50, 13.28, 13.99, 14.97, 15.99, 17.14)


#Developmental period of Diaphorina citri on Clas-infected Citrus sunki Hort. ex Tanaka
Dev.CLas.pos.sunki.Egg<-c(4.30,  4.20,  4.70,  8.03,  6.35,  8.12)
Dev.CLas.pos.sunki.1st<-c(1.91,  1.79,  2.04,  1.65,  2.29,  2.98)
Dev.CLas.pos.sunki.2nd<-c(2.07,  1.90,  2.37,  2.18,  2.27,  2.69)
Dev.CLas.pos.sunki.3rd<-c(2.23,  2.66,  1.90,  3.66,  2.10,  2.81)
Dev.CLas.pos.sunki.4th<-c(2.33,  2.04,  2.66,  3.14,  2.10,  3.00)
Dev.CLas.pos.sunki.5th<-c(2.71,  2.90,  2.66,  3.27,  3.40,  2.80)
Dev.CLas.pos.sunki.Total<-c(11.25, 11.28, 11.61, 13.90, 12.16, 14.28)


CLas.neg.reticulata.development<-data.frame(mget(ls(pattern="Dev.CLas.neg.reticulata")))
names(CLas.neg.reticulata.development)<-sub("Dev.CLas.neg.reticulata.","",names(CLas.neg.reticulata.development))

CLas.pos.reticulata.development<-data.frame(mget(ls(pattern="Dev.CLas.pos.reticulata")))
names(CLas.pos.reticulata.development)<-sub("Dev.CLas.pos.reticulata.","",names(CLas.pos.reticulata.development))

CLas.neg.sunki.development<-data.frame(mget(ls(pattern="Dev.CLas.neg.sunki")))
names(CLas.neg.sunki.development)<-sub("Dev.CLas.neg.sunki.","",names(CLas.neg.sunki.development))

CLas.pos.sunki.development<-data.frame(mget(ls(pattern="Dev.CLas.pos.sunki")))
names(CLas.pos.sunki.development)<-sub("Dev.CLas.pos.sunki.","",names(CLas.pos.sunki.development))

Development<-list(CLas.neg.reticulata.development,CLas.pos.reticulata.development, CLas.neg.sunki.development, CLas.pos.sunki.development)
names(Development)<-c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki")


#(2)Survival ratio of different life stages of Diaphorina citri fed on healthy and “Candidatus Liberibacter asiaticus”-infected hosts

#Survival ratio of Diaphorina citri on healthy Citrus reticulata cv. Mashuiju
Surv.CLas.neg.reticulata.Egg<-c(86.11, 83.33, 80.28, 92.98, 85.32, 82.35)
Surv.CLas.neg.reticulata.1st<-c(72.28, 78.75, 82.29, 83.25, 75.32, 83.60)
Surv.CLas.neg.reticulata.2nd<-c(80.28, 97.60, 89.58, 93.26, 84.26, 90.26)
Surv.CLas.neg.reticulata.3rd<-c(76.49, 92.00, 77.99, 88.65, 73.37, 76.69)
Surv.CLas.neg.reticulata.4th<-c(76.92, 79.86, 85.71, 89.68, 91.26, 82.70)
Surv.CLas.neg.reticulata.5th<-c(80.00, 87.50, 91.67, 86.32, 80.26, 88.36)
Surv.CLas.neg.reticulata.Total<-c(27.31, 49.41, 45.17, 53.29, 34.10, 42.28)

#Survival ratio of Diaphorina citri on Clas-infected Citrus reticulata cv. Mashuiju
Surv.CLas.pos.reticulata.Egg<-c(88.57, 79.43, 73.58, 75.87, 81.42, 77.36)
Surv.CLas.pos.reticulata.1st<-c(69.57, 76.67, 77.50, 74.26, 79.23, 68.37)
Surv.CLas.pos.reticulata.2nd<-c(87.50, 82.61, 74.19, 80.25, 79.37, 86.34)
Surv.CLas.pos.reticulata.3rd<-c(92.86, 84.21, 93.33, 87.37, 91.37, 92.25)
Surv.CLas.pos.reticulata.4th<-c(90.91, 91.00, 93.75, 98.50, 91.37, 93.35)
Surv.CLas.pos.reticulata.5th<-c(96.00, 98.63, 97.01, 91.35, 97.68, 88.56)
Surv.CLas.pos.reticulata.Total<-c(49.33, 47.87, 48.81, 46.84, 51.27, 45.02)


#Survival ratio of Diaphorina citri on healthy Citrus sunki Hort. ex Tanaka
Surv.CLas.neg.sunki.Egg<-c(89.74, 85.52, 83.55, 88.28, 88.89, 86.00)
Surv.CLas.neg.sunki.1st<-c(75.71, 86.44, 89.76, 85.71, 87.83, 83.33)
Surv.CLas.neg.sunki.2nd<-c(89.55, 96.08, 93.57, 83.08, 97.78, 76.79)
Surv.CLas.neg.sunki.3rd<-c(80.00, 86.73, 84.58, 83.78, 88.46, 78.14)
Surv.CLas.neg.sunki.4th<-c(87.50, 91.76, 89.66, 80.65, 79.98, 84.00)
Surv.CLas.neg.sunki.5th<-c(79.17, 83.33, 80.25, 96.00, 92.80, 92.89)
Surv.CLas.neg.sunki.Total<-c(37.57, 55.08, 51.11, 46.19, 56.38, 39.01)

#Survival ratio of Diaphorina citri on Clas-infected Citrus sunki Hort. ex Tanaka
Surv.CLas.pos.sunki.Egg<-c(78.41, 71.56, 78.35, 75.27, 80.00, 85.19)
Surv.CLas.pos.sunki.1st<-c(68.42, 72.22, 73.08, 77.50, 95.65, 80.77)
Surv.CLas.pos.sunki.2nd<-c(76.92, 89.74, 86.84, 75.64, 95.45, 85.71)
Surv.CLas.pos.sunki.3rd<-c(90.48, 85.71, 93.94, 94.02, 87.14, 76.00)
Surv.CLas.pos.sunki.4th<-c(100.00, 95.45, 96.77, 91.67, 100.00, 75.56)
Surv.CLas.pos.sunki.5th<-c(100.00, 95.24, 89.47, 100.00, 95.36, 100.00)
Surv.CLas.pos.sunki.Total<-c(47.62, 50.51, 51.62, 50.52, 75.87, 39.75)

CLas.neg.reticulata.survival<-data.frame(mget(ls(pattern="Surv.CLas.neg.reticulata")))
names(CLas.neg.reticulata.survival)<-sub("Surv.CLas.neg.reticulata.","",names(CLas.neg.reticulata.survival))

CLas.pos.reticulata.survival<-data.frame(mget(ls(pattern="Surv.CLas.pos.reticulata")))
names(CLas.pos.reticulata.survival)<-sub("Surv.CLas.pos.reticulata.","",names(CLas.pos.reticulata.survival))

CLas.neg.sunki.survival<-data.frame(mget(ls(pattern="Surv.CLas.neg.sunki")))
names(CLas.neg.sunki.survival)<-sub("Surv.CLas.neg.sunki.","",names(CLas.neg.sunki.survival))

CLas.pos.sunki.survival<-data.frame(mget(ls(pattern="Surv.CLas.pos.sunki")))
names(CLas.pos.sunki.survival)<-sub("Surv.CLas.pos.sunki.","",names(CLas.pos.sunki.survival))

Survival<-list(CLas.neg.reticulata.survival,CLas.pos.reticulata.survival, CLas.neg.sunki.survival, CLas.pos.sunki.survival)
names(Survival)<-c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki")


#(3)Fecundity
#Fecundity of D. citri from CLas- C.reticulata
Fec.CLas.neg.reticulata.5d <-c(69, 93, 108, 55, 45, 70)
Fec.CLas.neg.reticulata.10d <-c(145, 143, 171, 174, 188, 146)
Fec.CLas.neg.reticulata.15d <-c(197, 203, 244, 216, 254, 195)
Fec.CLas.neg.reticulata.20d <-c(234, 218, 278, 288, 271, 247)
Fec.CLas.neg.reticulata.25d <-c(301, 289, 322, 292, 287, 285)
Fec.CLas.neg.reticulata.30d <-c(312, 302, 335, 296, 326, 348)

#Fecundity of D. citri from Clas+ C.reticulata
Fec.CLas.pos.reticulata.5d <-c(73, 96, 102, 64, 76, 57)
Fec.CLas.pos.reticulata.10d<-c(156, 162, 197, 196, 145, 186)
Fec.CLas.pos.reticulata.15d<-c(198, 227, 254, 252, 262, 267)
Fec.CLas.pos.reticulata.20d<-c(295, 284, 357, 361, 316, 349)
Fec.CLas.pos.reticulata.25d<-c(363, 341, 383, 392, 349, 381)
Fec.CLas.pos.reticulata.30d<-c(384, 352, 413, 398, 349, 401)

#Fecundity of D. citri from CLas- C. sunki
Fec.CLas.neg.sunki.5d <-c(99, 69, 95, 73, 76, 64)
Fec.CLas.neg.sunki.10d<-c(154, 184, 132, 122, 198, 183)
Fec.CLas.neg.sunki.15d<-c(213, 196, 174, 228, 293, 291)
Fec.CLas.neg.sunki.20d<-c(294, 291, 248, 354, 433, 396)
Fec.CLas.neg.sunki.25d<-c(340, 353, 293, 395, 472, 482)
Fec.CLas.neg.sunki.30d<-c(346, 363, 391, 496, 482, 502)


#Fecundity of D. citri from Clas+ C. sunki
Fec.CLas.pos.sunki.5d <-c(57, 67, 102, 86, 78, 65)
Fec.CLas.pos.sunki.10d<-c(162, 108, 218, 244, 218, 205)
Fec.CLas.pos.sunki.15d<-c(221, 201, 234, 383, 297, 304)
Fec.CLas.pos.sunki.20d<-c(288, 341, 298, 494, 377, 401)
Fec.CLas.pos.sunki.25d<-c(367, 366, 390, 576, 490, 549)
Fec.CLas.pos.sunki.30d<-c(452, 387, 460, 625, 530, 553)


CLas.neg.reticulata.fecundity<-data.frame(mget(ls(pattern="Fec.CLas.neg.reticulata")))
names(CLas.neg.reticulata.fecundity)<-sub("Fec.CLas.neg.reticulata.","",names(CLas.neg.reticulata.fecundity))

CLas.pos.reticulata.fecundity<-data.frame(mget(ls(pattern="Fec.CLas.pos.reticulata")))
names(CLas.pos.reticulata.fecundity)<-sub("Fec.CLas.pos.reticulata.","",names(CLas.pos.reticulata.fecundity))

CLas.neg.sunki.fecundity<-data.frame(mget(ls(pattern="Fec.CLas.neg.sunki")))
names(CLas.neg.sunki.fecundity)<-sub("Fec.CLas.neg.sunki.","",names(CLas.neg.sunki.fecundity))

CLas.pos.sunki.fecundity<-data.frame(mget(ls(pattern="Fec.CLas.pos.sunki")))
names(CLas.pos.sunki.fecundity)<-sub("Fec.CLas.pos.sunki.","",names(CLas.pos.sunki.fecundity))

Fecundity<-list(CLas.neg.reticulata.fecundity,CLas.pos.reticulata.fecundity, CLas.neg.sunki.fecundity, CLas.pos.sunki.fecundity)
names(Fecundity)<-c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki")


#(4)Starvation ability
#Survival rate from Clas-plant
Star.CLas.neg.8h <-c(95, 100, 95, 100, 95, 95)
Star.CLas.neg.16h<-c(80, 75, 65, 75, 70, 75)
Star.CLas.neg.24h<-c(50, 55, 40, 50, 50, 55)
Star.CLas.neg.32h<-c(45, 35, 40, 30, 20, 40)
Star.CLas.neg.40h<-c(10, 15, 15, 25, 15, 25)
Star.CLas.neg.48h<-c(5, 0, 10, 10, 0, 10)
Star.CLas.neg.56h<-c(0, 0, 0, 0, 0, 0)

#Survival rate from Clas+plant
Star.CLas.pos.8h <-c(75, 95, 90, 90, 100, 90)
Star.CLas.pos.16h<-c(50, 90, 40, 50, 65, 70)
Star.CLas.pos.24h<-c(30, 50, 15, 35, 35, 35)
Star.CLas.pos.32h<-c(5, 40, 0, 20, 5, 15)
Star.CLas.pos.40h<-c(0, 10, 0, 0, 0, 5)
Star.CLas.pos.48h<-c(0, 0, 0, 0, 0, 0)
Star.CLas.pos.56h<-c(0, 0, 0, 0, 0, 0)


Star.CLas.neg<-data.frame(mget(ls(pattern="Star.CLas.neg.")))
names(Star.CLas.neg)<-sub("Star.CLas.neg.","",names(Star.CLas.neg))

Star.CLas.pos<-data.frame(mget(ls(pattern="Star.CLas.pos.")))
names(Star.CLas.pos)<-sub("Star.CLas.pos.","",names(Star.CLas.pos))

Starvation<-list(Star.CLas.neg, Star.CLas.pos)
names(Starvation)<-c("CLas-","CLas+")


#(5)Overwintering ability
#Survival rate of overwintering D. citri on Clas-plants
Over.CLas.neg.Nov <-c(100, 100, 100, 100, 100, 100)
Over.CLas.neg.Dec <-c(62, 58, 64, 80, 56, 60)
Over.CLas.neg.Jan_next<-c(38, 30, 34, 32, 38, 40)
Over.CLas.neg.Feb_next <-c(16, 26, 20, 22, 14, 24)

#Survival rate of overwintering D. citri on Clas+plants
Over.CLas.pos.Nov <-c(100, 100, 100, 100, 100, 100)
Over.CLas.pos.Dec <-c(48, 50, 60, 70, 52, 56)
Over.CLas.pos.Jan_next<-c(12, 10, 16, 14, 18, 22)
Over.CLas.pos.Feb_next <-c(2, 8, 6, 10, 14, 18)


Over.CLas.neg<-data.frame(mget(ls(pattern="Over.CLas.neg.")))
names(Over.CLas.neg)<-sub("Over.CLas.neg.","",names(Over.CLas.neg))

Over.CLas.pos<-data.frame(mget(ls(pattern="Over.CLas.pos.")))
names(Over.CLas.pos)<-sub("Over.CLas.pos.","",names(Over.CLas.pos))

Overwinter<-list(Over.CLas.neg, Over.CLas.pos)
names(Overwinter)<-c("CLas-","CLas+")

#You're left with 5 objects containing all data: Development, Fecundity, Survival, Starvation, Overwinter; last two contains only 1 host species.

#Preparing for plotting and analysis
require(plyr)
require(reshape2)
require(ggplot2)

long.Development<-melt(Development)
long.Fecundity<-melt(Fecundity)
long.Survival<-melt(Survival)
long.Starvation<-melt(Starvation)
long.Overwinter<-melt(Overwinter)

summary.Development<-ddply(long.Development, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))
summary.Fecundity<-ddply(long.Fecundity, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))
summary.Survival<-ddply(long.Survival, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))
summary.Starvation<-ddply(long.Starvation, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))
summary.Overwinter<-ddply(long.Overwinter, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))


long.Development$variable<-factor(long.Development$variable, levels=c("Egg", "1st","2nd", "3rd","4th","5th", "Total"))
long.Development$L1<-factor(long.Development$L1, levels=c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki"))
long.Fecundity$variable<-factor(long.Fecundity$variable, levels=c("5d", "10d","15d", "20d","25d","30d"))
long.Fecundity$L1<-factor(long.Fecundity$L1, levels=c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki"))
long.Survival$variable<-factor(long.Survival$variable, levels=c("Egg", "1st","2nd", "3rd","4th","5th", "Total"))
long.Survival$L1<-factor(long.Survival$L1, levels=c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki"))
long.Overwinter$variable<-factor(long.Overwinter$variable, levels=c("Nov","Dec","Jan_next","Feb_next"))
long.Overwinter$L1<-factor(long.Overwinter$L1, levels=c("CLas-","CLas+"))

summary.Fecundity$variable<-factor(summary.Fecundity$variable, levels=c("5d", "10d","15d", "20d","25d","30d"))
summary.Fecundity$L1<-factor(summary.Fecundity$L1, levels=c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki"))


#From this point all data is ready for analysis and plot

#Pushed line plots as suggested but using ggplot for better final image construction
line.colours<-c(rep("#E9695F",5), rep("#4F8FFA",5), rep("#00AC22",5),rep("#dde85d"))
line.colours2<-c(rep("#E9695F",5), rep("#4F8FFA",5))
line.colours3<-c(rep("#E9695F",5), rep("#4F8FFA",5))


#Development times
ggplot(long.Development, aes(x = variable, y =value, fill=(L1))) +
 geom_boxplot(width=0.6)+scale_fill_manual(values=c("#E9695F","#4F8FFA","#00AC22","#dde85d"))+
 theme_classic()+
 theme (legend.position=c(0.2,0.8)) +
 theme (legend.title=element_blank())+
 xlab("Life Stage Period") +
 ylab("Number of days")
x11()


 #number of eggs laid MADE INTO A CURVE
 Fig.1A<-ggplot(summary.Fecundity[seq(1,24,2),], aes(x = variable, y=mean, group=L1)) +
 geom_errorbar(aes(ymin=summary.Fecundity[seq(1,24,2),]$mean - summary.Fecundity[seq(1,24,2),]$se, ymax=summary.Fecundity[seq(1,24,2),]$mean + summary.Fecundity[seq(1,24,2),]$se), width = .1,, lwd=1.0) +
 geom_line(size=2, aes(colour=L1))+
 geom_point(aes(x = variable, y=mean, group=L1))+
 scale_y_continuous(limits=c(0,550), breaks=c(0,100,200,300,400,500))+
 theme_classic()+
 theme (legend.position=c(0.2,0.8)) +
 theme (legend.title=element_blank())+
 xlab("Days after finishing preoviposition period (d)") +
 ylab("Cumulative mean no. of eggs per female")


Fig.1B<-ggplot(summary.Fecundity[seq(2,24,2),], aes(x = variable, y=mean, group=L1)) +
 geom_errorbar(aes(ymin=summary.Fecundity[seq(2,24,2),]$mean - summary.Fecundity[seq(2,24,2),]$se, ymax=summary.Fecundity[seq(2,24,2),]$mean + summary.Fecundity[seq(2,24,2),]$se), width = .1,, lwd=1.0) +
 geom_line(size=2, aes(colour=L1))+
 geom_point(aes(x = variable, y=mean, group=L1))+
 scale_y_continuous(limits=c(0,550), breaks=c(0,100,200,300,400,500))+
 theme_classic()+
 theme (legend.position=c(0.2,0.8)) +
 theme (legend.title=element_blank())+
 xlab("Days after finishing preoviposition period (d)") +
 ylab("Cumulative mean no. of eggs per female")


# Multiple plot function obtained from R CookBook Website
#
# ggplot objects can be passed in, or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

#Plotting images with aesthetics corrections
multiplot(Fig.1A+theme(text = element_text(size = 10))+ theme(legend.title = element_blank()),Fig.1B+theme(text = element_text(size = 10))+ theme(legend.title = element_blank()),cols=2)


> wilcox.test(long.Fecundity$value[1:36], long.Fecundity$value[73:108])


> wilcox.test(long.Fecundity$value[37:72], long.Fecundity$value[108:144])
	##Script prepared by Eduardo G. P. Fox & Wu Fengnian for a paper analysing infection by a pest insect##############
	#Written on 2016 in P.R. China, Guangzhou

	#Host Plant-Mediated Interactions Between ‘Candidatus Liberibacter asiaticus’ and Its Vector Diaphorina citri Kuwayama (Hemiptera: Liviidae)
	#Published at https://doi.org/10.1093/jee/toy182

	#Order of itens presented according with appearance on published paper.

	#(1)Developmental period of different life stages of Diaphorina citri fed on healthy and “Candidatus Liberibacter asiaticus”-infected hosts

	#Developmental period of Diaphorina citri on healthy Citrus reticulata cv. Mashuiju
	Dev.CLas.neg.reticulata.Egg<-c(3.60, 3.40, 3.20, 4.10, 3.90, 3.50)
	Dev.CLas.neg.reticulata.1st<-c(1.18, 1.57, 1.03, 1.80, 1.20, 2.10)
	Dev.CLas.neg.reticulata.2nd<-c(1.50, 1.36, 1.57, 2.10, 1.48, 1.80)
	Dev.CLas.neg.reticulata.3rd<-c(1.71, 1.99, 1.57, 1.84, 2.67, 2.35)
	Dev.CLas.neg.reticulata.4th<-c(2.57, 2.46, 2.60, 2.68, 2.43, 2.85)
	Dev.CLas.neg.reticulata.5th<-c(3.64, 4.02, 3.22, 4.26, 3.87, 3.25)
	Dev.CLas.neg.reticulata.Total<-c(10.60, 11.38, 9.98, 12.68, 11.65, 12.35)


	#Developmental period of Diaphorina citri on Clas-infected Citrus reticulata cv. Mashuiju
	Dev.CLas.pos.reticulata.Egg<-c(4.00, 4.10, 3.90, 3.40, 4.10, 4.50)
	Dev.CLas.pos.reticulata.1st<-c(1.11, 1.36, 1.00, 1.30, 1.20, 1.90)
	Dev.CLas.pos.reticulata.2nd<-c(1.55, 1.06, 2.04, 1.80, 1.60, 2.30)
	Dev.CLas.pos.reticulata.3rd<-c(1.55, 1.99, 1.03, 1.60, 2.00, 1.80)
	Dev.CLas.pos.reticulata.4th<-c(1.77, 1.76, 1.80, 1.64, 1.86, 1.58)
	Dev.CLas.pos.reticulata.5th<-c(1.77, 2.37, 1.37, 2.32, 3.13, 2.30)
	Dev.CLas.pos.reticulata.Total<-c(7.75, 8.52, 7.22, 8.66, 9.79, 9.88)


	#Developmental period of Diaphorina citri on healthy Citrus sunki Hort. ex Tanaka
	Dev.CLas.neg.sunki.Egg<-c(4.20, 3.50, 5.20, 5.95, 5.91, 5.29)
	Dev.CLas.neg.sunki.1st<-c(2.05, 1.90, 2.26, 2.55, 2.93, 2.00)
	Dev.CLas.neg.sunki.2nd<-c(2.41, 2.24, 2.66, 1.04, 2.04, 2.86)
	Dev.CLas.neg.sunki.3rd<-c(2.34, 2.57, 2.13, 2.96, 2.72, 3.43)
	Dev.CLas.neg.sunki.4th<-c(2.87, 2.46, 3.27, 4.90, 4.17, 4.71)
	Dev.CLas.neg.sunki.5th<-c(3.82, 4.13, 3.69, 3.52, 4.13, 4.14)
	Dev.CLas.neg.sunki.Total<-c(13.50, 13.28, 13.99, 14.97, 15.99, 17.14)



	#Developmental period of Diaphorina citri on Clas-infected Citrus sunki Hort. ex Tanaka
	Dev.CLas.pos.sunki.Egg<-c(4.30, 4.20, 4.70, 8.03, 6.35, 8.12)
	Dev.CLas.pos.sunki.1st<-c(1.91, 1.79, 2.04, 1.65, 2.29, 2.98)
	Dev.CLas.pos.sunki.2nd<-c(2.07, 1.90, 2.37, 2.18, 2.27, 2.69)
	Dev.CLas.pos.sunki.3rd<-c(2.23, 2.66, 1.90, 3.66, 2.10, 2.81)
	Dev.CLas.pos.sunki.4th<-c(2.33, 2.04, 2.66, 3.14, 2.10, 3.00)
	Dev.CLas.pos.sunki.5th<-c(2.71, 2.90, 2.66, 3.27, 3.40, 2.80)
	Dev.CLas.pos.sunki.Total<-c(11.25, 11.28, 11.61, 13.90, 12.16, 14.28)



	CLas.neg.reticulata.development<-data.frame(mget(ls(pattern="Dev.CLas.neg.reticulata")))
	names(CLas.neg.reticulata.development)<-sub("Dev.CLas.neg.reticulata.","",names(CLas.neg.reticulata.development))

	CLas.pos.reticulata.development<-data.frame(mget(ls(pattern="Dev.CLas.pos.reticulata")))
	names(CLas.pos.reticulata.development)<-sub("Dev.CLas.pos.reticulata.","",names(CLas.pos.reticulata.development))

	CLas.neg.sunki.development<-data.frame(mget(ls(pattern="Dev.CLas.neg.sunki")))
	names(CLas.neg.sunki.development)<-sub("Dev.CLas.neg.sunki.","",names(CLas.neg.sunki.development))

	CLas.pos.sunki.development<-data.frame(mget(ls(pattern="Dev.CLas.pos.sunki")))
	names(CLas.pos.sunki.development)<-sub("Dev.CLas.pos.sunki.","",names(CLas.pos.sunki.development))

	Development<-list(CLas.neg.reticulata.development,CLas.pos.reticulata.development, CLas.neg.sunki.development, CLas.pos.sunki.development)
	names(Development)<-c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki")




	#(2)Survival ratio of different life stages of Diaphorina citri fed on healthy and “Candidatus Liberibacter asiaticus”-infected hosts

	#Survival ratio of Diaphorina citri on healthy Citrus reticulata cv. Mashuiju
	Surv.CLas.neg.reticulata.Egg<-c(86.11, 83.33, 80.28, 92.98, 85.32, 82.35)
	Surv.CLas.neg.reticulata.1st<-c(72.28, 78.75, 82.29, 83.25, 75.32, 83.60)
	Surv.CLas.neg.reticulata.2nd<-c(80.28, 97.60, 89.58, 93.26, 84.26, 90.26)
	Surv.CLas.neg.reticulata.3rd<-c(76.49, 92.00, 77.99, 88.65, 73.37, 76.69)
	Surv.CLas.neg.reticulata.4th<-c(76.92, 79.86, 85.71, 89.68, 91.26, 82.70)
	Surv.CLas.neg.reticulata.5th<-c(80.00, 87.50, 91.67, 86.32, 80.26, 88.36)
	Surv.CLas.neg.reticulata.Total<-c(27.31, 49.41, 45.17, 53.29, 34.10, 42.28)

	#Survival ratio of Diaphorina citri on Clas-infected Citrus reticulata cv. Mashuiju
	Surv.CLas.pos.reticulata.Egg<-c(88.57, 79.43, 73.58, 75.87, 81.42, 77.36)
	Surv.CLas.pos.reticulata.1st<-c(69.57, 76.67, 77.50, 74.26, 79.23, 68.37)
	Surv.CLas.pos.reticulata.2nd<-c(87.50, 82.61, 74.19, 80.25, 79.37, 86.34)
	Surv.CLas.pos.reticulata.3rd<-c(92.86, 84.21, 93.33, 87.37, 91.37, 92.25)
	Surv.CLas.pos.reticulata.4th<-c(90.91, 91.00, 93.75, 98.50, 91.37, 93.35)
	Surv.CLas.pos.reticulata.5th<-c(96.00, 98.63, 97.01, 91.35, 97.68, 88.56)
	Surv.CLas.pos.reticulata.Total<-c(49.33, 47.87, 48.81, 46.84, 51.27, 45.02)


	#Survival ratio of Diaphorina citri on healthy Citrus sunki Hort. ex Tanaka
	Surv.CLas.neg.sunki.Egg<-c(89.74, 85.52, 83.55, 88.28, 88.89, 86.00)
	Surv.CLas.neg.sunki.1st<-c(75.71, 86.44, 89.76, 85.71, 87.83, 83.33)
	Surv.CLas.neg.sunki.2nd<-c(89.55, 96.08, 93.57, 83.08, 97.78, 76.79)
	Surv.CLas.neg.sunki.3rd<-c(80.00, 86.73, 84.58, 83.78, 88.46, 78.14)
	Surv.CLas.neg.sunki.4th<-c(87.50, 91.76, 89.66, 80.65, 79.98, 84.00)
	Surv.CLas.neg.sunki.5th<-c(79.17, 83.33, 80.25, 96.00, 92.80, 92.89)
	Surv.CLas.neg.sunki.Total<-c(37.57, 55.08, 51.11, 46.19, 56.38, 39.01)

	#Survival ratio of Diaphorina citri on Clas-infected Citrus sunki Hort. ex Tanaka
	Surv.CLas.pos.sunki.Egg<-c(78.41, 71.56, 78.35, 75.27, 80.00, 85.19)
	Surv.CLas.pos.sunki.1st<-c(68.42, 72.22, 73.08, 77.50, 95.65, 80.77)
	Surv.CLas.pos.sunki.2nd<-c(76.92, 89.74, 86.84, 75.64, 95.45, 85.71)
	Surv.CLas.pos.sunki.3rd<-c(90.48, 85.71, 93.94, 94.02, 87.14, 76.00)
	Surv.CLas.pos.sunki.4th<-c(100.00, 95.45, 96.77, 91.67, 100.00, 75.56)
	Surv.CLas.pos.sunki.5th<-c(100.00, 95.24, 89.47, 100.00, 95.36, 100.00)
	Surv.CLas.pos.sunki.Total<-c(47.62, 50.51, 51.62, 50.52, 75.87, 39.75)

	CLas.neg.reticulata.survival<-data.frame(mget(ls(pattern="Surv.CLas.neg.reticulata")))
	names(CLas.neg.reticulata.survival)<-sub("Surv.CLas.neg.reticulata.","",names(CLas.neg.reticulata.survival))

	CLas.pos.reticulata.survival<-data.frame(mget(ls(pattern="Surv.CLas.pos.reticulata")))
	names(CLas.pos.reticulata.survival)<-sub("Surv.CLas.pos.reticulata.","",names(CLas.pos.reticulata.survival))

	CLas.neg.sunki.survival<-data.frame(mget(ls(pattern="Surv.CLas.neg.sunki")))
	names(CLas.neg.sunki.survival)<-sub("Surv.CLas.neg.sunki.","",names(CLas.neg.sunki.survival))

	CLas.pos.sunki.survival<-data.frame(mget(ls(pattern="Surv.CLas.pos.sunki")))
	names(CLas.pos.sunki.survival)<-sub("Surv.CLas.pos.sunki.","",names(CLas.pos.sunki.survival))

	Survival<-list(CLas.neg.reticulata.survival,CLas.pos.reticulata.survival, CLas.neg.sunki.survival, CLas.pos.sunki.survival)
	names(Survival)<-c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki")


	#(3)Fecundity
	#Fecundity of D. citri from CLas- C.reticulata
	Fec.CLas.neg.reticulata.5d <-c(69, 93, 108, 55, 45, 70)
	Fec.CLas.neg.reticulata.10d <-c(145, 143, 171, 174, 188, 146)
	Fec.CLas.neg.reticulata.15d <-c(197, 203, 244, 216, 254, 195)
	Fec.CLas.neg.reticulata.20d <-c(234, 218, 278, 288, 271, 247)
	Fec.CLas.neg.reticulata.25d <-c(301, 289, 322, 292, 287, 285)
	Fec.CLas.neg.reticulata.30d <-c(312, 302, 335, 296, 326, 348)

	#Fecundity of D. citri from Clas+ C.reticulata
	Fec.CLas.pos.reticulata.5d <-c(73, 96, 102, 64, 76, 57)
	Fec.CLas.pos.reticulata.10d<-c(156, 162, 197, 196, 145, 186)
	Fec.CLas.pos.reticulata.15d<-c(198, 227, 254, 252, 262, 267)
	Fec.CLas.pos.reticulata.20d<-c(295, 284, 357, 361, 316, 349)
	Fec.CLas.pos.reticulata.25d<-c(363, 341, 383, 392, 349, 381)
	Fec.CLas.pos.reticulata.30d<-c(384, 352, 413, 398, 349, 401)

	#Fecundity of D. citri from CLas- C. sunki
	Fec.CLas.neg.sunki.5d <-c(99, 69, 95, 73, 76, 64)
	Fec.CLas.neg.sunki.10d<-c(154, 184, 132, 122, 198, 183)
	Fec.CLas.neg.sunki.15d<-c(213, 196, 174, 228, 293, 291)
	Fec.CLas.neg.sunki.20d<-c(294, 291, 248, 354, 433, 396)
	Fec.CLas.neg.sunki.25d<-c(340, 353, 293, 395, 472, 482)
	Fec.CLas.neg.sunki.30d<-c(346, 363, 391, 496, 482, 502)


	#Fecundity of D. citri from Clas+ C. sunki
	Fec.CLas.pos.sunki.5d <-c(57, 67, 102, 86, 78, 65)
	Fec.CLas.pos.sunki.10d<-c(162, 108, 218, 244, 218, 205)
	Fec.CLas.pos.sunki.15d<-c(221, 201, 234, 383, 297, 304)
	Fec.CLas.pos.sunki.20d<-c(288, 341, 298, 494, 377, 401)
	Fec.CLas.pos.sunki.25d<-c(367, 366, 390, 576, 490, 549)
	Fec.CLas.pos.sunki.30d<-c(452, 387, 460, 625, 530, 553)


	CLas.neg.reticulata.fecundity<-data.frame(mget(ls(pattern="Fec.CLas.neg.reticulata")))
	names(CLas.neg.reticulata.fecundity)<-sub("Fec.CLas.neg.reticulata.","",names(CLas.neg.reticulata.fecundity))

	CLas.pos.reticulata.fecundity<-data.frame(mget(ls(pattern="Fec.CLas.pos.reticulata")))
	names(CLas.pos.reticulata.fecundity)<-sub("Fec.CLas.pos.reticulata.","",names(CLas.pos.reticulata.fecundity))

	CLas.neg.sunki.fecundity<-data.frame(mget(ls(pattern="Fec.CLas.neg.sunki")))
	names(CLas.neg.sunki.fecundity)<-sub("Fec.CLas.neg.sunki.","",names(CLas.neg.sunki.fecundity))

	CLas.pos.sunki.fecundity<-data.frame(mget(ls(pattern="Fec.CLas.pos.sunki")))
	names(CLas.pos.sunki.fecundity)<-sub("Fec.CLas.pos.sunki.","",names(CLas.pos.sunki.fecundity))

	Fecundity<-list(CLas.neg.reticulata.fecundity,CLas.pos.reticulata.fecundity, CLas.neg.sunki.fecundity, CLas.pos.sunki.fecundity)
	names(Fecundity)<-c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki")



	#(4)Starvation ability
	#Survival rate from Clas-plant
	Star.CLas.neg.8h <-c(95, 100, 95, 100, 95, 95)
	Star.CLas.neg.16h<-c(80, 75, 65, 75, 70, 75)
	Star.CLas.neg.24h<-c(50, 55, 40, 50, 50, 55)
	Star.CLas.neg.32h<-c(45, 35, 40, 30, 20, 40)
	Star.CLas.neg.40h<-c(10, 15, 15, 25, 15, 25)
	Star.CLas.neg.48h<-c(5, 0, 10, 10, 0, 10)
	Star.CLas.neg.56h<-c(0, 0, 0, 0, 0, 0)

	#Survival rate from Clas+plant
	Star.CLas.pos.8h <-c(75, 95, 90, 90, 100, 90)
	Star.CLas.pos.16h<-c(50, 90, 40, 50, 65, 70)
	Star.CLas.pos.24h<-c(30, 50, 15, 35, 35, 35)
	Star.CLas.pos.32h<-c(5, 40, 0, 20, 5, 15)
	Star.CLas.pos.40h<-c(0, 10, 0, 0, 0, 5)
	Star.CLas.pos.48h<-c(0, 0, 0, 0, 0, 0)
	Star.CLas.pos.56h<-c(0, 0, 0, 0, 0, 0)


	Star.CLas.neg<-data.frame(mget(ls(pattern="Star.CLas.neg.")))
	names(Star.CLas.neg)<-sub("Star.CLas.neg.","",names(Star.CLas.neg))

	Star.CLas.pos<-data.frame(mget(ls(pattern="Star.CLas.pos.")))
	names(Star.CLas.pos)<-sub("Star.CLas.pos.","",names(Star.CLas.pos))

	Starvation<-list(Star.CLas.neg, Star.CLas.pos)
	names(Starvation)<-c("CLas-","CLas+")



	#(5)Overwintering ability
	#Survival rate of overwintering D. citri on Clas-plants
	Over.CLas.neg.Nov <-c(100, 100, 100, 100, 100, 100)
	Over.CLas.neg.Dec <-c(62, 58, 64, 80, 56, 60)
	Over.CLas.neg.Jan_next<-c(38, 30, 34, 32, 38, 40)
	Over.CLas.neg.Feb_next <-c(16, 26, 20, 22, 14, 24)

	#Survival rate of overwintering D. citri on Clas+plants
	Over.CLas.pos.Nov <-c(100, 100, 100, 100, 100, 100)
	Over.CLas.pos.Dec <-c(48, 50, 60, 70, 52, 56)
	Over.CLas.pos.Jan_next<-c(12, 10, 16, 14, 18, 22)
	Over.CLas.pos.Feb_next <-c(2, 8, 6, 10, 14, 18)


	Over.CLas.neg<-data.frame(mget(ls(pattern="Over.CLas.neg.")))
	names(Over.CLas.neg)<-sub("Over.CLas.neg.","",names(Over.CLas.neg))

	Over.CLas.pos<-data.frame(mget(ls(pattern="Over.CLas.pos.")))
	names(Over.CLas.pos)<-sub("Over.CLas.pos.","",names(Over.CLas.pos))

	Overwinter<-list(Over.CLas.neg, Over.CLas.pos)
	names(Overwinter)<-c("CLas-","CLas+")

	#You're left with 5 objects containing all data: Development, Fecundity, Survival, Starvation, Overwinter; last two contains only 1 host species.

	#Preparing for plotting and analysis
	require(plyr)
	require(reshape2)
	require(ggplot2)

	long.Development<-melt(Development)
	long.Fecundity<-melt(Fecundity)
	long.Survival<-melt(Survival)
	long.Starvation<-melt(Starvation)
	long.Overwinter<-melt(Overwinter)

	summary.Development<-ddply(long.Development, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))
	summary.Fecundity<-ddply(long.Fecundity, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))
	summary.Survival<-ddply(long.Survival, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))
	summary.Starvation<-ddply(long.Starvation, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))
	summary.Overwinter<-ddply(long.Overwinter, c("variable","L1"), summarise, min=min(value, na.rm=T), max =max(value, na.rm=t),mean=mean(value, na.rm=T), sd=sd(value, na.rm=T), se = sd/sqrt(length(value)))



	long.Development$variable<-factor(long.Development$variable, levels=c("Egg", "1st","2nd", "3rd","4th","5th", "Total"))
	long.Development$L1<-factor(long.Development$L1, levels=c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki"))
	long.Fecundity$variable<-factor(long.Fecundity$variable, levels=c("5d", "10d","15d", "20d","25d","30d"))
	long.Fecundity$L1<-factor(long.Fecundity$L1, levels=c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki"))
	long.Survival$variable<-factor(long.Survival$variable, levels=c("Egg", "1st","2nd", "3rd","4th","5th", "Total"))
	long.Survival$L1<-factor(long.Survival$L1, levels=c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki"))
	long.Overwinter$variable<-factor(long.Overwinter$variable, levels=c("Nov","Dec","Jan_next","Feb_next"))
	long.Overwinter$L1<-factor(long.Overwinter$L1, levels=c("CLas-","CLas+"))

	summary.Fecundity$variable<-factor(summary.Fecundity$variable, levels=c("5d", "10d","15d", "20d","25d","30d"))
	summary.Fecundity$L1<-factor(summary.Fecundity$L1, levels=c("CLas- C. reticulata","CLas+ C. reticulata","CLas- C. sunki","CLas+ C. sunki"))


	#From this point all data is ready for analysis and plot

	#Pushed line plots as suggested but using ggplot for better final image construction
	line.colours<-c(rep("#E9695F",5), rep("#4F8FFA",5), rep("#00AC22",5),rep("#dde85d"))
	line.colours2<-c(rep("#E9695F",5), rep("#4F8FFA",5))
	line.colours3<-c(rep("#E9695F",5), rep("#4F8FFA",5))



	#Development times
	ggplot(long.Development, aes(x = variable, y =value, fill=(L1))) +
	geom_boxplot(width=0.6)+scale_fill_manual(values=c("#E9695F","#4F8FFA","#00AC22","#dde85d"))+
	theme_classic()+
	theme (legend.position=c(0.2,0.8)) +
	theme (legend.title=element_blank())+
	xlab("Life Stage Period") +
	ylab("Number of days")
	x11()


	#number of eggs laid MADE INTO A CURVE
	Fig.1A<-ggplot(summary.Fecundity[seq(1,24,2),], aes(x = variable, y=mean, group=L1)) +
	geom_errorbar(aes(ymin=summary.Fecundity[seq(1,24,2),]$mean - summary.Fecundity[seq(1,24,2),]$se, ymax=summary.Fecundity[seq(1,24,2),]$mean + summary.Fecundity[seq(1,24,2),]$se), width = .1,, lwd=1.0) +
	geom_line(size=2, aes(colour=L1))+
	geom_point(aes(x = variable, y=mean, group=L1))+
	scale_y_continuous(limits=c(0,550), breaks=c(0,100,200,300,400,500))+
	theme_classic()+
	theme (legend.position=c(0.2,0.8)) +
	theme (legend.title=element_blank())+
	xlab("Days after finishing preoviposition period (d)") +
	ylab("Cumulative mean no. of eggs per female")


	Fig.1B<-ggplot(summary.Fecundity[seq(2,24,2),], aes(x = variable, y=mean, group=L1)) +
	geom_errorbar(aes(ymin=summary.Fecundity[seq(2,24,2),]$mean - summary.Fecundity[seq(2,24,2),]$se, ymax=summary.Fecundity[seq(2,24,2),]$mean + summary.Fecundity[seq(2,24,2),]$se), width = .1,, lwd=1.0) +
	geom_line(size=2, aes(colour=L1))+
	geom_point(aes(x = variable, y=mean, group=L1))+
	scale_y_continuous(limits=c(0,550), breaks=c(0,100,200,300,400,500))+
	theme_classic()+
	theme (legend.position=c(0.2,0.8)) +
	theme (legend.title=element_blank())+
	xlab("Days after finishing preoviposition period (d)") +
	ylab("Cumulative mean no. of eggs per female")


	# Multiple plot function obtained from R CookBook Website
	#
	# ggplot objects can be passed in, or to plotlist (as a list of ggplot objects)
	# - cols: Number of columns in layout
	# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
	#
	# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
	# then plot 1 will go in the upper left, 2 will go in the upper right, and
	# 3 will go all the way across the bottom.
	#
	multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
	library(grid)

	# Make a list from the ... arguments and plotlist
	plots <- c(list(...), plotlist)

	numPlots = length(plots)

	# If layout is NULL, then use 'cols' to determine layout
	if (is.null(layout)) {
	# Make the panel
	# ncol: Number of columns of plots
	# nrow: Number of rows needed, calculated from # of cols
	layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
	ncol = cols, nrow = ceiling(numPlots/cols))
	}

	if (numPlots==1) {
	print(plots[[1]])

	} else {
	# Set up the page
	grid.newpage()
	pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

	# Make each plot, in the correct location
	for (i in 1:numPlots) {
	# Get the i,j matrix positions of the regions that contain this subplot
	matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

	print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
	layout.pos.col = matchidx$col))
	}
	}
	}

	#Plotting images with aesthetics corrections
	multiplot(Fig.1A+theme(text = element_text(size = 10))+ theme(legend.title = element_blank()),Fig.1B+theme(text = element_text(size = 10))+ theme(legend.title = element_blank()),cols=2)



	> wilcox.test(long.Fecundity$value[1:36], long.Fecundity$value[73:108])


	> wilcox.test(long.Fecundity$value[37:72], long.Fecundity$value[108:144])