Setting the working directory.

Length-based auximetric analysis

Growth data

We start importing the data. These data were extracted from Fishbase using the package rfishbase on June 2020. Only species from brackish and salted waters were included. Von Bertalanffy growth parameters, namely the rate coefficient \(K\) and asymptotic length \(L_\infty\) were imported from the Popgrowth table. As recommended by Fishbase, stocks for which the estimated \(L_\infty\) was outside the range 0.66-1.33 times the maximum length ever observed \(L_{max}\) were considered dubious and thus excluded. Likewise, stocks whose estimates were outside auximetric ellipses were removed. Notice though that a large number of stocks were not tested and were kept in the analysis. Total length was used for all species except for Scombridaes for which fork length was considered.

growthdata<-read.csv("VB_Growth_data_Fishbase_2020_June.csv")
summary(growthdata)
##     AutoCtr         SpecCode       Species            StockCode    
##  Min.   :  139   Min.   :    4   Length:5517        Min.   :    4  
##  1st Qu.: 2177   1st Qu.:  204   Class :character   1st Qu.:  408  
##  Median : 4703   Median : 1288   Mode  :character   Median : 1672  
##  Mean   : 7815   Mean   : 3828                      Mean   : 7248  
##  3rd Qu.:10235   3rd Qu.: 4449                      3rd Qu.: 4986  
##  Max.   :23914   Max.   :67314                      Max.   :58558  
##                                                                    
##    Locality          PopGrowthRef        Sex                 Loo         
##  Length:5517        Min.   :     9   Length:5517        Min.   :   2.10  
##  Class :character   1st Qu.:  2081   Class :character   1st Qu.:  25.90  
##  Mode  :character   Median : 13304   Mode  :character   Median :  41.00  
##                     Mean   : 38231                      Mean   :  62.46  
##                     3rd Qu.: 74354                      3rd Qu.:  77.20  
##                     Max.   :119692                      Max.   :1966.00  
##                                                                          
##    TLinfinity            K                 to            Winfinity       
##  Min.   :   2.10   Min.   : 0.0070   Min.   :-12.210   Min.   :       0  
##  1st Qu.:  26.80   1st Qu.: 0.1590   1st Qu.: -1.700   1st Qu.:     193  
##  Median :  41.97   Median : 0.2700   Median : -0.858   Median :     721  
##  Mean   :  63.92   Mean   : 0.4517   Mean   : -1.200   Mean   :   31059  
##  3rd Qu.:  79.00   3rd Qu.: 0.5300   3rd Qu.: -0.300   3rd Qu.:    4526  
##  Max.   :1966.00   Max.   :12.9570   Max.   :  1.100   Max.   :20000000  
##                                      NA's   :2883      NA's   :3091      
##      Type              Auxim              LinfLmax   LogKLogLoo     
##  Length:5517        Length:5517        Min.   :0   Min.   :-0.8340  
##  Class :character   Class :character   1st Qu.:0   1st Qu.:-0.4215  
##  Mode  :character   Mode  :character   Median :0   Median :-0.3149  
##                                        Mean   :0   Mean   :-0.2625  
##                                        3rd Qu.:0   3rd Qu.:-0.1297  
##                                        Max.   :0   Max.   : 0.9892  
##                                                    NA's   :2659     
##        tm                Lm            LmSex              TypeLm         
##  Min.   :  0.200   Min.   :  1.50   Length:5517        Length:5517       
##  1st Qu.:  2.000   1st Qu.: 16.50   Class :character   Class :character  
##  Median :  4.000   Median : 27.10   Mode  :character   Mode  :character  
##  Mean   :  6.315   Mean   : 41.77                                        
##  3rd Qu.:  7.200   3rd Qu.: 49.80                                        
##  Max.   :156.000   Max.   :922.00                                        
##  NA's   :5340      NA's   :4766                                          
##      LmMale          LmFemale      GrowthEnviron     
##  Min.   :  6.00   Min.   :  8.00   Length:5517       
##  1st Qu.: 19.85   1st Qu.: 20.40   Class :character  
##  Median : 29.00   Median : 34.70   Mode  :character  
##  Mean   : 48.76   Mean   : 49.44                     
##  3rd Qu.: 50.00   3rd Qu.: 62.50                     
##  Max.   :280.00   Max.   :257.00                     
##  NA's   :5410     NA's   :5360

There are more NAs for TLinfinity than for for K which means that for some stocks the former is missing. We thus keep only stocks with both parameters available.

nrow(growthdata[!is.na(growthdata$TLinfinity),])
## [1] 5517
length(unique(growthdata[!is.na(growthdata$TLinfinity),]$Species))
## [1] 1313

We have 5517 stocks for 1313 species with an estimate of of K and TLinfinity.

Auximetric analysis

We start by producing a log-log plot of K versus TLinfinity. According to Pauly and Munro (1984), \(log_{10}(K) \sim log_{10}(L_\infty)\) is expected to follow a regression line with a slope of 2, the intercept of which gives the growth index of the different species.

We start by adding an identifier, speciestype for 3 groups of species: tunas that are among the fish species with the higshest growth performances, deep-sea sharks that have life-history strategies relativeley close to that of coelancanths, and roughies that are among the fish with the lowest growth performances. For this, we need to load package rfishbase and use function species_list.

library(rfishbase)
growthdata$tunas<-growthdata$Species %in% species_list(Genus="Thunnus") | growthdata$Species %in% species_list(Genus="Katsuwonus") | growthdata$Species %in% species_list(Genus="Euthynnus") | growthdata$Species %in% species_list(Genus="Auxis") | growthdata$Species %in% species_list(Genus="Allothunnus")

growthdata$deepseasharks<- growthdata$Species %in% species_list(Order="Squaliformes")

growthdata$roughies<-growthdata$Species %in% species_list(Family="Trachichthyidae")

growthdata$speciestype[growthdata$tunas==T]<-"Tunas"
growthdata$speciestype[growthdata$deepseasharks==T]<-"Deep-sea sharks" 

growthdata$speciestype[growthdata$roughies==T]<-"Roughies"

growthdata$speciestype[growthdata$tunas==F & growthdata$deepseasharks==F & growthdata$roughies==F]<-"Other species"
growthdata$speciestype<-as.factor(growthdata$speciestype)
summary(growthdata)
##     AutoCtr         SpecCode       Species            StockCode    
##  Min.   :  139   Min.   :    4   Length:5517        Min.   :    4  
##  1st Qu.: 2177   1st Qu.:  204   Class :character   1st Qu.:  408  
##  Median : 4703   Median : 1288   Mode  :character   Median : 1672  
##  Mean   : 7815   Mean   : 3828                      Mean   : 7248  
##  3rd Qu.:10235   3rd Qu.: 4449                      3rd Qu.: 4986  
##  Max.   :23914   Max.   :67314                      Max.   :58558  
##                                                                    
##    Locality          PopGrowthRef        Sex                 Loo         
##  Length:5517        Min.   :     9   Length:5517        Min.   :   2.10  
##  Class :character   1st Qu.:  2081   Class :character   1st Qu.:  25.90  
##  Mode  :character   Median : 13304   Mode  :character   Median :  41.00  
##                     Mean   : 38231                      Mean   :  62.46  
##                     3rd Qu.: 74354                      3rd Qu.:  77.20  
##                     Max.   :119692                      Max.   :1966.00  
##                                                                          
##    TLinfinity            K                 to            Winfinity       
##  Min.   :   2.10   Min.   : 0.0070   Min.   :-12.210   Min.   :       0  
##  1st Qu.:  26.80   1st Qu.: 0.1590   1st Qu.: -1.700   1st Qu.:     193  
##  Median :  41.97   Median : 0.2700   Median : -0.858   Median :     721  
##  Mean   :  63.92   Mean   : 0.4517   Mean   : -1.200   Mean   :   31059  
##  3rd Qu.:  79.00   3rd Qu.: 0.5300   3rd Qu.: -0.300   3rd Qu.:    4526  
##  Max.   :1966.00   Max.   :12.9570   Max.   :  1.100   Max.   :20000000  
##                                      NA's   :2883      NA's   :3091      
##      Type              Auxim              LinfLmax   LogKLogLoo     
##  Length:5517        Length:5517        Min.   :0   Min.   :-0.8340  
##  Class :character   Class :character   1st Qu.:0   1st Qu.:-0.4215  
##  Mode  :character   Mode  :character   Median :0   Median :-0.3149  
##                                        Mean   :0   Mean   :-0.2625  
##                                        3rd Qu.:0   3rd Qu.:-0.1297  
##                                        Max.   :0   Max.   : 0.9892  
##                                                    NA's   :2659     
##        tm                Lm            LmSex              TypeLm         
##  Min.   :  0.200   Min.   :  1.50   Length:5517        Length:5517       
##  1st Qu.:  2.000   1st Qu.: 16.50   Class :character   Class :character  
##  Median :  4.000   Median : 27.10   Mode  :character   Mode  :character  
##  Mean   :  6.315   Mean   : 41.77                                        
##  3rd Qu.:  7.200   3rd Qu.: 49.80                                        
##  Max.   :156.000   Max.   :922.00                                        
##  NA's   :5340      NA's   :4766                                          
##      LmMale          LmFemale      GrowthEnviron        tunas        
##  Min.   :  6.00   Min.   :  8.00   Length:5517        Mode :logical  
##  1st Qu.: 19.85   1st Qu.: 20.40   Class :character   FALSE:5309     
##  Median : 29.00   Median : 34.70   Mode  :character   TRUE :208      
##  Mean   : 48.76   Mean   : 49.44                                     
##  3rd Qu.: 50.00   3rd Qu.: 62.50                                     
##  Max.   :280.00   Max.   :257.00                                     
##  NA's   :5410     NA's   :5360                                       
##  deepseasharks    roughies                speciestype  
##  Mode :logical   Mode :logical   Deep-sea sharks:  58  
##  FALSE:5459      FALSE:5511      Other species  :5245  
##  TRUE :58        TRUE :6         Roughies       :   6  
##                                  Tunas          : 208  
##                                                        
##                                                        
## 

We then add our coelacanth data to Fishbase Von Bertalanffy growth data.

colnames(growthdata)
##  [1] "AutoCtr"       "SpecCode"      "Species"       "StockCode"    
##  [5] "Locality"      "PopGrowthRef"  "Sex"           "Loo"          
##  [9] "TLinfinity"    "K"             "to"            "Winfinity"    
## [13] "Type"          "Auxim"         "LinfLmax"      "LogKLogLoo"   
## [17] "tm"            "Lm"            "LmSex"         "TypeLm"       
## [21] "LmMale"        "LmFemale"      "GrowthEnviron" "tunas"        
## [25] "deepseasharks" "roughies"      "speciestype"
coelacanth.data<-data.frame(AutoCtr=max(growthdata$AutoCtr)+1:3,SpecCode=rep(2063,3),Species=rep("Latimeria chalumnae",3),StockCode=rep(NA,3),Locality=rep("Comoro islands",3),PopGrowthRef=rep(NA,3),Sex=rep("mixed",3),Loo=c(203.8,296.1,286.4),TLinfinity=c(203.8,296.1,286.4),K=c(0.02,0.05,0.04),to=rep(NA,3),Winfinity=rep(NA,3),Type=rep("TL",3),Auxim=rep(NA,3),LinfLmax=rep(NA,3),LogKLogLoo=rep(NA,3),tm=rep(NA,3),Lm=rep(NA,3),LmSex=rep(NA,3),TypeLm=rep(NA,3),LmMale=rep(NA,3),LmFemale=rep(NA,3),GrowthEnviron=rep("open waters",3),tunas=rep(F,3),deepseasharks=rep(F,3),roughies=rep(F,3),speciestype=c("Coelacanth circuli","Coelacanth macrocirculi","Coelacanth previous study"))
head(growthdata)
##   AutoCtr SpecCode           Species StockCode
## 1     139        4 Engraulis ringens         4
## 2     140        4 Engraulis ringens         4
## 3     141        4 Engraulis ringens         4
## 4     142        4 Engraulis ringens         4
## 5     143        4 Engraulis ringens         4
## 6     144        4 Engraulis ringens         4
##                            Locality PopGrowthRef     Sex  Loo TLinfinity    K
## 1                           Iquique            9 unsexed 16.9       16.9 1.60
## 2 Northern/central stock, 4-14<b0>S            9 unsexed 18.2       18.2 0.75
## 3 Northern/central stock, 4-14<b0>S            9 unsexed 18.5       18.5 0.87
## 4 Northern/central stock, 4-14<b0>S            9 unsexed 19.0       19.0 0.62
## 5                           Iquique            9 unsexed 19.0       19.0 1.11
## 6 Northern/central stock, 4-14<b0>S            9 unsexed 19.3       19.3 0.75
##      to Winfinity Type          Auxim LinfLmax LogKLogLoo tm Lm LmSex TypeLm
## 1 -0.02        34   TL within ellipse        0    0.16624 NA NA  <NA>   <NA>
## 2    NA        41   TL within ellipse        0   -0.09915 NA NA  <NA>   <NA>
## 3    NA        43   TL within ellipse        0   -0.04773 NA NA  <NA>   <NA>
## 4    NA        46   TL within ellipse        0   -0.16235 NA NA  <NA>   <NA>
## 5    NA        NA   TL within ellipse        0    0.03544 NA NA  <NA>   <NA>
## 6    NA        48   TL within ellipse        0   -0.09719 NA NA  <NA>   <NA>
##   LmMale LmFemale GrowthEnviron tunas deepseasharks roughies   speciestype
## 1     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 2     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 3     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 4     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 5     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 6     NA       NA   open waters FALSE         FALSE    FALSE Other species
dim(growthdata)
## [1] 5517   27
growthdata<-rbind(growthdata,coelacanth.data)
dim(growthdata)
## [1] 5520   27

We can then compute the growth performance index \(\phi^\prime\) as proposed by Pauly and Munro (1984).

growthdata$gindex<-log10(growthdata$K)+2*log10(growthdata$TLinfinity)#Compute growth index log10(K)+2log10(Linfinity) for all stocks
head(growthdata)
##   AutoCtr SpecCode           Species StockCode
## 1     139        4 Engraulis ringens         4
## 2     140        4 Engraulis ringens         4
## 3     141        4 Engraulis ringens         4
## 4     142        4 Engraulis ringens         4
## 5     143        4 Engraulis ringens         4
## 6     144        4 Engraulis ringens         4
##                            Locality PopGrowthRef     Sex  Loo TLinfinity    K
## 1                           Iquique            9 unsexed 16.9       16.9 1.60
## 2 Northern/central stock, 4-14<b0>S            9 unsexed 18.2       18.2 0.75
## 3 Northern/central stock, 4-14<b0>S            9 unsexed 18.5       18.5 0.87
## 4 Northern/central stock, 4-14<b0>S            9 unsexed 19.0       19.0 0.62
## 5                           Iquique            9 unsexed 19.0       19.0 1.11
## 6 Northern/central stock, 4-14<b0>S            9 unsexed 19.3       19.3 0.75
##      to Winfinity Type          Auxim LinfLmax LogKLogLoo tm Lm LmSex TypeLm
## 1 -0.02        34   TL within ellipse        0    0.16624 NA NA  <NA>   <NA>
## 2    NA        41   TL within ellipse        0   -0.09915 NA NA  <NA>   <NA>
## 3    NA        43   TL within ellipse        0   -0.04773 NA NA  <NA>   <NA>
## 4    NA        46   TL within ellipse        0   -0.16235 NA NA  <NA>   <NA>
## 5    NA        NA   TL within ellipse        0    0.03544 NA NA  <NA>   <NA>
## 6    NA        48   TL within ellipse        0   -0.09719 NA NA  <NA>   <NA>
##   LmMale LmFemale GrowthEnviron tunas deepseasharks roughies   speciestype
## 1     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 2     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 3     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 4     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 5     NA       NA   open waters FALSE         FALSE    FALSE Other species
## 6     NA       NA   open waters FALSE         FALSE    FALSE Other species
##     gindex
## 1 2.659893
## 2 2.395204
## 3 2.473863
## 4 2.349899
## 5 2.602830
## 6 2.446176

We can now proceed with auximetric analysis per se. We start by producing the grid of regression lines \(log_{10}(K)\sim log_{10}(L_\infty)\) and associated labels corresponding to various intercepts \(\phi^\prime\) indicating various growth performance indices.

phitextposition<-function (phi,Klim,TLlim){
  phi<-phi+0.15
  Klim<-log10(Klim)
  TLlim<-log10(TLlim)
  (xup<-(phi-Klim[2])/2)
  if (TLlim[1]<xup & xup<TLlim[2]) {
    (positionLU<-10^c(xup,Klim[2]))
  } else if (xup<TLlim[1]) {
    yleft<-phi-2*TLlim[1]
    (positionLU<-10^c(TLlim[1],yleft))
  }
  
  (xdown<-(phi-Klim[1])/2)
  if (TLlim[1]<xdown & xdown<TLlim[2]) {
    (positionRD<-10^c(xdown,Klim[1]))
  } else if (xdown>TLlim[2]) {
    yright<-phi-2*TLlim[2]
    (positionRD<-10^c(TLlim[2],yright))
  }
  return(list(posLU=positionLU,posRD=positionRD))
}

phiindex<-0:7

Klim<-c(0.005,20)
TLlim<-c(2,2000)
gindexlabels<-lapply(phiindex,phitextposition,Klim=Klim,TLlim=TLlim)
gindexLU_xpos<-sapply(lapply(gindexlabels,'[[',1),'[[',1)
gindexLU_ypos<-sapply(lapply(gindexlabels,'[[',1),'[[',2)
gindexRD_xpos<-sapply(lapply(gindexlabels,'[[',2),'[[',1)
gindexRD_ypos<-sapply(lapply(gindexlabels,'[[',2),'[[',2)

We then produce the auximetric plot using packages ggplot2 for producing plots, scales for scaling plot axes and extrafont for importing system fonts.

library(ggplot2)
library(scales)
library(extrafont)
## Warning: package 'extrafont' was built under R version 4.0.3
loadfonts() #loading system fonts

s<-8 #setting size of points representing coelacanth data
p1<-ggplot(data=subset(growthdata,speciestype="Other species"), aes(x = TLinfinity, y = K)) + 
  geom_abline(intercept=phiindex,slope=-2,linetype="dashed") + #Growth performance index Phi isolines
  geom_point(aes(fill="Other species",shape="Other species"),color="white",size=6)+ #Data points for other species
  geom_point(data=subset(growthdata, tunas==T),aes(fill="Tunas",shape="Tunas"),color="white",size=6)+ #Data points for Tunas
  stat_ellipse(data=subset(growthdata, tunas==T),type="norm",na.rm=T,color="#CD9600",size=1.5)+ #Ellipsis for Tunas
  geom_point(data=subset(growthdata, deepseasharks==T),aes(fill="Deep-sea sharks",shape="Deep-sea sharks"),color="white",size=6)+ #Data points for Deep-sea sharks
  stat_ellipse(data=subset(growthdata, deepseasharks==T),type="norm",na.rm=T,color="#C77CFF",size=1.5)+ #Ellipsis for Deep-sea sharks
  geom_point(data=subset(growthdata, roughies==T),aes(fill="Roughies",shape="Roughies"),color="white",size=6)+ #Data points for Roughies
  stat_ellipse(data=subset(growthdata, roughies==T),type="norm",na.rm=T,color="#7CAE00",size=1.5)+ #Ellipsis for Roughies
  geom_point(data=subset(growthdata, speciestype=="Coelacanth circuli"),aes(fill="Coelacanth circuli",shape="Coelacanth circuli"),size=s,color="white")+ #Data point for coelacanth using circuli for ageing; present study
  geom_text(data=subset(growthdata, speciestype=="Coelacanth circuli"),aes(label="Coelacanth circuli"),hjust=-0.1,size=5)+
  geom_point(data=subset(growthdata, speciestype=="Coelacanth macrocirculi"),aes(fill="Coelacanth macrocirculi",shape="Coelacanth macrocirculi"),size=s,color="white")+ #Data point for coelacanth using macro-circuli for ageing; present study
  geom_text(data=subset(growthdata, speciestype=="Coelacanth macrocirculi"),aes(label="Coelacanth macrocirculi"),hjust=-0.1,size=5)+
  geom_point(data=subset(growthdata, speciestype=="Coelacanth previous study"),aes(fill="Coelacanth previous study",shape="Coelacanth previous study"),size=s,color="white")+ #Data point for coelacanth using macro-circuli for ageing; previous study by Froese & Palomares
  geom_text(data=subset(growthdata, speciestype=="Coelacanth previous study"),aes(label="Coelacanth previous study"),hjust=-0.1,size=5)+
  geom_text(data=data.frame(phiindex,phiindex,label=paste("phi*minute == ",phiindex)),aes(label=label, x=gindexLU_xpos ,y=gindexLU_ypos),vjust=2,hjust=-0.5,parse=T,size=5)+
  scale_shape_manual(name="Taxonomic\ngroup",values=c(
    "Coelacanth circuli"=23,
    "Coelacanth macrocirculi"=23,
    "Coelacanth previous study"=23,
    "Deep-sea sharks"=21,
    "Other species"=21,
    "Roughies"=21,
    "Tunas"=21))+
  scale_fill_manual(name="Taxonomic\ngroup",values=c(
    "Coelacanth circuli"="#F8766D",
    "Coelacanth macrocirculi"="#00BFC4",
    "Coelacanth previous study"="#619CFF",
    "Deep-sea sharks"="#C77CFF",
    "Other species"="darkgrey",
    "Roughies"="#7CAE00",
    "Tunas"="#CD9600"))+
  theme(legend.position="left",
        legend.key=element_blank(),
        legend.title=element_text(size=14),
        legend.text=element_text(size=12),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text=element_text(size=12),
        text=element_text(family="Helvetica"))+
  xlab(expression(bold("Asymptotic total length,"~TL[infinity]~(cm)))) + 
  ylab(expression(bold("Growth rate coefficient,"~K~(year^-1))))+
  scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x,n=3),
              labels = trans_format("log10", math_format(10^.x)),limits=c(2,2000),expand=c(0,0))+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),
              labels = trans_format("log10", math_format(10^.x)),limits=c(0.005,20),expand=c(0,0))+
  annotation_logticks(sides="trbl")+
  guides(shape = guide_legend(override.aes = list(size = 5)))

In addition to the auximetric plot, we produce a simple boxplot of the growth performance index \(\phi^\prime\) of the different species types (tunas, deep-sea sharks, roughies) to position coelacanth estimates along this univariate dimension.

We start by computing \(\phi^\prime\) for the different coelacanth estimates

gindex_circuli<-growthdata[growthdata$speciestype=="Coelacanth circuli",]$gindex
gindex_macrocirculi_pst<-growthdata[growthdata$speciestype=="Coelacanth macrocirculi",]$gindex
gindex_macrocirculi_FroesePalomares<-growthdata[growthdata$speciestype=="Coelacanth previous study",]$gindex
gindex_coelacanth<-c(gindex_circuli,gindex_macrocirculi_pst,gindex_macrocirculi_FroesePalomares)

quantilelabels<-sapply(growthdata[growthdata$speciestype %in% c("Coelacanth circuli","Coelacanth macrocirculi","Coelacanth previous study"),]$gindex,function (g) round(100*ecdf(growthdata$gindex)(g),0))

Then we produce the boxplot

p2<-ggplot(data=subset(growthdata,!(speciestype %in% c("Coelacanth circuli","Coelacanth macrocirculi","Coelacanth previous study"))), aes(x=speciestype, y=gindex)) +
  geom_boxplot(aes(fill=speciestype),show.legend=F,outlier.size=3)+ #Boxplot per species type
  stat_summary(fun=mean, geom="point", shape=21, fill="black",size=4) +
  stat_summary(data=data.frame(speciestype="Other species",gindex=gindex_circuli),fun=mean, geom="point", shape=23,fill="#F8766D",color="white", size=6) + #Adding data poitn for coelacanth aged using circuli; present study
  geom_hline(aes(yintercept=gindex_circuli),color="#F8766D") +
  stat_summary(data=data.frame(speciestype="Other species",gindex=gindex_macrocirculi_pst),fun=mean, geom="point", shape=23,fill="#00BFC4",color="white", size=6) + #Adding data poitn for coelacanth aged using macro-circuli; present study
  geom_hline(aes(yintercept=gindex_macrocirculi_pst),color="#00BFC4") +
  stat_summary(data=data.frame(speciestype="Other species",gindex=gindex_macrocirculi_FroesePalomares),fun=mean, geom="point", shape=23,fill="#619CFF",color="white", size=6) + #Adding data poitn for coelacanth aged using macro-circuli; previous study by Froese & Palomares
  geom_hline(aes(yintercept=gindex_macrocirculi_FroesePalomares),color="#619CFF") +
  geom_text(data=data.frame(label=paste(quantilelabels,"%",sep="")),aes(label=label, x=-Inf ,y=gindex_coelacanth),vjust=-0.25,size=5)+
  scale_fill_manual(values=c(
    "Deep-sea sharks"="#C77CFF",
    "Other species"="darkgrey",
    "Roughies"="#7CAE00",
    "Tunas"="#CD9600"))+
  labs(x="Taxonomic group", y=expression(bold("Growth index,"~phi*minute==log[10](K)+2*Log[10](TL[infinity])))) +
  coord_flip()+
  theme(
    axis.title.x = element_text(size=16, face="bold"),
    axis.title.y = element_text(size=16, face="bold"),
    axis.text=element_text(size=12),
    legend.key = element_blank(),
    text=element_text(family="Helvetica")#, legend.box = "horizontal"
  )

We finally combine the two plots using package ggpubr and save them.

library(ggpubr)

p_auxi<-ggarrange(p1,p2,labels=c("A","B"),ncol=1,nrow=2,heights=c(2,1),common.legend=T,legend="left")
p_auxi

ggsave("Auximetric plot_Helvetica.pdf",p_auxi,width=45,height=27,units="cm")

Slow-fast life-history continuum analysis

Life-history data

We start importing the data. These data were extracted from Fishbase using the package rfishbase on June 2020. Only species from brackish and salted waters were included.

Von Bertalanffy growth parameters were extracted as described in the previous section.

Maturity data, i.e. average age at maturity \(t_m\) and length at maturity \(L_m\), were extracted from the Maturity table. For stocks for whicha range of ages and/or lengths at maturity were available instead of a point estimate, we computed the mid-age and/or -length and used them as estimates of \(t_m\) and/or \(L_m\). If not a range but a minimum or maximum age or length at maturity was provided, we used it directly as an estimate for \(t_m\) or \(L_m\). We then extracted stocks for which we had at least an estimate of either \(t_m\) or \(L_m\), as \(t_m\) could be estimated from \(L_m\) by inverting the Von Bertalanffy growth model when available. When \(L_m\) was not given as total length (or fork length for Scombridae), it was converted using length-length relationships from the Length-Length table when the information was available.

Gestation time or spawning frequency data were extracted from the Fecundity table cross checked with the Spawning table.

Maximum longevity was extracted from the Species table.

Data were merged in a single table containing life-history traits with time dimensionality, namely Von Bertalanffy rate coefficient K, age at maturity tm, maximum longevity LongevityWild and spawning frequency SpawningFreq, plus asymptotic total length TLinfinity for correcting for size effects (see below).

LHdata<-read.csv("Life_history_traits_Fishbase_2020_June.csv")
summary(LHdata)
##    Species             SpecCode         StockCode         Sex           
##  Length:4443        Min.   :    6.0   Min.   :   38   Length:4443       
##  Class :character   1st Qu.:   30.0   1st Qu.:   85   Class :character  
##  Mode  :character   Median :  142.0   Median :  520   Mode  :character  
##                     Mean   :  945.3   Mean   : 7445                     
##                     3rd Qu.: 1327.0   3rd Qu.: 1948                     
##                     Max.   :16324.0   Max.   :56833                     
##                                                                         
##    TLinfinity           K                 to            tunas        
##  Min.   :  2.10   Min.   : 0.0200   Min.   :-8.7800   Mode :logical  
##  1st Qu.: 32.60   1st Qu.: 0.1390   1st Qu.:-1.0120   FALSE:3946     
##  Median : 54.70   Median : 0.2400   Median :-0.0300   TRUE :497      
##  Mean   : 77.57   Mean   : 0.3475   Mean   :-0.6913                  
##  3rd Qu.: 88.70   3rd Qu.: 0.4400   3rd Qu.: 0.0000                  
##  Max.   :651.00   Max.   :12.9570   Max.   : 0.7850                  
##                                                                      
##  deepseasharks    roughies       speciestype            gindex     
##  Mode :logical   Mode :logical   Length:4443        Min.   :1.402  
##  FALSE:4435      FALSE:4430      Class :character   1st Qu.:2.547  
##  TRUE :8         TRUE :13        Mode  :character   Median :2.826  
##                                                     Mean   :2.916  
##                                                     3rd Qu.:3.225  
##                                                     Max.   :4.788  
##                                                                    
##       TLm              tm         LongevityWild     SpawningFreq    
##  Min.   :  1.0   Min.   : 0.082   Min.   :  0.16   Min.   :  0.330  
##  1st Qu.: 20.2   1st Qu.: 1.500   1st Qu.:  9.00   1st Qu.:  1.000  
##  Median : 30.6   Median : 3.000   Median : 20.00   Median :  1.000  
##  Mean   : 45.8   Mean   : 4.156   Mean   : 20.94   Mean   :  1.856  
##  3rd Qu.: 43.5   3rd Qu.: 5.012   3rd Qu.: 25.00   3rd Qu.:  1.000  
##  Max.   :427.4   Max.   :86.949   Max.   :149.00   Max.   :150.000  
##  NA's   :1432

Then we add coelacanth data.

coelacanth.data[,c("Species","SpecCode","StockCode","Sex","TLinfinity","K","to","tm","deepseasharks","roughies","speciestype")]
##               Species SpecCode StockCode   Sex TLinfinity    K to tm
## 1 Latimeria chalumnae     2063        NA mixed      203.8 0.02 NA NA
## 2 Latimeria chalumnae     2063        NA mixed      296.1 0.05 NA NA
## 3 Latimeria chalumnae     2063        NA mixed      286.4 0.04 NA NA
##   deepseasharks roughies               speciestype
## 1         FALSE    FALSE        Coelacanth circuli
## 2         FALSE    FALSE   Coelacanth macrocirculi
## 3         FALSE    FALSE Coelacanth previous study
(combdata.ceolacanth<-growthdata[growthdata$Species=="Latimeria chalumnae",c("Species","SpecCode","StockCode","Sex","TLinfinity","K","to","tunas","deepseasharks","roughies","speciestype","gindex")])
##                  Species SpecCode StockCode   Sex TLinfinity     K   to tunas
## 3741 Latimeria chalumnae     2063      2258 mixed      219.0 0.058 -3.6 FALSE
## 5518 Latimeria chalumnae     2063        NA mixed      203.8 0.020   NA FALSE
## 5519 Latimeria chalumnae     2063        NA mixed      296.1 0.050   NA FALSE
## 5520 Latimeria chalumnae     2063        NA mixed      286.4 0.040   NA FALSE
##      deepseasharks roughies               speciestype   gindex
## 3741         FALSE    FALSE             Other species 3.444316
## 5518         FALSE    FALSE        Coelacanth circuli 2.919438
## 5519         FALSE    FALSE   Coelacanth macrocirculi 3.641847
## 5520         FALSE    FALSE Coelacanth previous study 3.516006
(combdata.ceolacanth<-combdata.ceolacanth[-1,])
##                  Species SpecCode StockCode   Sex TLinfinity    K to tunas
## 5518 Latimeria chalumnae     2063        NA mixed      203.8 0.02 NA FALSE
## 5519 Latimeria chalumnae     2063        NA mixed      296.1 0.05 NA FALSE
## 5520 Latimeria chalumnae     2063        NA mixed      286.4 0.04 NA FALSE
##      deepseasharks roughies               speciestype   gindex
## 5518         FALSE    FALSE        Coelacanth circuli 2.919438
## 5519         FALSE    FALSE   Coelacanth macrocirculi 3.641847
## 5520         FALSE    FALSE Coelacanth previous study 3.516006
combdata.ceolacanth$TLm<-147 #Median maturation length of males and females
combdata.ceolacanth$TL1<-c(20.14,46.91,21.98)
combdata.ceolacanth$tm<-with(combdata.ceolacanth,1-(1/K)*log((TLinfinity-TLm)/(TLinfinity-TL1)))
combdata.ceolacanth$LongevityWild<-c(85,17,20)
combdata.ceolacanth$SpawningFreq<-c(1/5,1/1.5,1/1.5)
combdata.ceolacanth<-combdata.ceolacanth[,-which(colnames(combdata.ceolacanth)=="TL1")]
head(LHdata)
##                 Species SpecCode StockCode    Sex TLinfinity     K    to tunas
## 1    Achoerodus gouldii    14420     14219 female       68.2 0.140  0.06 FALSE
## 2     Alopias pelagicus     5891      6200  mixed      328.0 0.120  0.00 FALSE
## 3 Alopias superciliosus     2534      2730   male      385.0 0.088 -4.24 FALSE
## 4 Alopias superciliosus     2534      2730 female      422.0 0.092 -4.21 FALSE
## 5 Alopias superciliosus     2534      2730 female      422.0 0.092 -4.21 FALSE
## 6 Alopias superciliosus     2534      2730   male      385.0 0.088 -4.24 FALSE
##   deepseasharks roughies   speciestype   gindex    TLm       tm LongevityWild
## 1         FALSE    FALSE Other species 2.813697  65.30 17.00000            70
## 2         FALSE    FALSE Other species 4.110929     NA 13.20000            29
## 3         FALSE    FALSE Other species 4.115404 276.00 18.57972            20
## 4         FALSE    FALSE Other species 4.214413     NA 11.50000            20
## 5         FALSE    FALSE Other species 4.214413 374.96 28.05790            20
## 6         FALSE    FALSE Other species 4.115404     NA 10.00000            20
##   SpawningFreq
## 1            1
## 2            1
## 3            1
## 4            1
## 5            1
## 6            1
dim(LHdata)
## [1] 4443   16
LHdata<-rbind(LHdata,combdata.ceolacanth)
dim(LHdata)
## [1] 4446   16

Slow-fast continuum analysis

Data inspection

We start by looking at the data after \(log_{10}\) transformation

panel.hist <- function(x, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks; nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...)
}
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y,use="complete.obs"))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste0(prefix, txt)
  if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = cex.cor * r)
}

pairs(log10(LHdata[c("TLinfinity","K","gindex","tm","LongevityWild","SpawningFreq")]),diag.panel = panel.hist,lower.panel=panel.smooth,upper.panel=panel.cor)

head(LHdata)
##                 Species SpecCode StockCode    Sex TLinfinity     K    to tunas
## 1    Achoerodus gouldii    14420     14219 female       68.2 0.140  0.06 FALSE
## 2     Alopias pelagicus     5891      6200  mixed      328.0 0.120  0.00 FALSE
## 3 Alopias superciliosus     2534      2730   male      385.0 0.088 -4.24 FALSE
## 4 Alopias superciliosus     2534      2730 female      422.0 0.092 -4.21 FALSE
## 5 Alopias superciliosus     2534      2730 female      422.0 0.092 -4.21 FALSE
## 6 Alopias superciliosus     2534      2730   male      385.0 0.088 -4.24 FALSE
##   deepseasharks roughies   speciestype   gindex    TLm       tm LongevityWild
## 1         FALSE    FALSE Other species 2.813697  65.30 17.00000            70
## 2         FALSE    FALSE Other species 4.110929     NA 13.20000            29
## 3         FALSE    FALSE Other species 4.115404 276.00 18.57972            20
## 4         FALSE    FALSE Other species 4.214413     NA 11.50000            20
## 5         FALSE    FALSE Other species 4.214413 374.96 28.05790            20
## 6         FALSE    FALSE Other species 4.115404     NA 10.00000            20
##   SpawningFreq
## 1            1
## 2            1
## 3            1
## 4            1
## 5            1
## 6            1

We can already see that the growth index \(\phi^\prime\) is very much determined by \(TL_\infty\). There is also a strong correlation between the von Bertalanffy rate \(K\), age at maturity \(t_m\) and longevity which is quite expected in case of a slow-fast continuum. The correlation with spawning cycles is less obvious.

Analysis without accounting for phylogeny

We proceeed with the first analysis that looks at the slow-fast continuum after removing the size effect but wihtout accounting for phylogenetic relationships bewteen species. For this, we regress each log10-transformed life-history traits with time dimensionality on log10-transformed \(TL_\infty\) and take the residuals. We then apply a centred and scaled PCA on the residuals. We use packages ade4 and factoextra.

library(ade4)
library(factoextra) #load package for ggplot-like plotting of pca results
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
partial.matrix<-sapply(LHdata[c("K","tm","LongevityWild","SpawningFreq")],function (x) resid(lm(log10(x)~log10(LHdata$TLinfinity))))  #Partialling out size effect through linear regression
  
pc.cor.partial<-dudi.pca(partial.matrix, scannf = F,nf=5) #Perfoming centred and scaled PCA
  
(var.expl.partial<-round(100*pc.cor.partial$eig/sum(pc.cor.partial$eig),1))
## [1] 54.6 24.5 12.5  8.5
cumsum(var.expl.partial) #The 2 first axes explain 79.1%
## [1]  54.6  79.1  91.6 100.1

We produce the customized biplot of the PCA.

pcresults.partial<-data.frame(pc.cor.partial$li,speciestype=LHdata$speciestype) #Gathering results

s<-8 #setting size of points representing coelacanth data
  
p3<-ggplot(data=subset(pcresults.partial,speciestype=="Other species"),aes(x=Axis1,y=Axis2))+
    geom_point(aes(fill="Other species",shape="Other species"),color="white",size=6)+ #Other species stocks projection on the plane of the 2 first PCs
    geom_point(data=subset(pcresults.partial,speciestype=="Tunas"),aes(fill="Tunas",shape="Tunas"),color="white",size=6)+ #Tunas' stocks projection on the plane of the 2 first PCs
    stat_ellipse(data=subset(pcresults.partial,speciestype=="Tunas"),type="norm",na.rm=T,color="#CD9600",size=1)+ #Ellipsis for tunas
    geom_point(data=subset(pcresults.partial,speciestype=="Deep-sea sharks"),aes(fill="Deep-sea sharks",shape="Deep-sea sharks"),color="white",size=6)+ #Deep-sea sharks' stocks projection on the plane of the 2 first PCs
    stat_ellipse(data=subset(pcresults.partial,speciestype=="Deep-sea sharks"),type="norm",na.rm=T,color="#C77CFF",size=1)+ #Ellipsis for deep-sea sharks
    geom_point(data=subset(pcresults.partial,speciestype=="Roughies"),aes(fill="Roughies",shape="Roughies"),color="white",size=6)+ #Roughies' stocks projection on the plane of the 2 first PCs
    stat_ellipse(data=subset(pcresults.partial,speciestype=="Roughies"),type="norm",na.rm=T,color="#7CAE00",size=1)+ #Ellipsis for roughies
    geom_point(data=subset(pcresults.partial,speciestype=="Coelacanth circuli"),aes(fill="Coelacanth circuli",shape="Coelacanth circuli"),size=s,color="white")+ #Coelacanth projection when aged using circuli; present study
    geom_text(data=subset(pcresults.partial, speciestype=="Coelacanth circuli"),aes(label="Coelacanth circuli"),hjust=0.5,vjust=-1.4,size=5)+
    geom_point(data=subset(pcresults.partial,speciestype=="Coelacanth macrocirculi"),aes(fill="Coelacanth macrocirculi",shape="Coelacanth macrocirculi"),size=s,color="white")+ #Coelacanth projection when aged using macro-circuli; present study
    geom_text(data=subset(pcresults.partial, speciestype=="Coelacanth macrocirculi"),aes(label="Coelacanth macrocirculi"),hjust=0.1,vjust=2,size=5)+
    geom_point(data=subset(pcresults.partial,speciestype=="Coelacanth previous study"),aes(fill="Coelacanth previous study",shape="Coelacanth previous study"),size=s,color="white")+ #Coelacanth projection when aged using macro-circuli; previosu study by Froese & Palomares
    geom_text(data=subset(pcresults.partial, speciestype=="Coelacanth previous study"),aes(label="Coelacanth previous study"),hjust=0.85,vjust=2,size=5)+
    geom_segment(data=pc.cor.partial$c1,aes(x=rep(0,4),y=rep(0,4),xend=5*CS1,yend=5*CS2),arrow=arrow(length=unit(0.015, "npc")),size=1)+ #Arrows representing life-history traits
    geom_text(data=pc.cor.partial$c1,aes(x=5*CS1,y=5*CS2,label=c("Growth coefficient","Maturation age","Longevity","Spawning frequency")),nudge_x=0.75*c(1.2,-0.9,-0.6,0),nudge_y=0.2*c(0,-0.15,1,1),size=5)+
    scale_shape_manual(name="Taxonomic\ngroup",values=c(
    "Coelacanth circuli"=23,
    "Coelacanth macrocirculi"=23,
    "Coelacanth previous study"=23,
    "Deep-sea sharks"=21,
    "Other species"=21,
    "Roughies"=21,
    "Tunas"=21))+
    scale_fill_manual(name="Taxonomic\ngroup",values=c(
    "Coelacanth circuli"="#F8766D",
    "Coelacanth macrocirculi"="#00BFC4",
    "Coelacanth previous study"="#619CFF",
    "Deep-sea sharks"="#C77CFF",
    "Other species"="darkgrey",
    "Roughies"="#7CAE00",
    "Tunas"="#CD9600"))+
    xlab(paste("PC1 (",var.expl.partial[1],"%)",sep=""))+
    ylab(paste("PC2 (",var.expl.partial[2],"%)",sep=""))+
    scale_x_continuous(breaks=c(-4,-2,0,2,4))+
    scale_y_continuous(breaks=c(-2.5,0,2.5,5,7.5))+
    theme(legend.position="left",
        legend.key=element_blank(),
        legend.title=element_text(size=14),
        legend.text=element_text(size=12),
        axis.title.x = element_text(size=16, face="bold"),
        axis.title.y = element_text(size=16, face="bold"),
        axis.text=element_text(size=12),
        text=element_text(family="Helvetica"))+
    geom_hline(yintercept=0,linetype=2)+
    geom_vline(xintercept=0,linetype=2)+
    guides(shape = guide_legend(override.aes = list(size = 5)))

In addition to the PCA biplot, we produce a simple boxplot of the stock scores along the first PC taken an index of the slow-fast continuum for the different species types (tunas, deep-sea sharks, roughies) to position coelacanth estimates along this univariate dimension.

We start by extracting the scores for the different coelacanth estimates

Axis1_circuli<-pcresults.partial[pcresults.partial$speciestype=="Coelacanth circuli",]$Axis1
Axis1_macrocirculi_pst<-pcresults.partial[pcresults.partial$speciestype=="Coelacanth macrocirculi",]$Axis1
Axis1_macrocirculi_FroesePalomares<-pcresults.partial[pcresults.partial$speciestype=="Coelacanth previous study",]$Axis1
Axis1_coelacanth<-c(Axis1_circuli,Axis1_macrocirculi_pst,Axis1_macrocirculi_FroesePalomares)
  
quantilelabels.sf<-sapply(pcresults.partial[pcresults.partial$speciestype %in% c("Coelacanth circuli","Coelacanth macrocirculi","Coelacanth previous study"),]$Axis1,function (g) round(100*ecdf(pcresults.partial$Axis1)(g),0))

Then, we produce the boxplot.

p4<-ggplot(data=subset(pcresults.partial,!(speciestype %in% c("Coelacanth circuli","Coelacanth macrocirculi","Coelacanth previous study"))), aes(x=speciestype, y=Axis1)) +
  geom_boxplot(aes(fill=speciestype),show.legend=F,outlier.size=3)+ #Boxplot per species type
  stat_summary(fun=mean, geom="point", shape=21, fill="black",size=4) +
  stat_summary(data=data.frame(speciestype="Other species",Axis1=Axis1_circuli),fun=mean, geom="point", shape=23,fill="#F8766D",color="white", size=6) + #Adding coelacanth aged using circuli; present study
  geom_hline(aes(yintercept=Axis1_circuli),color="#C77CFF") +
  stat_summary(data=data.frame(speciestype="Other species",Axis1=Axis1_macrocirculi_pst),fun=mean, geom="point", shape=23,fill="#00BFC4",color="white", size=6) + #Adding coelacanth aged using macro-circuli; present study
  geom_hline(aes(yintercept=Axis1_macrocirculi_pst),color="#CD9600") +
  stat_summary(data=data.frame(speciestype="Other species",Axis1=Axis1_macrocirculi_FroesePalomares),fun=mean, geom="point", shape=23,fill="#619CFF",color="white", size=6) + #Adding coelacanth aged using macro-circuli; previous study by Froese & Palomares
  geom_hline(aes(yintercept=Axis1_macrocirculi_FroesePalomares),color="#619CFF") +
  geom_text(data=data.frame(label=paste(quantilelabels.sf,"%",sep="")),aes(label=label, x=-Inf ,y=Axis1_coelacanth),vjust=-0.25,size=5)+
  scale_fill_manual(values=c(
    "Deep-sea sharks"="#C77CFF",
    "Other species"="darkgrey",
    "Roughies"="#7CAE00",
    "Tunas"="#CD9600"))+
  labs(x="Taxonomic group", y="Slow-fast life-history axis, PC1") +
  coord_flip()+
  theme(
    axis.title.x = element_text(size=16, face="bold"),
    axis.title.y = element_text(size=16, face="bold"),
    axis.text=element_text(size=12),
    legend.key = element_blank(),
    text=element_text(family="Helvetica")
  )

We finally combine the two plots and save them.

p_pca<-ggarrange(p3,p4,labels=c("A","B"),ncol=1,nrow=2,heights=c(2,1),common.legend=T,legend="left")  
p_pca