The evolution of Batesian mimicry – whereby harmless species avoid predation through their resemblance to harmful species – has long intrigued biologists. In rare cases, Batesian mimicry is linked to intraspecific colour variation, in which only some individuals within a population resemble a noxious ‘model’. Here, we assess intraspecific colour variation within a widespread New Zealand stonefly, wherein highly melanised individuals of Zelandoperla closely resemble a chemically defended aposematic stonefly, Austroperla cyrene. We assess convergence in the colour pattern of these two species, compare their relative palatability to predators, and use genome-wide association mapping to assess the genetic basis of this resemblance. Our analysis reveals that melanised Zelandoperla overlap significantly with Austroperla in colour space, but are significantly more palatable to predators, implying that they are indeed Batesian mimics. Analysis of 194,773 genome-wide SNPs reveals an outlier locus (ebony) strongly differentiating melanic versus non-melanic Zelandoperla. Genotyping of 338 specimens from a single Zelandoperla population indicates that ebony explains nearly 70% of the observed variance in melanism. As ebony has a well-documented role in insect melanin biosynthesis, our findings indicate this locus has a conserved function across deeply divergent hexapod lineages. Distributional records suggest a link between the occurrence of melanic Zelandoperla and the forested ecosystems where the model Austroperla is abundant, suggesting the potential for adaptive shifts in this system underpinned by environmental change.

Photographs were taken in RAW format using a Canon 40D camera with an EF-S 60 mm f/2.8 lens under fixed magnification, exposure time, aperture, and ISO sensitivity settings. The average red, green, and blue (RGB) colour channel values of each image were then linearised and equalised in reference to the nominal RGB values of the grey patches of the ColorGauge Pico target in Adobe Photoshop CC 2021. Colour measurements were taken from calibrated images in five regions of interest (ROIs) using the RGB Measure function in ImageJ. The first three chosen ROIs correspond to areas of the head (1), pronotum (2), and femorae (3), to provide a broad representation of overall body melanisation. The additional ROIs correspond to areas on the foreleg tibiae (4) and forewing bases (5) associated with warning colouration in the aposematic model Austroperla. Mean RGB values from each ROI were converted to HSB (hue, saturation, brightness) values using the colorsys module in Python, and only the achromatic brightness values of the five ROIs, within a range of 0% (black) to 100% (white), were used in subsequent analyses. A principal component analysis (PCA) was performed on brightness data from the three ROIs associated with overall body melanisation in a combined dataset (n = 360) containing Zelandoperla and the aposematic model Austroperla.