How gene expression diverges during adaptation might be strongly influenced by the geographic setting and timeframe over which species evolve. To contrast transcriptomic patterns of replicate trophic adaptations that evolved convergently during both allopatric and sympatric contexts, we conducted RNA-seq on the trophically important lower pharyngeal jaws of two sympatrically and four allopatrically diverged species pairs of cichlid fishes. We first show that all of these species pairs have convergently diverged along a crushing trophic axis and that the sympatric pairs are as phenotypically divergent as the allopatric pairs. Then, we found that distinct sets of genes were differentially expressed in the jaws of sympatrically diverging pairs as compared to jaws in older allopatric species pairs. The genes that were differentially expressed in the jaws of allopatric pairs also were more highly expressed on average than in the sympatric pairs. Finally, for genes that were differentially expressed, the magnitude of differences in expression between the jaws were greater for sympatrically diverging species pairs. The particular genes, their expression levels, and the magnitude of expression differences between sympatrically originating adaptations might all play an important role in generating and maintaining boundaries to gene flow during the rapid ecological divergence that often characterizes sympatric speciation.

To adjust for body size, we digitally photographed the lateral view of each fish with a ruler in the image prior to dissection of the LPJ and its storage at -20°C in RNAlater^TM (Sigma-Aldrich, St. Louis, USA). All fish were adult size and ranged from ~50 mm to ~80 mm based on subsequent digital measurements of standard length using ImageJ v2.0.0-rc-69/1.52i (Schindelin et al., 2012). Following photography of the whole fish, we euthanized them with an overdose of MS222 and removed their LPJs. After dissecting the LPJ from the branchial cavity, we immediately removed surrounding connective and muscle tissue from each jaw while they were submerged in RNAlater^TM. Then, we photographed the dorsal and ventral surfaces of the LPJs using a Leica MZ10F stereomicroscope with a size standard.

Using the LPJ digital images, multiple measurements were obtained with ImageJ.  Using the image of the toothed dorsal aspect of the LPJ, we measured the area of the six most posterior teeth along the midline which tend to be the six largest teeth on the LPJ. The tooth number of each LPJ was also determined. To linearize the average tooth areas for each individual, this measurement was square root transformed and adjusted for size by taking residuals from the regression of values on individual standard length for each species pair. Then, with the ventral image of the LPJ, we measured the total length of the LJP along its midline, the linear amount of interdigitation between the two halves of the LPJ, and the total length of the interdigitation along the suture. The total amount of suturing of each LPJ, P_L, was then determined according to. Effectively, the total LPJ length, L_j, the suture length L_s along the jaw midline, as well as the length of the interdigitation, L_i, were combined using the equation:

P_L = (Ls/Lj)/(L_i/L_s)

RNA was extracted from the jaws for high throughput sequencing of their trasncriptomes. Importantly, pharyngeal teeth are continuously replaced and the lower pharyngeal jaw bones capable of actively remodeling, so the toothed jaws of the size of fishes used should show substantial transcriptional activity. Furthermore, all fish were raised in aquaria on similar flake food diets. Therefore, the transcriptional differences observed should generally be constituitive differences and not phenotypic plastic responses to diet, For all LPJs used in this study, we extracted RNA using the ReliaPrep miRNA Cell and Tissue Miniprep System (Promega, Madison, USA) after flash freezing LPJs in liquid nitrogen and grinding them to a fine powder for further processing accordint to the manufacturers protocol. Concentration and integrity of RNA was measured using a Qubit^TM 2.0 Fluorometer (Invitrogen by Thermo Fischer Scientific, Oregon, USA) and a Bioanalyzer 2100 system (Agilent Technologies, Waldbronn, Germany), respectively.

We used 200 ng of total RNA per sample for the library preparations that employed the SENSE mRNA-Seq Library Prep Kit V2 (LEXOGEN, Vienna, Austria). To produce the libraries, we followed the manufacturers protocol for a sequencing length of 150bp and using 17 PCR cycles. For quality control, DNA concentration and quality were assessed using a Qubit^TM 2.0 and an Agilent Tape Station 4150, respectively. After equimolar pooling of RNA libraries, pools were size-selected in a range of 300-600bp and sequenced on an Illumina HiSeq X-ten platform (paired end 2x 150bp) at the Beijing Genomics Institute (BGI, Hong Kong). Raw reads were processed with Trimmomatic v0.36 (default parameters) for quality control and to remove remaining adapters. Using STAR v2.7, trimmed reads were first aligned to the reference genome of Amphilophus citrinellus and then quantified at the gene level for each of the 22495 annotated protein-coding genes.