Specimen sampling

For behavioral analysis, adult black chafer H. kiotonensis (Brenske, 1894; Fig. 1A), H. convexopyga (Moser, 1912; Fig. 1B), and H. picea (Waterhouse, 1875; Fig. 1C) were collected from the Toyonaka campus of Osaka University in Osaka, Japan (34°48′N, 135°27′E) from April to September 2016. Additional samples of H. convexopyga and H. picea were collected in May 2020. Adults of H. loochooana loochooana (Sawada, 1950) were collected from Nisi-Henna-Zaki (24°54′N, 125°15′E) and Hirara-Higashi-Nakasone (herein called Manabi-No-Mori, 24°48′N, 125°17′E) (Fig. 1E) Miyako Island, Okinawa, Japan, from June to August 2016 and 2018.

For phylogenetic analysis, in addition to the species mentioned above, H. loochooana okinawana (Nomura, 1964; Fig. 1D) were collected in Itoman City (26°7′N, 127°40′E), Okinawa, Japan, in June 2016; H. parallera (Motschulsky, 1854; Fig. 1F) on the bank of the Yodo River (34°43′N, 135°31′E), Osaka, Japan, in July 2016; and Miridiba castanea (Waterhouse, 1875) on the Toyonaka campus of Osaka University in May 2016. The collection sites, dates, and sex of all samples are shown in Table 1. Miridiba castanea was used as an outgroup. Holotrichia loochooana loochooana, H. loochooana okinawana, and their food, the beach cabbage, Scaevola taccada, were collected on Miyako Island, and transported by air. All collected beetles were kept in an air-conditioned room at 23–27°C under natural daylight before use. Holotrichia loochooana loochooana and H. loochooana okinawana were fed leaves of S. taccada and other species were fed the Japanese cherry Prunus yedoensis ‘Somei-yoshino’.

Table 1.

Collection and information of specimens used for phylogenetic analysis.

Activity recording

The activity recording method was adopted from Kawasaki et al. (2017) and Watanabe and Shiga (2020). The activity of the beetles on the ground was individually recorded under periods of 12-h light and 12-h darkness (LD) for 10 days and subsequent constant darkness (DD) for another 14 days at 25 ± 1°C. A fluorescent lamp with an intensity of 1.35 W/m² (Panasonic FL15W; Panasonic, Osaka) was used as the light source. During the recording period, a leaf of P. yedoensis was provided to H. kiotonensis, H. convexopyga, and H. picea, and a leaf of S. taccada to H. loochooana loochooana, every 5 days. Approximately 20–25 mL of water was sprayed on the soil surface after it became dry.

The top of the container was covered with a transparent glass plate. Color or monochromatic images of the soil surface were taken every 1 min with a web camera (DC-NCR13U; Hanwha Japan, Tokyo). Beetle movements were detected automatically through differences in the total pixel values between two serial photos and eye inspections plotted on a double-plotted actogram (Kawasaki et al., 2017; Watanabe and Shiga, 2020). The presence of rhythmicity in beetle activity was determined separately in LD and DD using a chi-square periodogram analysis (Enright, 1965; Sokolove and Bushell, 1978).

DNA sequencing

Collected beetles were sacrificed, stored in 99.5% ethanol, and placed in a freezer until subsequent experiments. Genomic DNA was extracted from the hind leg, and purification of total DNA from insects was performed using a DNeasy^® Blood & Tissue Kit (QIAGEN, Hilden, Germany) following the manufacturer's instructions. Polymerase chain reaction (PCR) was conducted to amplify the cytochrome c oxidase subunit 1 (CO1) (Former et al., 1994) and 16S ribosomal RNA (16S rRNA) (Hosoya et al., 2001) genes in the mitochondrial DNA and the histone 3 gene (H3) (Ogden and Whiting, 2003) in the nuclear DNA using a GoTaq^® Green Master Mix (Promega, Madison, Wisconsin, United States). The primer sequences are listed in Table 2. The PCR amplification temperatures were as follows: for CO1, initial denaturation at 94°C for 2 min, followed by 30 cycles at 94°C for 15 s, 45°C for 30 s, and 72°C for 1 min; for 16S rRNA and H3, initial denaturation at 94°C for 2 min, and then 30 cycles at 94°C for 15 s, 50°C for 30 s, and 72°C for 1 min. The product yield was monitored via electrophoresis using a 1.5% agarose gel (code No. 312-01193, Nippon Gene, Toyama), and the amplicons were purified using a Monarch PCR & DNA Cleanup Kit (New England Biolabs Inc., Ipswich, MA).

Table 2.

Primer sequences used in the study.

Purified samples were prepared as 10–40 ng of CO1 and 16S rRNA PCR products and 6–20 ng of H3 PCR products using a NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA). Sequencing was outsourced (Fasmac Co., Ltd., Kanagawa).

Phylogenetic analysis

We used three datasets for phylogenetic analysis: H3, CO1, and 16S rRNA. After removing data of beetles that lacked any of the three sequences, the H3, CO1, and 16S rRNA sequences were concatenated and aligned using T-Coffee ( http://tcoffee.crg.cat/apps/tcoffee/do:mcoffee; Notredame et al., 2000) and trimmed using TrimAl (Capella-Gutiérrez et al., 2009) to obtain 29 1527-bp long sequences. Aligned sequences were checked manually to find the 1^st to 3^rd codon positions of H3 and CO1. The optimal partitioning strategy was assessed using PartitionFinder (Lanfear et al., 2017) to obtain seven partitioned datasets for maximum likelihood estimation and the Bayesian inference (see   Supplementary Table S1 (zs210091_TableS1.pdf)).

Using partitioned datasets, a maximum likelihood estimation was conducted using the IQ-Tree online version ( http://iqtree.cibiv.univie.ac.at/; Trifinopoulos et al., 2016) with 1000 bootstrap replications (Efron, 1982; Felsenstein, 1985). Models were selected using the Bayesian Information Criterion with the ModelFinder subprogram in IQ-Tree for partitioned sequences (see   Supplementary Table S1 (zs210091_TableS1.pdf)). Phylogenetic trees and bootstrap values were then obtained. The Bayesian inference method was performed using Bayesian MrBayes 3.2.7 (Ronquist et al., 2012). The selection of rate-setting models was made using PartitionFinder (Lanfear et al., 2017) (see   Supplementary Table S1 (zs210091_TableS1.pdf)). Markov chain Monte Carlo analyses were run for 1,000,000 generations with tree sampling every 5000 generations, and convergence was obtained. The first 25% of the sampled trees were discarded as burn-in trees.

Statistical analyses

The Shapiro-Wilk test was used to examine if the data set was drawn from a normal distribution. Bartlett's test was used to test whether the samples had equal variances. For nonparametric tests, the Exact Wilcoxon rank-sum test (for data sets with homoscedasticity), Fligner-Policello test (for datasets with heteroscedasticity), or the Steel-Dwass test was used. All tests were conducted using R software ver. 4.0.3 ( http://cran.r-project.org, Retrieved 10 October 2020) with an additional package (“ http://aoki2.si.gunma-u.ac.jp/R/src/Steel-Dwass.R”, encoding = “euc-jp”, Retrieved 26 November 2021).

Activity rhythm

We recorded the activity rhythms of 13 H. kiotonensis, 14 H. convexopyga, and 25 H. loochooana loochooana collected at Manabi-No-Mori, 17 H. loochooana loochooana collected at Nishi-Henna-Zaki, and 14 H. picea. The collection sites, sex, and behavior types of the beetles used in activity recording are listed in   Supplementary Table S2 (zs210091_TableS2.pdf). In H. loochooana loochooana, there was no significant difference in the free-running period between males and females (H. loochooana loochooana Manabi-No-Mori: Fligner–Policello test, P = 0.369; H. loochooana loochooana Nishi-Henna-Zaki: Wilcoxon–Mann–Whitney test, P = 0.832). In other species, no noticeable difference of the activity rhythm between sexes was observed, although statistical tests were not applicable because of the small sample size in this study. Therefore, we showed the results for males and females together in Figs. 2 and 3.

Holotrichia kiotonensis were mainly active during the night period, but their rhythmicity was unclear compared with other species. Sometimes their activity occurred during the light period of the LD cycles (Fig. 2A). During DD only 38% of beetles showed circadian rhythms, with a free-running period of 22.6 ± 0.8 h (mean ± S.D.) (Fig. 3). The other beetles exhibited arrhythmicity in DD.

Holotrichia convexopyga showed nocturnal activity in LD cycles. They walked around and ate leaves on the ground during the dark period but stayed mostly underground in the light period. All beetles exhibited circadian rhythm in DD with a period of 23.2 ± 0.6 h (Figs. 2B, 3).

Holotrichia loochoana loochoana collected from Manabi-No-Mori and Nishi-Henna-Zaki both exhibited circadian rhythms in DD (Fig. 3A). However, differences between the two populations were found in the light period of the LD cycles. Holotrichia loochoana loochoana from Manabi-No-Mori started major activities at the onset of darkness (Fig. 2C₁, and see   Supplementary Figure S1A (zs210091_FigS1.pdf)), whereas those from Nishi-Henna-Zaki started to increase activities approximately 4–5 hours before light-off (Fig. 2C₂, and see  Supplementary Figure S1B (zs210091_FigS1.pdf)). The free-running period under DD of the Manabi-No-Mori population (23.4 ± 0.4 h) was slightly but significantly longer than that of the Nishi-Henna-Zaki population (22.8 ± 0.8 h) (Steel-Dwass test, P = 0.038; Fig. 3B).

In H. picea, most beetles exhibited nocturnal activity every 2 days under LD, with 92.9% of them (n = 14) showing circabidian rhythm in DD (Figs. 2D, 3A). In beetles with circabidian rhythms in DD, some switched appearance from the odd days (day 1, 3, 5, and so on) to the even days (day 2, 4, 6, and so on), or vice versa, as counted from the start of the recording (day-switching). One beetle switched appearance from odd to even days by emerging in two consecutive nights (day 15 and day 16) in DD (Fig. 2D₂). Another beetle switched twice in DD (Fig. 2D₃). During DD, 50.0% of the beetles showed circabidian rhythms with day-switching (Fig. 3A). Day-switching was also observed during the LD cycles in four of the 14 beetles (see   Supplementary Table S2 (zs210091_TableS2.pdf)). Interestingly, among the beetles showing circabidian rhythms, two showed a circadian-like activity component. In one beetle (Fig. 2D₄), major activity was observed in the even days of the LD cycles and in DD until day 16, while a minor activity component appeared on the odd days in the same period. In the minor activity, the beetle emerged on the ground briefly at the beginning of the subjective nights, followed by burrowing back into the soil. From day 17, the occurrence of major and minor activities switched (Fig. 2D₄). This beetle was classified as circabidian with day-switching (Fig. 3A). In the other beetle, a minor component appeared only in the last four circabidian cycles during DD and was classified as circabidian without day-switching in Fig. 3A. The free-running period of the circabidian rhythm without day-switching (n = 6) was 47.7 ± 0.3 h, which was significantly longer than that of other species (Steel-Dwass test, P < 0.05, Fig. 3B).

Fig. 2.

Double-plotted actogram (left) and the associated chi-square periodogram analysis (right) showing locomotor activity rhythm on the ground in adult Holotrichia species. (A) Holotrichia kiotonensis (female) showing circadian rhythm. (B) Holotrichia convexopyga (male) showing circadian rhythm. (C₁) Holotrichia loochooana loochooana (male) from Manabi-No-Mori and (C₂) Holotrichia loochooana loochooana (male) from Nishi-Henna-Zaki showing circadian rhythm. (D₁) Holotrichia picea (female) showing regular circabidian rhythm, (D₂) Holotrichia picea (female) showing circabidian rhythm with day switching (once), (D₃) Holotrichia picea (sex not determined) showing circabidian rhythm with day switching (twice), (D₄) Holotrichia picea (female) showing circabidian rhythm with a circadian-like pattern. Periodogram analysis was conducted for the activity during the DD period indicated with a gray vertical bar on the actogram. The peak values are indicated.

Fig. 3.

Summary of activity rhythms in Holotrichia kiotonensis, Holotrichia convexopyga, Holotrichia loochooana loochooana Manabi-No-Mori (Manabi), Holotrichia loochooana loochooana Nishi-Henna-Zaki (Nishi), and Holotrichia picea. (A) Proportion of beetles exhibiting circadian rhythm, circabidian rhythm without day-switching, circabidian rhythm with day-switching, and arrhythmic patterns in their locomotor activity during constant darkness. (B) Average and standard deviation of the free-running period during constant darkness. Columns with different letters show significant differences in the free-running period (Steel-Dwass test, P < 0.05).

Phylogenetic analysis

We obtained a 682.83 ± 0.35 (average ± S.E.)-bp sequence for CO1, 521.24 ± 0.26-bp sequence for 16S rRNA, and 344.14 ± 0.35-bp sequence for H3 from 29 individuals (Table 1, Accession nos. LC677229–LC677315). Lee et al. (2015) showed that the nodal supporting values on the phylogenetic tree, bootstrap values, and Bayesian posterior probabilities were improved using the concatenated sequences of CO1 and 16S rRNA compared with those based on their respective single sequences. Therefore, we employed concatenated sequences of the H3, CO1, and 16S rRNA genes for phylogenetic analysis. There was distinct separation into two major clades, clades A and B, in both the maximum likelihood estimation analysis and Bayesian inference analysis with strong supporting values in the Holotrichia group (Fig. 4). Six ingroup taxa representatives were clustered into six monophyletic groups corresponding to nominal species, except for H. picea, according to the maximum likelihood estimation (Fig. 4A). In the Bayesian inference method, H. picea was clustered in one group (96%, Fig. 4B). Clade A included H. loochooana loochooana, H. loochooana okinawana, H. kiotonensis, and H. convexopyga. Two subspecies of H. loochooana were separated into distinct groups. They formed one cluster with moderate supporting values, a bootstrap value of 89, and a Bayesian inference probability of 66. The Manabi-No-Mori and Nishi-Henna-Zaki groups were not separated in the H. loochooana loochooana cluster (Fig. 4). Holotrichia kiotonensis clustered with the H. loochooana clade (bootstrap value of 91 and Bayesian inference probability of 66). Clade B was composed of H. picea and H. paralella. High nodal supporting values in clade B were obtained using the Bayesian inference method.

Fig. 4.

Phylogenetic relationships among Holotrichia species reconstructed using the maximum likelihood estimation (A) and Bayesian inference method (B) from concatenated sequences of the histone H3, cytochrome c oxidase subunit 1, and 16S ribosomal RNA genes. Numbers above the branches show bootstrap values in (A) and Bayesian posterior inferences in (B). Holotrichia loochooana loochooana with and without an asterisk are collected at Manabi-No-Mori and Nishi-Henna-Zaki, respectively. Scale bars indicate 5% nucleotide mutation rate.