Evolution of the intermuscular bones in the Cyprinidae (Pisces) from a phylogenetic perspective

Abstract Intermuscular bones (IBs) are widely present in morphologically generalized teleost fishes and are commonly found in the Cyprinidae. Intermuscular bones are small, hard spicules of bone that are formed by ossification in the myosepta between neighboring myomeres. Why fish have IBs, and whether there is any evolutionary pattern to their occurrence, has been poorly understood. However, the presence of IBs does substantially affect the meat quality and commercial values of many cyprinid fishes in aquaculture. In this study, we sampled 592 individuals of cyprinid fishes to systematically investigate the evolution of IBs from a phylogenetic point of view. We found that the total number of IBs in the Cyprinidae ranged from 73 to 169, and we clarified that only two categories of IBs (epineural and epipleural) were present in all examined cyprinids. Most of the IBs were distributed in the posterior region of the fish, which might be an optimal target for selecting fewer IB strains in aquaculture. There was a positive correlation between IBs and the number of vertebrae, thus making it possible to predict the approximate number of IBs by counting the number of vertebrae. Although the IBs displayed some correlation with phylogenetic relationships in some lineages and to ecological factors such as diet (especially carnivore), in an overall view the variations of IBs in cyprinids were extremely diverse. The number and patterns of IBs in these fishes may reflect their phylogenetic history, but have been shaped by multiple environment factors. In this study, we also confirmed that X‐ray photography remains an optimal and reliable method for the study of IBs.

. The habitats of cyprinids are very diverse, covering most freshwater areas, including extreme regions such as some of the world's highest place (e.g., the Qinghai-Tibet plateau), karst caves, and the saline conditions in coastal waters. The members of this family not only exhibit diverse features of evolutionary interest, but also have considerable applied importance to mankind (Winfield & Nelson, 1991). Representative examples include the zebrafish (Danio rerio), a model organism used in genomic and developmental biology, important aquaculture species like the common carp (Cyprinus carpio), and many popular aquarium species like goldfish (Carassius auratus; Stout et al., 2016).
One feature of the Cyprinidae, the intermuscular bones (IBs), has received very little attention from researchers. Intermuscular bones are small bony spicules that are formed by ossifications within the myosepta among neighboring myomeres (Meng, Su, & Li, 1987). They are thus a part of the axial skeleton of the endoskeleton system. Intermuscular bones, which ossify segmentally and are serially homologous in the myosepta, occur only in teleosts among recent vertebrates (Patterson & Johnson, 1995).
Almost all species of cyprinids have IBs. According to the position of attachment in fish muscle, IBs have been subdivided into three categories, namely epineural, epicentral, and epipleural (Owen, 1866). Epineurals are those IBs attached proximally to the neural arches above the horizontal septum, whereas epipleurals are attached to the hemal arches or the ribs that lie below the horizontal septum. Epicentrals attach to the central vertebrae in the plane of the horizontal septum (Owen, 1866). The IBs are of variable occurrence; for instance, in some teleosts, only epineurals and epipleurals are present, and in some Perciformes, IBs have completely disappeared (Patterson & Johnson, 1995). According to the nature of the branch patterns, IBs have been divided into seven types as follows: nonforked, one-end-unequal-bi-forked, oneend-equal-bi-forked, one-end-multi-forked, two-end-bi-forked, two-end-multi-forked, and tree-branch forking (Lv, Bao, Jiang, Yang, & Li, 2007).
Techniques for studying IBs are mainly from anatomy (dissection), clearing and staining, and radiography. All these methods have both advantages and disadvantages and are applicable for different purposes. Dissection has been the most commonly used method for studying IBs in the past few years, as it is an easy and straightforward method of getting results (Andria-Mananjara, Rasamoelina, & Vandeputte, 2016;Dong et al., 2006;Patterson & Johnson, 1995).
However, the tiny IBs in some small individuals can be easily missed.
Staining the fish bones using Alizarin red is an alternative method to observe IBs. It is suitable for most of the individuals including small ones as it can retain as much information as possible (Yao, Lv, Gong, Wu, & Bao, 2015); however, the complicated procedure of staining fish bones requires considerable times and is awkward for large individuals. In contrast to the above two more traditional methods, radiographing (usually X-ray) of specimens to observe IBs is a relatively new method which has a lot of advantages, for example, it is easy to operate and time-saving, and there is no damage to the specimen; it is especially applicable to living fishes (Perazza et al., 2017;Vallod & Arthaud, 2009;Xu, Zheng, Qian, & Luo, 2015). One of the supposed drawbacks for radiographing is the presence of overlapping areas which could result in misreading the number of IBs when taking radiograph from a single angle; however, this effect has not been evaluated as yet.
There are not very many studies involving IBs, and most of these involve only a few aspects. Many of the papers simply count and describe the IBs based on a few different species (e.g., Hypophthalmichthys molitrix, H. nobilis, and Megalobrama amblycephala in Dong et al., 2006). Although seven types of IBs based on shapes have been identified, all six of the more complex shapes were thought to be evolved from the simple nonforked shape.
Intermuscular bones found in the more anterior parts of the body were usually more complex than those present in the posterior section of the fish (Li et al., 2013). Some studies reported that the number of IBs did not increase with growth as a function standard (SL) or total (TL) length and that the number of IBs between the left and right sides of the body fish was not equal in certain individuals (Cao et al., 2014;Li et al., 2013). A recent study pointed out that IBs were probably related to the ploidy of a species; a fish with an increased ploidy would have an increased number of IBs (Li et al., 2013). An interesting pattern concerning the ossification process of IBs has also been reported: The ossification of IBs of most fishes progresses from the posterior to the anterior regions of the body (e.g., Danio rerio, Yao et al., 2015, H. molitri, Ke, Zhang, Jiang, & Bao, 2009and M. amblycephala, Wan et al., 2014; however, in Anguilla japonica the opposite situation is found (Yao et al., 2015).
More recent studies focused on the genetic basis of IBs through molecular methods. For instance, Wan et al. (2015) compared the miRNA profiles in M. amblycephala and found 20 miRNAs upregulated in IBs relative to connective tissues. A further analysis integrated with the miRNA and mRNA expression profiles has shown that many miRNA-mRNA interaction pairs might regulate the bone development and differentiation involved in IB formation (Wan et al., 2016).
The real function of IBs is not understood, although several hypotheses have been proposed (e.g., bilateral support for the muscles; Deng, 1959, transmitting muscle strength between sarcomeres; Dong et al., 2006;Lv et al., 2007, or enhancing the strength of herbivorous fishes to some extent; Danos & Staab, 2010). Although the functions of IBs are still largely unknown, the presence of IBs does substantially affect the flesh quality and the commercial values of many aquaculture fishes. Recent studies have shown that IB-less strains could be obtained by artificial breeding (Perazza et al., 2017;Xu et al., 2015), which shed light on the genetic improvement programs of selecting IB-less or even IB-null fishes.
As there are so many questions about IBs awaiting investigations, this study expands our understanding of IBs using representative cyprinid fishes. The main purposes of this study were as follows: (a) to evaluate the reliability of radiographs in counting the number and shape of IBs; (b) to study the number of IBs from a phylogenetic perspective; and (c) to find whether there are any relationships between IBs and potential ecological factors.

| Sample identifications and selections
A total of 592 specimens of cyprinid fishes were examined in the present study (Table 1). Among them, 88 individuals were fresh samples collected from Yunnan province of China, and the others were ethanol-preserved samples deposited in the Kunming Institute of Zoology (KIZ), Chinese Academy of Sciences. We identified samples according to the 12-subfamily classification system of Chen (1998), which is one of the most useful and popular systems for taxonomic studies involving Chinese cyprinids (Jiang et al., 2018).
We selected samples for a large coverage of all the 12 subfamilies, and we ensured that each subfamily included at least one genus and two species. For each diagram, we counted the total number of IBs, the epineurals, the epipleurals, and the values for AFPD, BFPDA, and PFPA, using Adobe Photoshop CS4 (Adobe Systems Inc). In order to evaluate the reliability of radiographing in counting IBs, we used 20 individuals of Sinocyclocheilus grahami based on the method described above, and then, we dissected these same fishes and counted the IBs for comparisons ( Figure 1b).

| Determination of fitting curve
In order to find which function model has the best fit to the IB data collected, a fitting curve estimation was employed in SPSS 19 (SPSS for Windows). We selected all function models on the "Models" option to include as much information as possible and then used the parameter of Rsq, a representative of the statistic value of R 2 , to distinguish the best-fitting curve. The Rsq closest to 1 was identified as the optimal regression model from these various fitting curve categories (such as linear, quadratic, cubic, or inverse). Using the optimal regression model (fitting curve) that was identified from SPSS as the equation category, we drew scatter plots in SigmaPlot 12.5 (SPSS Inc.) and then obtained the specific equation by choosing the optimal curve with R 2 closest to 1.

| Reconstruction of phylogenetic relationship
In consideration of the many polyploid species in Cyprinidae (Yang, Sado, et al., 2015), we reconstructed the phylogenetic relationship using only mitochondrial genes. For each species that we used for IB counting, sequences of three mitochondrial genes (CYT b, COI, and ND4) belonging to the species were downloaded from NCBI (Appendix S1). We aligned and trimmed sequences in MEGA6 (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013) and then combined the three gene fragments into a concatenated dataset. The phylogenetic tree was inferred using the Bayesian inference (BI) method with partitioned models in MrBayes 3.1 (Huelsenbeck & Ronquist, 2001). The best model of each gene fragment was calculated using jModelTest version 2.1.6 (Darriba, Taboada, Doallo, & Posada, 2012), under the Akaike information criteria (AIC) as recommended by Posada and Buckley (2004). Posterior probabilities (PP) were based on two independent MCMC runs for 4 × 10 6 generations, with sampling at every 100 generations yielding 40,000 trees.
The first 25% of trees were discarded as burn-in.

| Ancestral reconstruction of the number of IBs
Each species in the phylogenetic tree of the cyprinids was assigned with the total number of IBs that we obtained from radiography. The maximum parsimony method was selected to reconstruct the ancestral state of the total number of IBs using Mesquite 3.10 (Maddison & Maddison, 2016). Then, five separately parallel analyses were performed similar to the ancestral state reconstruction of the total number of IBs, by using the number of two categories of IBs (epineurals and epipleurals) and three regions of IBs (AFPD, BFPDA, and PFPA) as the inputting states.

| Correlations of IBs and ecological factors
In order to find whether any correlations between IBs and ecological factors, we further mapped diets, depth of occurrence, and water velocity of these fishes onto the phylogenetic tree. We collected the ecological information from some unpublished sampling records in KIZ, as well as from the literature (Chen, 1998;Chu & Chen, 1989;Cui, Li, & Chu, 2011;Jia, Chen, Tao, & He, 2013;Yan et al., 2012;Yue, Shan, & Lin, 2000;Zeng & Liu, 2009;Zhang, 1998;Zhang & Dai, 2011;Zheng, Chen, & Yang, 2009). According to the food composition of these fishes as adults, we identified three categories of diet in our samples: carnivorous, omnivorous, and herbivorous. The carnivorous fishes predominantly fed on invertebrates and vertebrates, while the omnivores fed on plants, invertebrates, and vertebrates, and the herbivorous fishes mostly took aquatic plants and algae. We also recognized three levels of water depth among these fishes. They are either epipelagic, mesopelagic, or benthic fishes. In some papers, species called pelagic were classified as epipelagic fishes, and middle-lower ones were reclassified as benthic fishes.
Considering that the velocity of water is also an important factor of habitats, these cyprinids were further divided into three classes: lentic (predominantly found in lentic habitats such as lakes and pools), lotic (mostly occurring in streams or rivers), and lentic/lotic fishes (roughly with similar frequency of occurrence in both lentic and lotic habitats).

| Comparative numbers of IBs from radiograph and anatomy
The total number of IBs counted from X-ray photographs was largely consistent with those obtained using anatomy (Figure 2a). The differences between two counting methods were generally <4, and the variations of 80% of the samples were from 0 to 2 ( Figure 2b). The median values of IBs for all samples was 104.6 using X-ray photograph and 104.3 using anatomy, and the maximum number of IBs was the same (113) in both methods (Figure 2c). In addition, the number of IBs showed no relationship to standard length in either method in the fishes sampled ( Figure 2d).

| Summary of the number of IBs in Cyprinidae
In the X-ray photographs of 592 individuals, only two categories of IBs, the epineural and epipleural, were observed in our specimens.

| Correlation between IBs and vertebrae in Cyprinidae
The

| Phylogenetic relationships of the Cyprinidae in this study
A phylogenetic tree of the Cyprinidae based on three mitochondrial genes (CYT b, COI, and ND4) revealed three major subclades (Clades I, II, and III) in our analysis ( Figure 5).
According to Chen (1998), each of the 12 morphologically putative subfamilies was hypothesized to be a monophyletic group based on the skeletal characters elaborated by Chen, Yue,

| Ancestral reconstruction of IBs
The ancestral state of the total number of IBs, the two categories of

| Correlations of IBs and ecological factors
The states of reconstructed IBs for two categories (epineurals and epipleurals; panel B-C) and three regions (AFPD, BFPDA, and PFPA; panel D-F), as well as three ecological factors (diet, water depth, and water velocity; panel G-I), were mapped onto the character state tree of IBs based on the phylogenetic relationship ( Figure 6).

| D ISCUSS I ON
The presence of IBs in fish muscles is something of a mystery to ichthyologists, as well as being a troublesome issue in aquaculture-especially in countries like China where cultivation of cyprinids plays an important role in people's nutrition (Tang, Han, Mao, Zhang, & Shan, 2016 (Li et al., 2013;Perazza et al., 2017;Vallod & Arthaud, 2009;Xu et al., 2015). In this study, we compared IB counts between anatomy and X-ray photographs based on 20 medium-sized individuals of S. grahami. We confirmed that the results of X-ray photographs were not significantly different from those of anatomy ( Figure 2).
In consideration of the many advantages of X-rays, such as ease of operation, time-saving, no specimen damaging, and especially its applicability to living fishes, we propose that X-ray photographs are superior to the other two methods. The only drawback is the availability of an X-ray machine.
In our work, based on the large number of specimens and species that we examined, we confirm that epicentrals are not present in cyprinid fishes.
The IB counts in this study reveal that the total number of IBs in cyprinids ranged from 73 to 169 (Figure 3). This range is greater than that reported in a previous study by Lv et al. (2007), which we attribute to the more extensive material examined in our study. The IB counts of most species were similar to those reported by previous researches (Dong et al., 2006;Jiang, Yang, & Bao, 2008;Ke et al., 2009;Li et al., 2013;Lv et al., 2007), which suggests that the IBs are relatively stable at the species level. For instance, in this study, the sample size of each species was generally more than 3 and the proportion of intraspecific deviations of IBs were usually <8% (Appendix S1). However, the IBs of artificially bred species seem to be more complicated. The only species sampled from fish farm in this study, S. grahami, showed a wide range of intraspecific deviation (21%). This intraspecific instability was also found in the different strains of artificially bred common carp: There was a significant difference (p < 0.001) among three examined strains, and these differences were found mainly in the posterior regions of the fish (Cao et al., 2014). In dividing the fish body into three regions (Figure 1), we also found that the number of IBs in PFPA was generally greater than those in AFPD and BFPDA, and in some species, these differences could be as great as five times the normal number. Therefore, we suggest that the PFPA region might be the optimal target for selective breeding for fewer IBs in cyprinids.
The number of sarcomeres reflects, to some extent, the number of vertebrae in fishes (Chen et al., 2017). To the authors' knowledge, this paper is the first attempt to investigate the possible evolutionary patterns of IBs from a phylogenetic perspective. To achieve this goal, a strongly corroborated phylogenetic tree was indispensable. In this study, phylogenetic relationship reconstruction of the Cyprinidae was obtained based on three mitochondrial DNA fragments (CYT b, COI, and ND4). The tree topology ( Figure 5) was largely consistent with most of the previous genetic studies derived from different DNA fragments (Wang et al., 2012;Yang, Sado, et al., 2015). The three major subclades (Clades I, II, and III) were also identical with those found by Jiang et al. (2018), who provided a well-resolved inter-subfamily relationship of cyprinids at genome level (although the monophyly of each subfamily was beyond their aims for that paper). In this study, discussion of the deep phylogenetic relationship of the Cyprinidae is not our main focus; however, the tree we present here indicates that most of the morphological-based subfamilies proposed by Chen (1998) are not monophyletic groups ( Figure 5), thus implying this popular 12-subfamily classification (Chen, 1998) needs to be revised. Based on the phylogenetic tree in this study (Figure 5), the number of IBs seemed to have phylogenetic relevance. For instance, the species in Clade III-2, including the nominal subfamilies Hypophthalmichthyinae, Xenocyprinae, and Cultrinae and some species of Leuciscinae and Gobiobotinae, turned out to have similar numbers of total IBs, epineurals, and epipleurals ( Figure 6).
Tinca tinca, a species previously placed in the Leuciscinae, has long been treated as incertae sedis in both morphological and molecular studies (see more discussion in Jiang et al., 2018). Here, this species had a similar IB composition to Acheilognathus, but distinct from most other Leuciscinae species (Figure 6). The close relationship of Tinca tinca and the Acheilognathinae revealed in this study ( Figure 5) was also demonstrated by Yang, Sado, et al. (2015) based on both mitochondrial and nuclear DNA fragments. Percocypris, previously placed in the Barbinae (Chen, 1998), was revealed in our tree as the sister group of Schizothorax, which is consistent with a previous molecular study that suggested Percocypris originated from a common ancestor with Schizothorax (Wang, Yang, & Chen, 2013 (Danos & Staab, 2010;Dong et al., 2006;Lv et al., 2007). However, when we looked at omnivorous fishes, the situation was more complicated. For example, the closely phylogenetically related species, P. parva, S. dabryi, and A.
rivularis, displayed significantly divergent IB numbers. We did not observe any correlation between IBs and water depth or velocity.
We infer that the evolution of IBs in the Cyprinidae would probably not be affected by single ecological factor, but rather be a product of multiple factors. For instance, we compared the number and shape of IBs in Sinocyclocheilus, a group that is characterized by containing many cave species. We found the typical cave-dwelling species (known as troglobites) had fewer IBs than those atypical cave species (known as troglophiles, see also in Yang, Pan, et al., 2015). Caves are an extreme environment without light, which embody a series of different ecological factors relative to open water (e.g., scarce of food supply, low water fluctuation, and low predation pressure).
In summary, the variations of IBs in the modern Cyprinidae were full of diversity, which may reflect their phylogenetic history, but has then been shaped by multiple environment factors. Botany, Chinese Academy of Sciences) for taking X-ray photograph of some of the samples.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflicts of interest to this work.

AUTH O R CO NTR I B UTI O N S
KFY, WSJ, and JXY designed the work and scientific objectives. KFY, XAW, YWZ, and XFP collected and prepared the fish samples. KFY analyzed the data. KFY and WSJ prepared the manuscript. JXY revised the manuscript. All authors have read and approved the final manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
DNA sequences: All sequences used in this study were downloaded from NCBI (see accession numbers in Appendix S1). The following datasets have been deposited in Dryad http://doi.org/10.5061/ dryad.3jt83k4. Dataset S1. DNA sequence alignment for the concatenated dataset with three genes in NEXUS format. Dataset S2.
Input file for Bayesian inference (BI) analyses of the concatenated dataset in NEXUS format. Dataset S3. Input file for Mesquite analysis in NEXUS format.