Macroecology of North American suckers (Catostomidae): tests of Bergmann's and Rapoport's rules

Abstract Discerning spatial macroecological patterns in freshwater fishes has broad implications for community assembly, ecosystem dynamics, management, and conservation. This study explores the potential interspecific covariation of geographic range (Rapoport's rule) and body size (Bergmann's rule) with latitude in North American sucker fishes (Cypriniformes: Catostomidae). While numerous tests of Rapoport's and Bergmann's rules are documented in the literature, comparatively few of these studies have specifically tested for these patterns, and none have incorporated information reflecting shared ancestry into analyses of North American freshwater fish through a hierarchical model. This study utilized a hierarchical modeling approach with Bayesian inference to evaluate the role that evolution has played in shaping these distributional corollaries. Rapoport's rule was supported at the tribe level but not across family and subfamily groupings. Particularly within the Catostominae subfamily, two tribes reflected strong support for Rapoport's rule while two suggested a pattern was present. Conversely, Bergmann's rule was not supported in Catostomidae. This study provides additional information regarding the pervasiveness of these “rules” by expanding inferences in freshwater fishes and specifically addressing the potential for these macroecological patterns to play a role in the distribution of the understudied group Catostomidae.


Introduction
Large-scale macroecological patterns, or "rules", provide essential information for understanding distribution (Brown 1995), providing management recommendations (Fowler et al. 2013), and aid in refining conservation efforts (Jennings and Blanchard 2004) for populations, species, and higher order taxonomic groups. The covariation of geographic range (Rapoport's rule, Rapoport 1975;Stevens 1989) and body size (Bergmann's rule, Bergmann 1847) with latitude are among the most wellstudied macroecological patterns. These patterns have been explored in both terrestrial and aquatic systems at different taxonomic scales (e.g., intraspecific, interspecific); however, results have been mixed (reviewed in Gaston et al. 1998;Blackburn et al. 1999).
Collectively, Rapoport's and Bergmann's rules have been the subject of much debate, primarily resulting from a lack of any consistent mechanism to explain their occurrences. Explanations for Rapoport's rule include latitudinal correlations with climate variation, geologic history (e.g., glaciation), watershed area, species richness trends (e.g., competition), and species nichegeographic relationships (Gaston et al. 1998;Arita et al. 2005). Explanations for Bergmann's rule primarily invoke temperature clines concurrent with latitude that coincide with development and maturation times (Bergmann 1847;Ray 1960;Sibly and Atkinson 1994). Irrespective of mechanism, however, these "rules" still serve as useful abstractions to better understand large-scale distribution patterns.
Rapoport's rule has been documented across all North American freshwater fishes (Stevens 1989;Rohde et al. 1993). Both Stevens (1989) and Rohde et al. (1993) used geographic range data from over 700 species' (Lee et al. 1980) and identified increasing range sizes concurrent with northern latitudes. Further interpretation of these studies indicates that this pattern seems to be relegated to the Nearctic and Palearctic zoogeographic regions (i.e., above 35-40 degrees). This conclusion provides strong evidence that the rule may be a by-product of the Pleistocene glacial history of these regions. Specific to Bergmann's rule, while this hypothesis was developed in the context of interspecific body size variation, the application has been primarily in studies of intraspecific variation (Rensch 1938;reviewed in Blackburn et al. 1999). Despite the breadth of literature on the topic, comparatively few of these studies have tested for Bergmann's rule in fishes, particularly in North American freshwater fish (Belk and Houston 2002;Rypel 2014), and fewer still have explored interspecific variation in North American fishes (Knouft 2004). Belk and Houston (2002) used a dataset including length at age and maximum length data from 18 species representing 10 families. Their results did not indicate any uniform relationship between maximum length and latitude (although several species exhibited inverse relationships at particular age lengths). More recently, Rypel (2014) tested maximum lengths obtained from record angling records of 29 species representing 14 families and found results contrary to Belk and Houston (2002). Consistent with thermal niche, Rypel (2014) found that certain taxa demonstrated Bergmann's rule while others either exhibited inverse relationships or no body size trend with latitude. Specific to Catostomidae, Knouft (2004) parsed out significant positive family level relationships between latitudinal variation and mean regional community body size distributions using least squares linear regression in an analysis of North American freshwater fish.
The use of these types of comprehensive datasets provides overarching evidence for all North American freshwater fishes; however, the large taxonomic scales of these analyses also creates the potential problem of signal loss in a particular family or group that diverges from the overall pattern. For example, whole assemblage tests of Rapoport's rule have the potential to obscure patterns in particular genera or families and intraspecific tests of Bergmann's rule do not address variation between individuals or within higher clades. The relationship between range size and body size as a function of dispersal potential may also generate spurious patterns related to latitude, particularly in the recently glaciated Nearctic and Palearctic zoogeographic regions (Blackburn et al. 1999). Furthermore, few tests of Rapoport's or Bergmann's rules account for phylogeny (none in freshwater fish studies), which results in taxonomic independence issues that have the potential to also change signal or lead to invalid conclusions entirely (Clauss et al. 2013).
The taxonomic richness and phylogenetic resolution in the freshwater fish family Catostomidae (Suckers), coupled with the variation in body size and geographic range size, provides a unique case study opportunity to assess these two long standing tenets of macroecology, Bergmann's and Rapoport's rules, in an understudied group of fishes. Collectively, Catostomidae includes over 70 recognized species that occupy important niches in both lentic and lotic aquatic food webs across North America. Functionally, Catostomidae utilize their modified fleshy lips with protrusible mouth, pharyngeal arches, teeth, and pads to feed on benthic algae and invertebrates including aquatic insect larvae and mollusks (Boschung Jr. & Mayden 2004). Their importance as a basal consumer is compounded in aquatic ecosystems as a result of their abundance, size distribution, life-history patterns, and geographic distribution where in many aquatic systems Catostomidae comprise more biomass than any other group of fishes (Becker 1983), occupy a wide range of size classes (Page and Burr 2011), and exhibit the capability to link extensive reaches within systems or between streams, lakes, and rivers via extensive spring spawning migration runs (Cooke et al. 2005;Reid 2006). Traditionally, these taxa have seen little management focus; however, their roles in aquatic ecosystems have generated recent conservation interest, particularly in efforts to better understand their ecology and evolution (Cooke et al. 2005). The objective of this study was to test for the covariation of geographic range (Rapoport's rule) and body size (Bergmann's rule) with latitude in the North American freshwater fish family Catostomidae at multiple taxonomic scales to better understand these fishes and extend our understanding of the prevalence of these general ecological tenets.

Methods
Catostomidae is comprised of 72 recognized species arranged in four subfamilies and several tribes; Myxocyprininae -1 species, Ictiobinae -8 species, Cycleptinae -2 species, and Catostominae -61 species (Nelson 2006;76 species cited in Harris et al. 2014) that range in body size (TL) from about 16 cm (Roanoke Hogsucker Hypentelium roanokense) to 100 cm (Bigmouth Buffalo Ictiobus cyprinellus) and are distributed across North America occupying a wide variety of habitats (Lee et al. 1980;Page and Burr 2011). Catostominae has been further subdivided into 4 tribes: Catostomini, Thoburnini, Moxostomatini, and Erimyzonini (Doosey et al. 2010). This study used species traits (latitude, maximum body size, and areal geographic range size) compiled for 62 Catostomidae taxa from Page and Burr (2011). Taxa were selected based on data availability and to ensure taxonomic coverage of the family. Latitude was assigned using the midpoint method (Rohde et al. 1993) wherein each species' latitude was treated as an individual point rather than a band (Stevens 1989). The midpoint method was used to specifically denote latitudinal variation instead of band methods to reduce nonindependent variation in mean range size at a given latitude. However, despite these two methodological differences, these two methods most frequently result in identical conclusions (Gaston et al. 1998). All geographic information was extracted from GIS occurrence maps arranged in Page and Burr (2011) using Quantum GIS 2.0.1-Dufour (QGIS Development Team 2009). Geographic centroid (latitudinal and longitudinal in decimal degrees) was determined using the polygon centroid tool in Quantum GIS. Body size and range size were standardized to z-score.

Statistical analysis
Latitudinal midpoint (lat i ) of species i was modeled as a linear function of areal geographic range size and maximum body size. Here, lat i is modeled as a normal distribution where the mean is a linear function of areal range size (rs i ) and maximum body size (bs i ) for species i.
where l i is the mean latitudinal centroid of each species from the normal distribution (norm), a jk is the intercept and represents the hypothetical mean latitudinal centroid with a areal range size and body size of zero for subfamily j and tribe k, and b 1jk and b 2jk are model coefficients representing the effect of areal range size and body size for subfamily j and tribe k, thus representing the tribe level coefficients. To estimate the effect at different levels of species classification, subfamily and tribe (as delineated in the phylogenetic analysis of Catostomidae of Doosey et al. (2010)) were treated as random effects with tribe nested within subfamily for the intercept and effect of body size and areal range size. Thus, a jk , b 1jk and b 2jk are given a hierarchical prior: where l ak , l 1k , and l 2k represent the subfamily level intercept and effect of areal range size and body size; and r 2 aj ; r 2 1j , and r 2 2j represent the subfamily variance for the effect of areal range size and body size. The next level of the model specified global level coefficients, h: l ak $ normalðh a ; r 2 a Þ l 1k $ normalðh 1 ; r 2 1 Þ l 2k $ normalðh 2 ; r 2 2 Þ where h a , h 1 , and h 2 represent the global intercept and effect of areal range size and body size; and r 2 a ; r 2 1 , and r 2 2 represent the overall standard deviation for the effect of areal range size and body size. As areal range size and body size are known to be correlated (Gaston and Blackburn 1996), we used a Bayesian Lasso approach to include both variables in the model. The Bayesian Lasso is a variable selection technique that uses a double-exponential prior on the coefficients (Tibshirani 1996;Park and Casella 2008). The Bayesian Lasso will pull the weakest parameter to 0 thus providing a variable selection method with correlated predictors. Thus, the hyperpriors, l 1 and l 2 , were given a double-exponential prior: Further, lambda and mu.tau were given noninformative gamma priors. Uncertainty due to natural individual variation from phylogenetic relationships was accounted for in our analysis by treating phylogenetic classification (e.g., subfamily and tribe) as a random effect. This method makes it possible to directly test relationships at multiple phylogenetic classification scales. While other methods of accounting for phylogenetic uncertainty exist (e.g., de Villemereuil et al. 2012; Jacquemin and Doll 2014) they preclude the ability to assess relationships at multiple scales. For example, de Villemereuil et al. (2012) describe a method of using information from a phylogenetic tree as a variance-covariance matrix in a multivariate normal model. While this method directly incorporates the correlation of traits with closely related species, it does not allow detection of a relationship between latitudinal centroid with areal range size and body size at multiple classification scales. Further, phylogenetic classification could not be used as a random effect and phylogenetic tree information as a variance-covariance matrix in the same model because it would be using similar information multiple times, potentially biasing parameter estimates. Nevertheless, we attempted to fit a model without random effects for subfamily and tribe following the methods of de Villemereuil et al. (2012) and Jacquemin and Doll (2014) to determine an overall effect and compared the two methods using the penalized deviance information criterion (Plummer 2008). The modeling approach using phylogenetic classification as random effects was found to be the best model. Bayesian inference was used to estimate parameters of the model. We used vague (i.e., noninformative) priors for all model parameters except the correlation between slopes to indicate we presume no strong a priori knowledge of the model parameters. Independent univariate normal distributions with a mean of 0 and precision of 0.0001 were used for the individual components of h and a noninformative gamma prior with shape and scale parameter set to 0.001 was used for individual r 2 , lambda and mu.tau. To generate posterior distributions, we used JAGS 3.4 (Plummer 2003) implemented in R 3.1.3 (R Development Core Team 2015) using the rjags package (Plummer 2013). We ran 3 MCMC chains for a total of 3,850,000 steps, saving every 15 steps, and discarding the first 100,000 steps as a burn-in period, resulting in 250,000 saved steps. The burnin period is necessary to reduce the effect of the starting values on the MCMC results (Gelman et al. 2004). Convergence of the MCMC algorithm was assessed using the Brooks-Gelman-Rubin (BGR) scale-reduction factor (Brooks & Gelman 1998). The BGR factor is the ratio of between chain variability to within chain variability. Convergence is obtained when the upper limit of the BGR factor is below 1.10, indicating there is not more variability between chains compared to within chains. JAGS code to implement the model is located in the appendix.
The global coefficients for the effect of geographic range size and body size (h 1 and h 2 ) were positive; however, these did not credibly differ (95% CI) from zero, suggesting no relationship at the family level. The median estimate for the effect of areal range size was 0.033 (95% Credible Intervals = À0.525 to 4.292) and body size was 0.006 (95% Credible Intervals = À0.847 to 2.003).
Interestingly, subfamily level coefficients for the effect of geographic range size (l 1k ) were not consistent across subfamilies (Fig. 2). All three subfamilies resulted in 95% credible intervals that overlapped zero (Fig. 2). However, the subfamily Catostominae resulted in 90% credible intervals (0.014-4.234) that did not overlap zero, suggesting a positive effect.
Tribe level coefficients for the effect of geographic range size (b 1jk ) were not consistent across tribes of the subfamily Catostominae (Fig. 3). Two tribes, Catostomini and Moxostomatini, resulted in 95% credible intervals that were positive and did not overlap zero suggesting a significant positive effect (Fig. 3). However, the remaining tribes were positively skewed, suggesting a weak but positive relationship between geographic range size and latitudinal centroid (Fig. 3). Tribe level coefficients for the remaining subfamilies are not shown due to only one subfamily being present. Thus, the posterior of these tribes were similar to their subfamily.
Posterior predicted values for latitudinal centroid for the Catostomini tribe consistently increased with geographic range size (Fig. 4). There is a predicted 16% increase in the median latitudinal centroid as areal range size increased from one standard deviation below average to one standard deviation above average. This change is equivalent to a geographic distance of 657 km.
Posterior predicted values for latitudinal centroid for the Moxostomatini tribe consistently increased with geographic range size (Fig. 5). There is a predicted 18.4% increase in the median latitudinal centroid as areal range size increased from one standard deviation below average to one standard deviation above average. This change is equivalent to a geographic distance of 542 km.
Subfamily level coefficients for the effect of body size (b 2k ) were consistent across subfamilies (Fig. 6). The posterior distribution is peaked over zero, which is similar to the double-exponential prior we specified, suggesting no credible effect of body size across tribes.
Tribe level coefficients for the effect of body size (b 2k ) were consistent across tribes of the subfamily Catostominae (Fig. 7). The tribe level effects mimic those of the subfamily and were peaked at zero. Tribe levels have not been defined for the remaining subfamilies (Doosey et al. 2010).

Discussion
This study indicated corollaries in range size consistent with Rapoport's rule for the Catostomidae family. At a finer scale, the strongest corollaries occurred in tribes arranged in the Catostominae subfamily. However, no subfamily or tribe of Catostomidae supported Bergmann's rule. The lack of support for Bergmann's rule also precludes an overall interaction between body size and range size, which indicates that there is not a cumulative effect whereby larger fish are not expected to exhibit even larger range sizes with increasing latitude. The present study increases our knowledge on an understudied yet functionally important group representing a large portion of the North American freshwater fish assemblage (~8% of ichthyofauna; Harris et al. 2014).

Evolution
Catostomidae occupy one of the largest geographic distributions among freshwater fish families globally. The family exhibits a disjunct contemporary and paleo distribution between North America and Asia. This distribution pattern extends from the Yangtze River Basin to Siberia and throughout North America (Berra 2007). The most widely accepted hypothesis for the evolutionary divergence and dispersion of the Catostomidae is from Darlington (1957), who hypothesized that the group originated in Asia (Eocene epoch 35-55 mya) and radiated across North America via Beringia (and in one case, Catostomus catostomus, moved back into Siberia; Bachevskaya et al. 2014). Despite only preliminary fossil evidence when formulated, the vicariancedispersal hypothesis of Darlington (1957) has garnered recent support from expanded fossil (Cavender 1986;Chang et al. 2001) and molecular (Bachevskaya et al. 2014) records. Given the evolutionary history of Catostomidae and the role that range expansion and distribution have played in their diversification, the present study provides specific evidence as to the importance of geographic location in understanding range size variation, irrespective of body size. Interestingly, while native Catostomidae are all but extirpated from Asia (except for Myxocyprinus), they have flourished in North America. This may be the result of increased competition with Cyprinidae in Asia and the timely availability of open niches in North America (Chang et al. 2001), particularly those in smaller stream systems. Knouft and Page (2003; using a phylogenetically based analysis) and Smith (1992; using a qualitative approach) suggested that the majority of speciation events in Catostomidae have occurred as a result of smaller bodied individuals involved in smaller stream vicariance events. This coincides with the evolutionary trend of body size and habitat preference (stream size) found in the fossil record whereby deeper bodied taxa that occupy large bodies of water tend to be evolutionarily basal to more recent taxa exhibiting increasingly fusiform body shapes and occupying smaller streams (e.g., Ictiobus vs. Catostomus). From an ecological perspective, larger bodied Catostomidae have also been shown to occupy larger ranges (Pyron 1999). The results of this study, however, indicate that there is not a relationship between established range/body size corollaries and geographic position whereby smaller or larger taxa do not tend to occur further north than opposite ends of the size spectrum, irrespective of evolutionary history. The lack of any relationship with body size is surprising given the vicariance hypotheses of Smith (1992) and the increased diversity in smaller streams in the American Southeast.

Ecology
Recent macroecology literature (Knouft 2004;Griffiths 2010) has summarized several trends that tend to emerge    for all North American fishes when observed as a whole. For example, species richness tends to decline with increasing latitude concurrent with an increased proportion of larger body-sized individuals that also tend to exhibit larger geographic range sizes. However, these patterns seem to be a likely artifact of increasingly large, mobile, migratory, and generalist species acting in a colonizing fashion following Pleistocene glacial events (Knouft 2004;Griffiths 2010). Related to Catostomidae, the lack of Bergmann's rule seems to refute the body size component and visual estimation of published niche breadth data (see Pyron 1999) does not seem to suggest a relationship with latitude. Latitudinal macroecological analyses incorporating migration information, body size, and niche are necessary to formally test this multifaceted hypothesis. However, previous analyses (Pyron 1999) have indicated that Catostomidae with larger geographic range sizes do tend to exhibit higher local abundances, occupy wider ecological niches, and have larger body sizes, after accounting for phylogeny. The use of phylogenetic information in analyses of Rapoport's and Bergmann's rules in recent studies (Cruz et al. 2005;Clauss et al. 2013) represents an important step in the understanding of spatial distribution patterns. Coupling comparative methods with large-scale distribution and life-history information may ultimately help to parse out the potential contributions of ecology vs. phylogeny in shaping understanding of species distribution. Cruz et al. (2005) demonstrated improved detectability of macroecological trends such as Bergmann's rule at lower taxonomic scales (e.g., genera compared with family) and suggested that decreasing scale could better elicit specific underlying mechanisms. This conclusion is supported by Clauss et al. (2013), who identified Bergmann's rule in phylogenetic analyses but not in conventional statistics, particularly among closely related species. While our results indicated a similar trend at the family level and lower order tribe level groupings, a stronger effect was identified at the tribal level, suggesting that while Catostomidae respond similarly with respect to these macroecological patterns, there are taxonomic differences in relative effect.

Conclusion
Ultimately, the implications of identifying macroecological patterns are relevant for further disentangling evolutionary trends, community assembly ecology, and improving conservation efforts for populations, species, and higher order taxonomic groups. Due to their high biomass, variable life history, and relative abundance in aquatic ecosystems, Catostomidae serve as important functional components and indicators of ecological integrity (Harris et al. 2014). However, as their status does not relate to game fisheries, their study has not historically been emphasized to the degree of some other stocks. This study provides insight into their distribution patterns while outlining a potential template that could be applied to other taxonomic scales and groups.