Diverse clades of fishes adapted to feeding on the benthos repeatedly converge on steep craniofacial profiles and shorter, wider heads. But in an incipient radiation, to what extent is this morphological evolution measurable and can we distinguish the relative genetic vs. plastic effects? We use the Trinidadian guppy (Poecilia reticulata) to test the repeatability of adaptation and the alignment of genetic and environmental effects shaping poecilid craniofacial morphology. We compare wild-caught and common garden lab-reared fish to quantify the genetic and plastic components of craniofacial morphology across four populations from two river drainage systems (n=56 total). We first use microCT to capture 3D morphology, then place both landmarks and semilandmarks to perform size-corrected 3D morphometrics and quantify shape space. We find a measurable, significant, and repeatable divergence in craniofacial shape between high predation invertivore and low predation detritivore populations. As predicted from previous examples of piscine adaptive trophic divergence, we find increases in head slope and craniofacial compression among the benthic detritivore foragers. Furthermore, the effects of environmental plasticity among benthic detritivores produces exaggerated craniofacial morphological change along a parallel axis to genetic morphological adaptation from invertivore ancestors. Overall, many of the major patterns of benthic-limnetic craniofacial evolution appear convergent among disparate groups of teleost fishes.

Collection and animal husbandry

We collected female and male guppies, Poecilia reticulata, from four distinct populations in the Northern Range mountains on the island of Trinidad in spring 2013 and spring 2014. Guppies were collected from two independent drainages: Aripo river (southern slope) and Yarra river (northern slope) (Figure 1A); these watersheds are completely isolated from one another, with no admixture occurring over the mountain divide (but some can occur within each drainage).

The specific populations collected were from benthic detritivore, or low predation, and invertivore, or high predation, sites within the Aripo and Yarra drainages (four populations total). Twenty females and five males were collected at each of the four sites (n = 100 total); the fish were transported in bottles filled halfway with fresh stream water and treated with 2mg Tetracycline (1mg/L) and StressCoat (according to manufacturer directions) to minimize infection during transportation. The fish were housed in an open-air field laboratory in Verdant Vale, Arima, Trinidad, within 2-liter tanks each equipped with air stone at (~18+3 degrees C) temperatures and exposed to ambient light, which maintained a light:dark schedule at approximately 12:12 hours. Fish were fed twice daily on algae flakes (TetraMin Tropical Flakes) in rations that amounted to two minutes of feeding time per meal (Arendt & Reznick, 2005; D. N. Reznick, 1982). 30% water changes were performed weekly throughout the husbandry period.

We randomly selected 10 female guppies from each population to represent the “wild-caught” phenotype and euthanized them as quickly as possible after capture to minimize plastic response to laboratory settings. We did not record the exact duration that we had the fish in captivity, but it was often less than a week and never more than 30 days.  These fish were euthanized via overdose of buffered 1 g/L tricaine methanesulfonate (Tricaine-S, Western Chemical Inc., Ferndale, WA, USA), fixed for 24 hours in 4% paraformaldehyde and preserved in 70% ethanol at -20 degrees C.

A subset of live female (n = 10) and male (n = 5) guppies from each of the four populations (n=60 total) were exported from Trinidad with a permit issued by the Fisheries Division of the Ministry of Agriculture, Land and Fisheries, Republic of Trinidad and Tobago, and imported into the United States with a valid U.S. Fish and Wildlife Service Declaration of Importation permit (USFWS Form 3-177). These fish were bred to produce the “lab-reared” treatment of the study.

Lab-reared fish were housed in community aquaria, which are 38-liter tanks with open stocks of fish. Population density was regulated to reduce overcrowding. Open stock temperature was maintained at 25 to 26C and fish were fed newly hatched Artemia, Repashy Community Crave (a gel that is prepared in small amounts every few days), and occasional Omega one flakes, Daro worms and Daphnia, which were cultured in the laboratory. This mixed diet was used to minimize strong directional plasticity. It should be noted that we did not conduct a reciprocal feeding experiment; fish from all four represented populations were fed the same food throughout their life in the lab, which provides a common garden environment from which to compare the genetic basis of craniofacial traits. Lab-reared fish were housed for a minimum of two generations to remove any epigenetic effects. We selected 5-10 individual lab-reared female guppies from each of the four populations. These lab-reared fish were euthanized, fixed and preserved following the above procedure.

All research reported here followed strict ethical guidelines and complied with the

US federal government’s regulations. Procedures were approved by the Institutional Animal Care and Use Committee at Brown University (protocol: 1211035 to E. L. Brainerd), Harvard University (protocol: 20–03-2 to G. V. Lauder) and University of California Riverside (#:20200003 to D. N. Reznick).

CT scanning and segmentation

            We scanned all specimens using a micro Computed Tomography scanner (Bruker SkyScan 1173, Kontich, Belgium). To scan a group of fish, we first wrapped 1-7 individuals in Kimwipes and placed them within 5 mL transport tubes with a small amount of ethanol to keep the specimens moist and stabilized within the tube. We µ-CT scanned each tube using 0.5° increments over 180° with the x-ray source set at 38kV and 190µA. This resulted in image stacks with an isotropic resolution between 6.8µm and 9.0µm that were then reconstructed using NRecon (Bruker, Kontich, Belgium). We used Mimics (v. 22.0; Materialise NV, Leuven, Belgium) to segment individuals from each scan and isolate the pre-pectoral skeleton, then exported scans as 3D models for morphometric analysis. We removed any scans where the heads showed obvious distortion such as a lateral skew of the craniofacial morphology. Finally, we used MeshLab (v. 2020.07; Cignoni et al., 2008) to simplify our meshes (Quadratic Edge Collapse Decimation filter) to reduce necessary computing power for visualization. We used Python to batch process the meshes and reduced each scan to 1x10⁶ faces.

3D Geometric Morphometric Analysis

We placed 11 fixed landmarks and 21 sliding semilandmarks across the craniofacial apparatus to capture its three-dimensional shape (Figure 2; supplemental figure 1). We use these landmarks to capture elongation and compression across all three spatial axes, as well as the craniofacial slope, orientation of the mouth (including the articulation of the lower jaw and the quadrate), and shape of the suspensorium. While many past studies have included landmarks on the oral jaws, we exclude these points due to confounding variation induced from inconsistent positioning during fixation. However, each of the traits that we measure has been used to infer functional differences between populations or species of fish that have diverged along a trophic axis (Albertson et al., 2003b, 2005; Cooper et al., 2010, 2011). We restricted landmarks to the left side of the head and the midline to avoid capturing asymmetries within individuals. In cases where the left side of the head was damaged we instead landmarked the right side and reflected the points about the x axis to correct the chirality. We collected sliding semilandmarks by first placing an arbitrary number of points along either the maxilla or preopercle to fully capture the curvature of the structure, then resampling these curves to a fixed number of equally spaced semilandmarks (8 on the maxilla, 13 on the preopercle) bounded by two fixed landmarks on each structure. All landmarking was done in 3D Slicer (v. 4.11; Fedorov et al., 2012).

            We created three data sets, one with all the data points, one with just the lab-reared fish, and one with just the wild-caught fish. Using Geomorph (v. 4.0.0; D. Adams et al., 2021) we conducted a geometric morphometric analysis on each dataset. First, we subjected each coordinate data set to a generalized Procrustes analysis using bending energy. This method compares the shapes of each individual after accounting for size, rotation, and position. Next, we regressed the shape variable on log-transformed geometric centroid size for each of our data sets. In all cases there was a significant allometric effect, so we took the residuals of these regressions to attain an allometry-free measure of shape (Outomuro & Johansson, 2017). For each data set we then used the function procD.lm to run multivariate regressions of allometry-corrected shape on feeding niche, drainage, and the interaction effect between the two. For the data set containing both lab-reared and wild-caught fish we also included the rearing condition as an explanatory variable. To visualize the distribution of each data set in shape space the shape data were subjected to a principal component analysis (PCA) using the gm.prcomp function and plotted in the morphospace defined by the first two principal component axes. Finally, we subjected the lab-reared and wild-caught data sets to a trajectory analysis (RRPP, V1.3.1)(Collyer & Adams, 2022) to determine whether the evolved path through morphospace from high predation to low predation populations differed between river drainage systems. We plotted population mean shapes and the linear trajectories between them to visualize their evolutionary paths.

The regression and PCA analyses on the full dataset revealed fish differed based on feeding niche and drainage, but that the rearing condition primarily affected the magnitude of the divergence and not the direction of the divergence in shape space. Therefore, we only visualized the shape differences between benthic detritivore and invertivore individuals in the wild-caught individuals since these were further diverged than were the lab-reared fish. We did this by plotting the size-corrected and allometry-corrected mean position of each landmark for benthic detritivore and invertivore fish separately for each drainage. We placed lines between landmarks on the same bone and filled in the area between homologous lines to aid in visualization. All geometric morphometric analyses were run using Geomorph (v. 4.0.0)(D. Adams et al., 2021) in RStudio (v.1.3.1056).

Although these landmarks inform craniofacial elongation among these populations, they do not allow us to measure head length since this is typically measured from the oral jaws to the caudal end of the neurocranium. Instead, we infer head length based on the distance along the anterior-posterior axis from the dorsal most point on the posterior end of the preopercle (Supp. Figure 2) to the dorsal tip of the maxilla. To compare our results using this nontraditional metric to past studies on head length we additionally measured the distance from the dorsal tip of the maxilla to the posteroventral margin of the basioccipital. We ran a linear regression between these two measurements to check if our metric of craniofacial elongation was positively correlated with head length. Similarly, we compared our measurement of head width, the distance from the widest point of the preopercle to the midline plane, to a more common measurement of width, the distance from the lateral most point on the operculum (Supp. Figure 2) to the midline plane of the fish.