rates of molecular evolution vary dramatically among organismal lineages and across genomes ( bromham and penny 2003 ) and understanding what causes this rate variation is a fundamental topic in evolutionary biology ( lanfear et al . 2010 ) .
past studies reporting on rate variation usually focused on establishing a correlation between substitution rates and those organismal traits that might potentially be affecting the supply of mutations ( e.g. , generation time or metabolic rate ; wu and li 1985 ; martin and palumbi 1993 ; smith and donoghue 2008 ; korall et al . 2010 ; gaut et al . 2011 ; lanfear et al .
2013 ) , or the rate at which available mutations are fixed in the population ( e.g. , weakened purifying selection due to small population sizes ; woolfit and bromham 2003 , 2005 ) . establishing such correlations is important for understanding broad evolutionary processes in relation to certain life history traits .
few studies , on the other hand , have investigated whether there are other factors that might also influence evolutionary rates .
plant chloroplast genomes ( plastomes ) generally comprise a pair of inverted repeat ( ir ) regions and two single - copy ( sc ) regions .
the two ir copies are identical in sequence but run in opposite directions ( fig .
1a ) , and their sequence identity is maintained by gene conversion ( birky and walsh 1992 ) .
biased gene conversion , whereby new mutations are preferentially corrected back to ancestral states ( birky and walsh 1992 ) , is hypothesized to be responsible for a significantly lower rate of nucleotide substitution in the ir compared to the sc region ( clegg et al .
1984 ; wolfe et al . 1987 ; wu and chaw 2015 ; zhu et al . 2016 ) .
if the ir region does indeed suppress nucleotide substitutions , it is expected that genes moving into or out of the ir ( via expansion or contraction of the region ) would experience corresponding substitution rate shifts  decelerating after entering the ir and accelerating after exiting it . 

( a ) the typical plant chloroplast genome ( plastome ) comprises a pair of inverted repeat ( ir ) regions separating a single copy region ( sc ) .
the red arrowheads indicate the direction in which the rrna genes in the ir are transcribed .
( b ) in some fern chloroplast genomes , the direction of transcription in the ir is inverted .
the phylogeny on the left shows the relationships among the sampled fern chloroplast genomes ; a part of their genome organization is shown on the right .
gene lengths are not to scale , gene arrow tips indicate the direction of transcription , and a few genes are omitted for clarity .  *  indicates that chln is not always present .
note that rps12 in ophioglossales , psilotales , and equisetales lacks the second intron , and therefore there is no exon 3 .
( a ) the typical plant chloroplast genome ( plastome ) comprises a pair of inverted repeat ( ir ) regions separating a single copy region ( sc ) .
the red arrowheads indicate the direction in which the rrna genes in the ir are transcribed .
( b ) in some fern chloroplast genomes , the direction of transcription in the ir is inverted .
the phylogeny on the left shows the relationships among the sampled fern chloroplast genomes ; a part of their genome organization is shown on the right .
gene lengths are not to scale , gene arrow tips indicate the direction of transcription , and a few genes are omitted for clarity . 
note that rps12 in ophioglossales , psilotales , and equisetales lacks the second intron , and therefore there is no exon 3 .
only a few studies have examined whether the ir region can influence gene evolutionary rates .
two of these focused on situations where genes in the ir were translocated to the sc .
perry and wolfe ( 2002 ) discovered that in legumes , where the ir region has disappeared , the former ir genes experienced accelerated substitution rates . however , lin et al .
( 2012 ) found that when ycf2 was translocated out of the ir in ginkgo biloba , its substitution rate did not increase ( although it is possible that this event was too recent for a rate difference to be detectable ) .
( 2016 ) tested whether genes translocated into the ir from the sc exhibited rate changes by examining such translocation events across a representative sample of angiosperms , gymnosperms , and ferns , and found that substitution rates did indeed decrease after gene translocation into the ir .
these three studies used pairwise synonymous substitutions to measure rate variation , which ignores the role of phylogenetic history and is unable to incorporate nucleotide substitution models to yield more accurate rate estimates . here
we estimate rate variation in a phylogenetic context , and employ a likelihood ratio test to identify significant rate shifts .
in addition , we examine whether translocation into the ir alters the selective environment experienced by plastid loci , and whether plastome structure changes cause shifts in gc content .
we focus on ferns ( the sister group to seed plants ) , where plastome rearrangements have inverted the ir orientation ( fig .
1b ) and translocated five loci  trnh , psba , ycf2 , rps7 , and exons 23 of rps12 ( henceforth simply 
this expanded and inverted ir occurs in extant schizaeales , salviniales , cyatheales , and polypodiales , whereas the ancestral ir configuration is found in gleicheniales , osmundales , marattiales , psilotales , ophioglossales , and equisetales , indicating that the ir reorganization occurred early in the evolution of leptosporangiate ferns ( fig .
independent of this major ir inversion and expansion , ir boundaries also expanded in the eusporangiate fern orders marattiales and psilotales by incorporating rps7 and rps12 into the ir ( fig .
the only other locus that transitioned into and out of the ir is ndhb ( gao et al .
2013 ) ; however , because ndhb always resides next to , or right on , the ir boundary , its inclusion or exclusion from the ir appears to be labile and difficult to resolve phylogenetically .
we compiled psba , ycf2 , rps7 , and rps12 sequence data from published plastomes ( trnh and other trna genes were not included because of their short sequence lengths ; table 1 ) , and tested the fit of two clock models : a 1-rate model and a 2-rate model .
the 1-rate model applies a single substitution rate across all branches of the phylogeny , whereas the 2-rate model allows one  in - ir  rate , when these loci are in the ir ( blue branches in fig .
2a and b ) , and another  out - ir  rate , when they are in the sc ( orange branches in fig .
, the 1-rate model was rejected in favor of the 2-rate model ( p < 0.0000001 ) , and within the 2-rate model
the  in - ir  rates were 2- to 3-fold slower than the  out - ir  rates ( fig . 2c  f ) .
this model - testing result indicates that substitution rates decreased following gene translocation into the ir , and is consistent with the findings of zhu et al .
psba intergenic spacer has reduced levels of sequence variation in polypodiales ( where it is located in the ir ) , compared to the other studied fern orders ( where it is in the sc ) . because the translocation of rps7 and rps12 into the ir occurred three times ( fig .
2b ) within ferns , we also tested a 4-rate model that allows three separate  in - ir  rates and one  out - ir  rate .
for both loci , the 4-rate model fit significantly better than the 2-rate and 1-rate models ( p < 0.0000001 ) , and all three  in - ir  rates were slower than the  out - ir  rate ( fig .
2d ) , a result that is again consistent with a rate deceleration for genes translocated into the ir . 
 fig .
( a , b ) the resultant rate assignment schemes of the 2-rate models for psba and yfc2 ( a ) and rps7 and rps12 ( b ) .
the symbols on ( b ) branches mark the rate assignment scheme for the 4-rate model ( see d ) .
( c , d ) rate differences of all sampled chloroplast genes observed for the 2-rate or 4-rate model .  *  indicates genes in which the 1-rate model was rejected in favor of the 2-rate or 4-rate model ( p < 0.0000001 ) .
the rate estimates from the 4-rate model of rps7 and rps12 are shown at the bottom of ( d ) .
the rate deceleration that is observed in ycf2 , psba , rps7 , and rps12 , is unique among chloroplast genes , and is due to their translocation into the ir . 

( a , b ) the resultant rate assignment schemes of the 2-rate models for psba and yfc2 ( a ) and rps7 and rps12 ( b ) .
the symbols on ( b ) branches mark the rate assignment scheme for the 4-rate model ( see d ) .
( c , d ) rate differences of all sampled chloroplast genes observed for the 2-rate or 4-rate model .  *  indicates genes in which the 1-rate model was rejected in favor of the 2-rate or 4-rate model ( p < 0.0000001 ) .
the rate estimates from the 4-rate model of rps7 and rps12 are shown at the bottom of ( d ) .
the rate deceleration that is observed in ycf2 , psba , rps7 , and rps12 , is unique among chloroplast genes , and is due to their translocation into the ir .
fern plastomes sampled in this study to rule out the possibility that the rate deceleration observed in psba , ycf2 , rps7 , and rps12 is part of a plastome - wide phenomenon ( e.g. , a change in polymerase proof - reading efficiency ; parkinson et al . 2005 ) , or that perhaps the improved fit is due simply to the extra parameter accommodating noise in the data rather than being related to the translocation ( lanfear 2011 ) , we examined other chloroplast genes that have not translocated into or out of the ir . based on the 1-rate model
, we first confirmed that fern ir genes have significantly lower substitution rates than the non - ir genes ( p < 10 
 ; two - tailed t - test ) , a result consistent with past studies on other plant groups ( clegg et al .
1984 ; wolfe et al . 1987 ; wu and chaw 2015 ; zhu et al . 2016 ) . next we compared the 1-rate and 2-rate models as above , but for genes that did not move into or out of the ir
the majority of these genes showed no rate shift ( i.e. , the 1-rate model could not be rejected in favor of the 2-rate model ) , and if they did , the direction of the rate change was predominantly opposite to that seen in psba , ycf2 , rps7 , and rps12 ( fig . 2c  f and 
 supplementary table s1 , supplementary material online ) .
the only exception is atpb ; although it experienced no translocation , it nonetheless showed a significant deceleration ( fig .
however , the rate change in atpb is only 1.2-fold , much lower than the rate change observed in the ir - translocated loci .
this comparison demonstrates that the rate deceleration in our focal loci is exceptional among chloroplast genes , and is best explained by their translocation into the ir .
next , we tested whether the degree of rate deceleration differs between synonymous and nonsynonymous substitutions , which would suggest a change in selection pressures upon translocation into the ir .
2015 ) , a codon - based , branch - site random effects method , to test for any change in selective strength .
we found that between the  in - ir  and  out - ir 
branches , there is no significant difference in the nonsynonymous / synonymous substitution rate ratio (  ) , suggesting that protein - level selection has little or no effect on the rate decrease .
the exception is ycf2 , in which selection relaxed along the  in - ir  branches ( p < 0.0001 ; table 2 ) .
however , because relaxation of selection would tend to accelerate nucleotide substitution rates , and we still observed ycf2 deceleration , this indicates its rate decrease is not due to selection . 

2015)locusselection intensity kp valuepsba1.020.8223ycf20.42<0.0001rps71.960.2013rps120.770.1178only ycf2 showed a significant signature of selection relaxation ( k  1 ) .
test for shifts in selective pressures using relax ( wertheim et al . 2015 ) only ycf2 showed a significant signature of selection relaxation ( k  1 ) .
another hallmark of the ir is its high gc content relative to the sc region . recently
wu and chaw ( 2015 ) showed that in cycads there are more a / t to g / c substitutions in the ir ( compared to the sc ) , and they proposed a gc - biased gene conversion mechanism for the ir .
if this model is correct , we would expect to see an increase in gc content once a genomic segment translocates into the ir , particularly at the third - codon positions and in noncoding regions , which are under reduced selective constraint .
our analysis of gc content corroborates this prediction  the third codon positions of psba , ycf2 , rps7 , and rps12 , as well as the trnh 
psba intergenic spacer , have a higher gc content when they are in the ir than when in the sc region ( fig .
interestingly , gc variation is less obvious ( or opposite , in psba ) in the first and second codon positions ( fig .
3 ) , implying that selection reduces this bias by limiting the substitutions that get fixed at these positions ( since they usually result in amino acid changes ) .
it should be noted that ycf2 is the only locus with a significant gc increase at the second codon position , a result consistent with the relaxed selection we detected for it ( table 2 ) . to investigate whether the rate deceleration found above was an artifact of the gc content increase , we reanalyzed our focal loci with the third codon position removed .
we found similar rate decreases , with no significant difference in the fold rate change ( p = 0.75 , paired two - tailed t - test ; supplementary table s2 , supplementary material online ) , suggesting that the shifts in rate and in gc content are decoupled in the ir . 

relative gc content calculated as gc content of the locus divided by the whole - genome gc content .
relative gc content calculated as gc content of the locus divided by the whole - genome gc content .
we compiled psba , ycf2 , rps7 , and rps12 sequence data from published plastomes ( trnh and other trna genes were not included because of their short sequence lengths ; table 1 ) , and tested the fit of two clock models : a 1-rate model and a 2-rate model .
the 1-rate model applies a single substitution rate across all branches of the phylogeny , whereas the 2-rate model allows one  in - ir  rate , when these loci are in the ir ( blue branches in fig .
2a and b ) , and another  out - ir  rate , when they are in the sc ( orange branches in fig .
, the 1-rate model was rejected in favor of the 2-rate model ( p < 0.0000001 ) , and within the 2-rate model
the  in - ir  rates were 2- to 3-fold slower than the  out - ir  rates ( fig . 2c  f ) .
this model - testing result indicates that substitution rates decreased following gene translocation into the ir , and is consistent with the findings of zhu et al .
psba intergenic spacer has reduced levels of sequence variation in polypodiales ( where it is located in the ir ) , compared to the other studied fern orders ( where it is in the sc ) . because the translocation of rps7 and rps12 into the ir occurred three times ( fig .
2b ) within ferns , we also tested a 4-rate model that allows three separate  in - ir  rates and one  out - ir  rate .
for both loci , the 4-rate model fit significantly better than the 2-rate and 1-rate models ( p < 0.0000001 ) , and all three  in - ir  rates were slower than the  out - ir  rate ( fig .
2d ) , a result that is again consistent with a rate deceleration for genes translocated into the ir . 
 fig .
( a , b ) the resultant rate assignment schemes of the 2-rate models for psba and yfc2 ( a ) and rps7 and rps12 ( b ) .
the symbols on ( b ) branches mark the rate assignment scheme for the 4-rate model ( see d ) .
( c , d ) rate differences of all sampled chloroplast genes observed for the 2-rate or 4-rate model .  *  indicates genes in which the 1-rate model was rejected in favor of the 2-rate or 4-rate model ( p < 0.0000001 ) .
the rate estimates from the 4-rate model of rps7 and rps12 are shown at the bottom of ( d ) .
the rate deceleration that is observed in ycf2 , psba , rps7 , and rps12 , is unique among chloroplast genes , and is due to their translocation into the ir . 

( a , b ) the resultant rate assignment schemes of the 2-rate models for psba and yfc2 ( a ) and rps7 and rps12 ( b ) .
the symbols on ( b ) branches mark the rate assignment scheme for the 4-rate model ( see d ) .
( c , d ) rate differences of all sampled chloroplast genes observed for the 2-rate or 4-rate model .  *  indicates genes in which the 1-rate model was rejected in favor of the 2-rate or 4-rate model ( p < 0.0000001 ) .
the rate estimates from the 4-rate model of rps7 and rps12 are shown at the bottom of ( d ) .
the rate deceleration that is observed in ycf2 , psba , rps7 , and rps12 , is unique among chloroplast genes , and is due to their translocation into the ir .
fern plastomes sampled in this study to rule out the possibility that the rate deceleration observed in psba , ycf2 , rps7 , and rps12 is part of a plastome - wide phenomenon ( e.g. , a change in polymerase proof - reading efficiency ; parkinson et al . 2005 ) , or that perhaps the improved fit is due simply to the extra parameter accommodating noise in the data rather than being related to the translocation ( lanfear 2011 ) , we examined other chloroplast genes that have not translocated into or out of the ir . based on the 1-rate model
, we first confirmed that fern ir genes have significantly lower substitution rates than the non - ir genes ( p < 10 
 ; two - tailed t - test ) , a result consistent with past studies on other plant groups ( clegg et al .
1984 ; wolfe et al . 1987 ; wu and chaw 2015 ; zhu et al . 2016 ) . next we compared the 1-rate and 2-rate models as above , but for genes that did not move into or out of the ir
the majority of these genes showed no rate shift ( i.e. , the 1-rate model could not be rejected in favor of the 2-rate model ) , and if they did , the direction of the rate change was predominantly opposite to that seen in psba , ycf2 , rps7 , and rps12 ( fig . 2c  f and 
 supplementary table s1 , supplementary material online ) .
the only exception is atpb ; although it experienced no translocation , it nonetheless showed a significant deceleration ( fig .
however , the rate change in atpb is only 1.2-fold , much lower than the rate change observed in the ir - translocated loci .
this comparison demonstrates that the rate deceleration in our focal loci is exceptional among chloroplast genes , and is best explained by their translocation into the ir .
next , we tested whether the degree of rate deceleration differs between synonymous and nonsynonymous substitutions , which would suggest a change in selection pressures upon translocation into the ir .
2015 ) , a codon - based , branch - site random effects method , to test for any change in selective strength .
we found that between the  in - ir  and  out - ir  branches , there is no significant difference in the nonsynonymous / synonymous substitution rate ratio (  ) , suggesting that protein - level selection has little or no effect on the rate decrease .
the exception is ycf2 , in which selection relaxed along the  in - ir  branches ( p < 0.0001 ; table 2 ) .
however , because relaxation of selection would tend to accelerate nucleotide substitution rates , and we still observed ycf2 deceleration , this indicates its rate decrease is not due to selection . 

2015)locusselection intensity kp valuepsba1.020.8223ycf20.42<0.0001rps71.960.2013rps120.770.1178only ycf2 showed a significant signature of selection relaxation ( k  1 ) .
test for shifts in selective pressures using relax ( wertheim et al . 2015 ) only ycf2 showed a significant signature of selection relaxation ( k  1 ) .
another hallmark of the ir is its high gc content relative to the sc region . recently
wu and chaw ( 2015 ) showed that in cycads there are more a / t to g / c substitutions in the ir ( compared to the sc ) , and they proposed a gc - biased gene conversion mechanism for the ir .
if this model is correct , we would expect to see an increase in gc content once a genomic segment translocates into the ir , particularly at the third - codon positions and in noncoding regions , which are under reduced selective constraint .
our analysis of gc content corroborates this prediction  the third codon positions of psba , ycf2 , rps7 , and rps12 , as well as the trnh 
psba intergenic spacer , have a higher gc content when they are in the ir than when in the sc region ( fig .
3 ) . interestingly , gc variation is less obvious ( or opposite , in psba ) in the first and second codon positions ( fig .
3 ) , implying that selection reduces this bias by limiting the substitutions that get fixed at these positions ( since they usually result in amino acid changes ) .
it should be noted that ycf2 is the only locus with a significant gc increase at the second codon position , a result consistent with the relaxed selection we detected for it ( table 2 ) . to investigate whether the rate deceleration found above was an artifact of the gc content increase , we reanalyzed our focal loci with the third codon position removed .
we found similar rate decreases , with no significant difference in the fold rate change ( p = 0.75 , paired two - tailed t - test ; supplementary table s2 , supplementary material online ) , suggesting that the shifts in rate and in gc content are decoupled in the ir . 

 *  denotes p < 0.05 in 2-tailed t - test . relative gc content calculated as gc content of the locus divided by the whole - genome gc content .
 *  denotes p < 0.05 in 2-tailed t - test . relative gc content calculated as gc content of the locus divided by the whole - genome gc content .
in this study , we demonstrated that when genes are translocated into the ir , their nucleotide substitution rates dropped significantly ( 2- to 3-fold ) , and that this deceleration is not shared with other nontranslocated chloroplast genes .
in addition to rate deceleration , gc content increases following translocation , indicating that the ir affects both substitution rates and gc content .
our finding also points to plastome rearrangements that can result in rate heterogeneity among lineages .
this has important implications for studies trying to use genomic data to correlate rate heterogeneity together with life history traits , as well as for dating evolutionary events ( e.g. , schuettpelz and pryer 2006 ; rothfels and schuettpelz 2014 ) . without the knowledge of genome structure , or modeling for possible hidden rate shifts
we sampled at least one plastome from each of the 11 fern orders ( except for hymenophyllales ; table 1 ) .
the genome data were downloaded from genbank , and individual gene sequences were extracted . because some plastomes are incomplete , we focused on the genes that are present across all the sampled plastomes ; a total of 48 loci were included .
we inferred multiple sequence gene alignments using muscle ( edgar 2004 ) followed by manual inspection and adjustment .
we used baseml ( implemented in paml ; yang 2007 ) to estimate nucleotide substitution rates under the 1-rate ( global clock ) and 2-rate models ( see fig .
the baseml runs , we used a gtr + g substitution model , with a fixed topology and divergence times derived from the 25-nuclear - locus study of rothfels et al .
the 2-rate models were tested against the 1-rate model using likelihood - ratio tests . to test for differences in selection pressure among the  in - ir  and  out - ir  branches , we used a codon - based model - testing framework implemented in relax ( wertheim et al .
the key parameter in relax is k , which controls the degree of nonsynonymous / synonymous substitution rate ratio (  ) differences between the  in - ir  and  out - ir  branches ( in = out ) . in the null model ,
k is fixed at 1 ( i.e. , no difference in selective pressures between the in - ir versus out - ir branches ) , whereas in the alternative model , k is a free parameter .
values of k > 1 make  for the in - ir branches more extremestronger positive selection or stronger purifying selectionwhereas values of k < 1 make  closer to 1 ( i.e. , relaxation of selection ) .
both null and alternative models assume three site categories of . likelihood ratio tests are then used to compare the fit of the free - k versus fixed - k models .
we calculated relative gc content for each sequence as the gc content of the focal sequence divided by the average gc content of the plastome , in order to control for genomic background variation .
we used a two - tailed t - test to investigate whether there is any significant gc difference between the in - ir and out - ir sequences .
all the alignments and the control files for running baseml are available on figshare ( dx.doi.org/10.6084/m9.figshare.3483056.v2 ; dx.doi.org/10.6084/m9.figshare.3483059.v1 ) .
supplementary tables s1 and s2 are available at genome biology and evolution online ( http://www.gbe.oxfordjournals.org/ ) .
 plant chloroplast genomes ( plastomes ) are characterized by an inverted repeat ( ir ) region and two larger single copy ( sc ) regions . 
 patterns of molecular evolution in the ir and sc regions differ , most notably by a reduced rate of nucleotide substitution in the ir compared to the sc region . 
 in addition , the organization and structure of plastomes is fluid , and rearrangements through time have repeatedly shuffled genes into and out of the ir , providing recurrent natural experiments on how chloroplast genome structure can impact rates and patterns of molecular evolution . here 
 we examine four loci ( psba , ycf2 , rps7 , and rps12 exon 23 ) that were translocated from the sc into the ir during fern evolution . 
 we use a model - based method , within a phylogenetic context , to test for substitution rate shifts . 
 all four loci show a significant , 2- to 3-fold deceleration in their substitution rate following translocation into the ir , a phenomenon not observed in any other , nontranslocated plastid genes . 
 also , we show that after translocation , the gc content of the third codon position and of the noncoding regions is significantly increased , implying that gene conversion within the ir is gc - biased . taken together 
 , our results suggest that the ir region not only reduces substitution rates , but also impacts nucleotide composition . 
 this finding highlights a potential vulnerability of correlating substitution rate heterogeneity with organismal life history traits without knowledge of the underlying genome structure . 
Introduction
Results and Discussions
Rate Deceleration Following Translocation into IR
Rate Deceleration Is Not Likely Due to Selection at the Protein Level
GC Content Increases upon Translocation into the IR
Conclusion
Materials and Methods
Supplementary Material