Ultraconserved elements resolve the phylogeny and corroborate patterns of molecular rate variation in herons (Aves: Ardeidae)

Abstract

Thoroughly sampled and well-supported phylogenetic trees are essential to taxonomy and to guide studies of evolution and ecology. Despite extensive prior inquiry, a comprehensive tree of heron relationships (Aves: Ardeidae) has not yet been published. As a result, the classification of this family remains unstable, and their evolutionary history remains poorly studied. Here, we sample genome-wide ultraconserved elements (UCEs) and mitochondrial DNA sequences (mtDNA) of >90% of extant species to estimate heron phylogeny using a combination of maximum likelihood, coalescent, and Bayesian inference methods. The UCE and mtDNA trees are mostly concordant with one another, providing a topology that resolves relationships among the 5 heron subfamilies and indicates that the genera Gorsachius, Botaurus, Ardea, and Ixobrychus are not monophyletic. We also present the first genetic data from the Forest Bittern Zonerodius heliosylus, an enigmatic species of New Guinea; our results suggest that it is a member of the genus Ardeola and not the Tigrisomatinae (tiger herons), as previously thought. Finally, we compare molecular rates between heron clades in the UCE tree with those in previously constructed mtDNA and DNA–DNA hybridization trees. We show that rate variation in the UCE tree corroborates rate patterns in the previously constructed trees—that bitterns (Ixobrychus and Botaurus) evolved comparatively faster, and some tiger herons (Tigrisoma) and the Boat-billed Heron (Cochlearius) more slowly, than other heron taxa.

Resumen

Los árboles filogenéticos cuidadosamente muestreados y bien respaldados son esenciales para la taxonomía y para guiar los estudios de evolución y ecología. A pesar de una extensa investigación previa, aún no se ha publicado un árbol completo de las relaciones de las garzas (Aves: Ardeidae). Como resultado, la clasificación de esta familia sigue siendo inestable y su historia evolutiva sigue siendo poco estudiada. Aquí, tomamos muestras de elementos ultraconservados (EUCs) de todo el genoma y secuencias de ADN mitocondrial (ADNmt) de >90% de las especies existentes para estimar la filogenia de las garzas usando una combinación de métodos de máxima verosimilitud, coalescencia e inferencia bayesiana. Los árboles de EUC y ADNmt son en su mayoría concordantes entre sí, lo que proporciona una topología que resuelve las relaciones entre las cinco subfamilias de garzas e indica que los géneros Gorsachius, Botaurus, Ardea e Ixobrychus no son monofiléticos. También presentamos los primeros datos genéticos de Zonerodius heliosylus, una enigmática especie de Nueva Guinea; nuestros resultados sugieren que es un miembro del género Ardeola y no de Tigrisomatinae (garzas tigre), como se pensaba anteriormente. Por último, comparamos las tasas moleculares entre los clados de garzas en el árbol de EUC con aquellas de los árboles de ADNmt y de hibridación ADN–ADN construidos previamente. Mostramos que la variación de las tasas en el árbol de EUC corrobora los patrones de las tasas en los árboles construidos previamente—que Ixobrychus y Botaurus evolucionaron comparativamente más rápido y algunas garzas tigre (Tigrisoma) y Cochlearius más lento que otros taxones de garzas.

Graphical Abstract

Lay Summary

• We use genetic data from across the genome and produce a robust family tree for herons, which clarifies the relationships among subfamilies and genera.

• A comprehensive phylogeny of herons is lacking. As a result, their taxonomy is unstable and their evolutionary history is poorly known.

• Several species were found to be incorrectly classified, and we recommend appropriate taxonomic revisions.

• Comparisons of molecular evolution support previous studies. Bitterns have evolved comparatively faster, with some tiger herons and the Boat-billed Heron having evolved comparatively slower.

INTRODUCTION

Herons (Ardeidae) are wading birds in the order Pelecaniformes (Hackett et al. 2008, Jarvis et al. 2014, Prum et al. 2015) and consist of 61–66 species and 18–19 genera depending on classification (Kushlan and Hancock 2005, Clements et al. 2019, Gill et al. 2021). Herons are typically found in aquatic habitats, including lakes, marshes, rivers, forested streams, and along coastlines. Because they are relatively easy to observe, their behavior has been well studied, and they have been found to exhibit marked variation in foraging techniques, nesting, and migration (Meyerriecks 1960, Mock 1976, Green and Leberg 2005). Heron vagility is also notable, allowing them to move and proliferate across all continents, except Antarctica, and to colonize distant archipelagoes (e.g., Hawaii and Galapagos) that have never been connected to continental landmasses.

Currently, 5 subfamilies of herons are recognized (Kushlan and Hancock 2005): Tigrisomatinae (tiger herons), Botaurinae (bitterns), Ardeinae (typical herons), Agamiinae (Agami Heron Agamia agami), and Cochleariinae (Boat-billed Heron Cochlearius cochlearius). The Ardeinae is further divided into 3 tribes: Ardeini (day herons), Egrettini (egrets), and Nycticoracini (night herons). Although the monotypic Boat-billed Heron was originally placed in a separate family (Wetmore 1960), the monophyly of herons is now well established (Payne and Risley 1976, Sheldon 1987a, McCracken and Sheldon 1998, Sheldon et al. 2000, Chang et al. 2003, Huang et al. 2016). Nonetheless, substantial uncertainty regarding several aspects of heron phylogeny remains, including relationships between and within subfamilies and tribes, and the placement of several enigmatic taxa.

Herons are constrained to a wading, fishing body type, which is prone to convergent evolution, especially with respect to leg, neck, and bill length (Sheldon 1987a, McCracken and Sheldon 1998). They have also adapted to night feeding several times, and all night-feeding herons share distinctly similar features (e.g., relatively large eyes, broad bills, and short legs) whether they are closely related or not. Therefore, the use of morphology to reconstruct phylogenetic relationships has resulted in a series of phylogenetically inconsistent classifications. Bock (1956) was the last systematist to apply classic eclectic methodology (Mayr 1981) to the problem of heron relationships, using a combination of morphological and ecological traits to establish the composition of heron groups. He placed Cochlearius (a night-feeding heron) within the Nycticoracini, which included other night as well as some day herons (Gorsachius, Nycticorax, Nyctanassa, Pilherodius, and Syrigma), and he recognized 2 subfamilies, the Botaurinae and Ardeinae, the latter consisting of the tribes Tigriornithini (tiger herons), Ardeini, and Nycticoracini. Because herons have what appear to be consistent, phylogenetically influenced behavioral patterns, Curry-Lindahl (1971) tried using ethological characters to determine heron relationships. Payne and Risley (1976) were the first systematists to apply quantitative methods to the reconstruction of heron phylogeny. Using both cladistic and phenetic analysis methods, they compared osteological characters of 53 species and concluded that herons were best split into 4 subfamilies: Ardeinae, Nycticoracinae, Tigrisomatinae, and Botaurinae, with Ardeinae as sister to the rest. Cochlearius was placed within the Nycticoracinae but was elevated to tribal status (Cochlearini).

The first molecular study of heron relationships was by Sheldon (1987a). He employed DNA–DNA hybridization to compare the single-copy nuclear genomes of 27 species and found support for (1) the inclusion of the night herons Nycticorax and Nyctanassa within the Ardeinae, contrary to Payne and Risley (1976); (2) the paraphyly of Egretta and Ardea; (3) a sister relationship between the monotypic Whistling Heron (Syrigma sibilatrix) and Egretta; (4) a sister relationship between the monotypic Cattle Egret (Bubulcus ibis) and Ardea; (5) a sister relationship of Botaurinae and Ardeinae; and (6) placement of Cochleariinae and Tigrisomatinae as outgroups to the rest of the herons. While reconstructing the heron phylogeny, Sheldon (1987b) also discovered that different groups of herons evolved at different rates of sequence evolution; bitterns (Ixobrychus and Botaurus) evolved faster, and Boat-billed Heron and tiger herons (Tigrisoma) evolved more slowly than day and night herons. Subsequent DNA-DNA hybridization comparisons of the monotypic Zigzag Heron (Zebrilus undulatus) and White-crested Tiger Heron (Tigriornis leucolopha) (Sheldon et al. 1995; Figure 1A) placed Zebrilus within the Botaurinae and Tigriornis as a sister of the Neotropical tiger herons Tigrisoma. The precise position of Cochlearius, on the other hand, remained unresolved near the base of the heron tree. Despite these improvements, the DNA-DNA hybridization comparisons covered only about half the species in the heron family, leaving many relationships unresolved. Questions remained concerning the relationships of several enigmatic genera (e.g., Gorsachius, Agamia, Pilherodius, and Ardeola) and the identification of the sister clade of Ardeinae + Botaurinae.

FIGURE 1.

Previous phylogenetic hypotheses of Ardeidae, including (A) Sheldon et al. (1995), estimated using DNA-DNA hybridization distances; (B) Sheldon et al. (2000), estimated from a maximum likelihood analysis of one mitochondrial gene (cytochrome b); (C) Chang et al. (2003), estimated from a neighbor-joining analysis of one mitochondrial gene (12S rRNA); and (D) Huang et al. (2016), estimated from neighbor-joining analysis of one mitochondrial gene (cytochrome c oxidase subunit I). Trees were recreated as they appear in their respective manuscripts, while the classification follows that of Gill et al. (2021).

A quantitative comparison of Payne and Risley’s (1976) cladistic-osteological tree and the DNA-DNA hybridization tree (McCracken and Sheldon 1998) indicated that the osteological data employed by Payne and Risley (1976), particularly those pertaining to crania were homoplastic and prone to recovering ecological, not phylogenetic, relatedness. As a result, McCracken and Sheldon (1998) argued in favor of the DNA-DNA hybridization tree shown in Figure 1A as the better representation of heron phylogeny. They also noted that this tree was supported by vocal data (McCracken and Sheldon 1997). Subsequently, Sheldon et al. (2000) reconstructed the phylogeny using mtDNA cytochrome b (cytb) sequences of 15 species and 13 genera (Figure 1B). Their cytb tree was congruent with the DNA-DNA hybridization tree. In the process, they also discovered that cytb rates of evolution matched those of the single-copy nuclear genome, with bitterns again having faster rates of sequence evolution and Boat-billed Heron and tiger herons having slower rates of sequence evolution than typical herons.

At that time, although knowledge of heron phylogeny was substantially advanced, the existing data were still unable to determine the monophyly of the Nycticoracini, the relative positions of Tigrisomatinae and Cochlearius, and the relationships of several unexamined genera. Subsequent mtDNA studies using 12s rRNA or cytochrome oxidase I (COI) gene sequences (Chang et al. 2003, Huang et al. 2016; Figure 1C and D) provided insight into some of these relationships, particularly at the generic level. But, because the rates of evolution of these genes were, respectively, too slow and too fast, and the sequences were relatively short in length, the data could not resolve deeper relationships within the heron family. Despite improved taxonomic sampling (32 species and 17 genera), the phylogeny estimated by Huang et al. (2016) included several nodes that were not well resolved. In addition, the placement of Zebrilus as sister to the rest of the herons in their tree deviated substantially from previous phylogenetic estimates. Zhou et al. (2014, 2016) improved genetic sampling by sequencing and comparing complete mitochondrial genomes, but their limited taxonomic sampling left the positions of the Nycticoracini, the Ardeinae + Botaurinae clade, and certain genera (e.g., Agamia) unresolved. Additional studies have clarified some intrageneric relationships (Päckert et al. 2014, Tu et al. 2017, Duan et al. 2018), but have done little to provide a thorough picture of heron phylogeny or clarify some of the aforementioned outstanding questions.

Here, we sample genome-wide ultraconserved elements (UCEs) from all heron genera and >90% of species in the family to estimate heron phylogeny. UCEs are evolutionarily conserved markers that occur across eukaryotic lineages, permitting the reconstruction of phylogenies that span long timeframes (Faircloth et al. 2012, McCormack et al. 2012). UCEs have been useful for estimating phylogenies of several eukaryotic groups, including Hymenoptera (Faircloth et al. 2015) and Testudines (Crawford et al. 2015), and have been widely used in higher-level phylogenetic studies of birds (Moyle et al. 2016, White et al. 2017, Andersen et al. 2018, Oliveros et al. 2019, Salter et al. 2020). In the present work, we use UCEs to produce a significantly improved and more comprehensive estimate of heron phylogeny. We also use UCEs to reevaluate prior hypotheses with respect to among-lineage rate variation.

METHODS

Taxonomy

Throughout the manuscript, we follow the specific and generic classification of Gill et al. (2021). When referring to subfamilies and tribes we follow Kushlan and Hancock (2005).

Sampling and DNA Extraction

We sampled 55 species of herons and 3 Pelecaniformes outgroups (Hackett et al. 2008, Jarvis et al. 2014, Prum et al. 2015) (Table 1). For samples derived from muscle tissue or blood (hereafter referred to as “tissue samples”), we extracted DNA using the Qiagen DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA, USA), except for Ardea cinerea and Tigrisoma lineatum, which were extracted using the Maxwell RSC Blood DNA Kit (Promega, Fitchburg, WI, USA). Manufacturers’ protocols were followed in both instances. When extracting DNA from museum specimen toepads (collection date range: 1890–1989; Supplementary Material Table S4), we took additional measures to reduce contamination. In these instances, all subsampling was carried out under a laminar flow hood. DNA from toepads was extracted using the Maxwell RSC Blood DNA Kit (Promega, Fitchburg, WI, USA), following the manufacturer’s protocol. We quantified extracts of tissues using a Qubit 2.0 Fluorometer (Life Technologies) and a Qubit dsDNA BR Assay Kit following the manufacturer’s protocol. We quantified extractions of toepads using a Qubit dsDNA HS Assay Kit following the manufacturer’s protocol.

Table 1.

Genomic and mitochondrial summary statistics of samples. Species for which toepad samples were used are in bold. All non-toepad samples are derived from ethanol and/or liquid nitrogen preserved muscle tissues, with the exception of the Nycticorax nycticorax sample, which was derived from blood (but also preserved). Museum abbreviations are as follows: Louisiana State University Museum of Natural Science (LSUMNS), University of Kansas Natural History Museum (KUNHM), American Museum of Natural History (AMNH), Yale Peabody Museum (YPM), Museum of Southwestern Biology (MSB), University of Washington Burke Museum (UWBM), University of Alaska Museum (UAM), Field Museum of Natural History (FMNH), Denver Museum of Nature and Science (DMNS), San Diego Natural History Museum (SDNHM), Florida Museum of Natural History (FLMNH), University of Michigan Museum of Zoology (UMMZ), Bernice P. Bishop Museum (BPBM), National Museum of Natural History, Smithsonian Institution (USNM)

Taxon	Museum	Catalog no./Tissue no.	Locality	Number of cleaned reads	Mean trimmed read length	Number of UCEs	UCE base pairs	Average UCE contig length	Contigs > 1 kb	mtDNA reads	Average mtDNA coverage
Agamia agami	LSUMNS	B12815	Bolivia	1,004,098	142.7	4,782	3,653,798	764.1	409	-	-
Ardea alba	LSUMNS	B1343	Louisiana, USA	5,851,674	141.9	4,797	5,469,349	1140.2	3,871	14,313	155
Ardea cinerea	KUNHM	21788	Spain	5,544,480	143.7	4,782	5,768,590	1206.3	4,157	-	-
Ardea goliath	LSUMNS	B10361	Zoo/captive	5,355,170	143.7	4,776	5,363,246	1123.0	3,780	284	6
Ardea herodias	LSUMNS	B4095	Louisiana, USA	1,150,280	143.1	4,792	4,246,306	886.1	1,160	-	-
Ardea humbloti	AMNH	410700	Madagascar	1,561,798	123.7	1,866	574,304	307.8	0	322	6
Ardea insignis	YPM	041974	India	5,120,714	102.7	4,321	1,481,256	342.8	0	94	4
Ardea intermedia	MSB	177136	South Africa	4,225,662	140.6	4,826	5,335,651	1105.6	3,674	29,967	319
Ardea melanocephala	LSUMNS	B39300	Ghana	3,657,918	143.1	4,827	4,975,147	1030.7	3,066	1,391	19
Ardea pacifica	UWBM	62925	Australia	9,845,794	142.9	4,815	5,709,612	1,185.8	4,051	12,003	145
Ardea purpurea	LSUMNS	B39468	Ghana	8,174,222	142.7	4,781	5,665,514	1185.0	4,109	4,474	54
Ardeola bacchus	UAM	26000	Alaska, USA	5,189,426	143.0	4,814	5,590,940	1161.4	3,967	782	12
Ardeola idae	FMNH	282641	Mozambique	2,675,490	118.7	4,737	2,067,259	436.4	1	1,780	22
Ardeola ralloides	LSUMNS	B34283	South Africa	5,401,972	143.5	4,789	4,936,340	1030.8	3,058	1,829	23
Ardeola rufiventris	FMNH	282638	Mozambique	1,177,996	104.1	3,327	1,024,990	308.1	0	674	10
Ardeola speciosa	AMNH	DOT17256	Singapore	3,089,320	141.8	4,809	5,114,118	1063.4	3,336	8,533	99
Balaeniceps rex	LSUMNS	B19208	Zoo/captive	6,996,434	142.6	4,692	6,314,596	1345.8	4,378	-	-
Botaurus lentiginosus	LSUMNS	B18981	Louisiana, USA	5,395,564	143.4	4,801	5,085,418	1059.2	3,296	3,902	45
Botaurus poiciloptilus	UWBM	80401	Zoo/captive	4,993,778	141.4	4,788	4,736,274	989.2	2,602	-	-
Bubulcus ibis	LSUMNS	B19756	Louisiana, USA	4,210,744	141.7	4,823	5,285,277	1095.8	3,573	9,748	111
Butorides striata	LSUMNS	B12810	Bolivia	8,973,414	144.1	4,790	5,154,609	1,076.1	3,464	532	9
Butorides sundevalli	DMNS	36637	Ecuador	5,988,844	114.6	4,792	2,181,924	455.3	4	654	10
Butorides virescens	KUNHM	9507	El Salvador	11,855,224	141.6	4,772	5,515,136	1155.7	3,916	2,668	33
Cochlearius cochlearius	LSUMNS	B1339	Zoo/captive	4,534,388	143.6	4,817	5,090,094	1056.7	3,333	1,712	22
Egretta ardesiaca	SDNHM	51906	Zoo/captive	2,823,156	142.5	4,763	5,557,186	1166.7	3,829	-	-
Egretta caerulea	LSUMNS	B5283	Louisiana, USA	3,976,038	143.9	4,784	4,937,309	1,032.0	3,054	1,855	24
Egretta eulophotes	UAM	2805	Alaska, USA	2,820,300	80.4	3,659	1,177,743	321.9	1	135	4
Egretta garzetta	LSUMNS	B62605	Kuwait	8,184,010	143.8	4,807	5,221,173	1086.2	3,551	1,069	15
Egretta gularis	LSUMNS	B62603	Kuwait	7,979,534	143.6	4,825	5,325,234	1103.7	3,684	1,458	20
Egretta rufescens	LSUMNS	B6449	Louisiana, USA	4,618,372	142.8	4,784	5,425,143	1134.0	3,830	2,279	29
Egretta sacra	UAM	17951	Australia	2,767,046	142.1	4,730	5,981,141	1264.5	4,125	266	6
Egretta thula	LSUMNS	B6385	Louisiana, USA	6,613,380	142.3	4,799	5,380,174	1121.1	3,686	4,386	50
Egretta tricolor	LSUMNS	B19408	Louisiana, USA	8,833,596	143.6	4,832	5,638,645	1166.9	4,098	110	4
Egretta vinaceigula	LSUMNS	82383	Botswana	9,037,762	91.8	4,822	1,797,441	372.8	1	-	-
Gorsachius goisagi	USNM	572224	Indonesia	4,406,396	119.4	4,705	2,351,913	499.9	9	1,043	14
Gorsachius leuconotus	LSUMNS	B45084	Ghana	5,468,886	142.0	4,760	5,549,869	1165.9	3,956	1,929	25
Gorsachius melanolophus	KUNHM	10441	China	2,938,268	142.5	4,750	5,397,411	1136.3	3,674	947	13
Ixobrychus cinnamomeus	AMNH	DOT17237	Singapore	7,669,466	143.4	4,781	6,337,813	1325.6	4,384	48	3
Ixobrychus dubius	AMNH	DOT17848	Australia	4,142,108	143.1	4,771	5,099,404	1,105.7	3,280	107	4
Ixobrychus eurhythmus	AMNH	DOT17239	Singapore	1,986,654	142.9	4,773	6,368,615	1068.4	4,374	-	-
Ixobrychus exilis	LSUMNS	B3882	Louisiana, USA	5,947,416	141.9	4,737	5,508,084	1344.4	3,888	141	5
Ixobrychus flavicollis	UWBM	67898	Solomon Islands	6,057,482	141.1	4,830	5,336,042	1104.8	3,645	7,033	83
Ixobrychus involucris	LSUMNS	B35927	Trinidad and Tobago	3,592,406	142.6	4,781	5,275,144	1152.1	3,574	559	9
Ixobrychus sinensis	FLMNH	44361	Guam, USA	11,320,276	143.1	4,831	5,535,968	1145.9	4,006	-	-
Ixobrychus sturmii	UWBM	104503	Malawi	4,296,436	142.1	4,835	5,298,894	1095.9	3,610	5,361	64
Nyctanassa violacea	LSUMNS	B15549	Louisiana, USA	4,371,786	141.3	4,798	5,159,995	1075.4	3,357	1,142	16
Nycticorax caledonicus	KUNHM	10686	Australia	21,460,896	143.4	4,768	6,464,602	1355.8	4,434	1,603	20
Nycticorax nycticorax*	-	CHU006	Mozambique	5,027,066	142.9	4,735	5,651,243	1193.5	4,082	-	-
Pelecanus occidentalis	LSUMNS	B36186	Louisiana, USA	7,510,006	142.4	4,765	6,469,829	1357.8	4,421	147	5
Pilherodius pileatus	KUNHM	1247	Guyana	7,238,798	143.0	4,798	5,645,817	1176.7	4,032	4,431	54
Plegadis falcinellus	LSUMNS	B41207	Louisiana, USA	7,078,128	142.6	4,757	6,421,704	1349.9	4,457	34	3
Syrigma sibilatrix	LSUMNS	B6613	Bolivia	5,004,082	142.5	4,743	6,121,785	1290.7	4,272	742	12
Tigriornis leucolopha	UMMZ	235185	Gambia	6,735,734	141.9	4,831	5,269,006	1090.7	3,501	833	12
Tigrisoma fasciatum	LSUMNS	B4456	Peru	3,403,594	142.4	4,770	5,345,586	1120.7	3,724	1,605	20
Tigrisoma lineatum	KUNHM	3145	Paraguay	3,632,136	142.7	4,786	5,685,951	1188.0	4,014	55	4
Tigrisoma mexicanum	LSUMNS	B46531	Panama	5,860,178	139.9	4,813	5,320,783	1105.5	3,623	12,326	145
Zebrilus undulatus	LSUMNS	B12873	Bolivia	4,749,396	142.8	4,777	5,522,329	1156.0	3,959	756	11
Zonerodius heliosylus	BPBM	110742	Papua New Guinea	1,335,340	112.1	1,656	493,099	297.8	0	4,657	51

^*This sample is from a blood sample without an associated museum specimen that was sampled for Cumming et al. 2011. doi: 10.1007/s10393-011-0684-z. CHU006 is the sample identifier, and refers to Lake Chuali, Malawi.

Tissue and Toepad Sample Library Preparation

We standardized our non-toepad extractions by including 250 ng of DNA in 50 µL aliquots, which were subsequently fragmented with a Covaris S220 Focused-Ultrasonicator with the following settings: 175 W peak incident power, a duty factor of 2%, and 200 cycles per burst for 44–45 s. We targeted fragments of 500–600 base pairs (bp) in length. We then prepared libraries using Kapa Biosystems Library Hyper Prep Kits (#KK2602), following the manufacturer’s protocol, with some minor modifications: 25 µL per sample, at a concentration of 10 ng µL^–1, was subjected to end repair and A-tailing on a thermal cycler, followed by the ligation of 2 universal iTru stubs (iTru Stub Oligo 1 & iTru Stub Oligo 2) and incorporation of iTru dual-indexes (Glenn et al. 2019). Following ligation, we purified samples with Agencourt AMPure XP beads, at a volume of 0.8X. We then amplified libraries using the manufacturer’s recommended polymerase chain reaction (PCR) protocol, which included 10 cycles for the denaturation, annealing, and extension steps. We then performed an additional 1X Agencourt AMPure XP bead cleanup, followed by a quantification of the amplified libraries using a Qubit dsDNA HS Assay Kit on a Qubit 2.0 Fluorometer (Life Technologies). Pools of 8 amplified libraries were then grouped and enriched for 5,060 UCE loci using 5,742 baits in the myBaits UCE Tetrapods 5kv1-96 Library Capture Kit (sequences available at ultraconserved.org) (Faircloth et al. 2012), following the Mycroarray myBaits 3.01 Kit manufacturer’s protocol. After enrichment, we subsequently amplified capture-reaction products using 18 PCR cycles, followed by a 1.2X Agencourt AMPure XP magnetic bead cleanup. Finally, we quantified enriched, pooled libraries using a Qubit dsDNA HS Assay Kit on a Qubit 2.0 Fluorometer (Life Technologies). Prior to sequencing, pools consisting of equimolar libraries were combined.

We conducted toepad library preparation in a similar fashion to the tissue samples, with the following exceptions: (1) no sonication was performed due to the already degraded nature of the samples; (2) 1.5 mL lo-bind tubes were used in place of strip tubes when not conducting PCR; (3) 2.5 µM of iTru primer mix was added; (4) 3X Agencourt AMPure XP beads were used for initial, post-ligation, and post-amplification cleanups; (5) PCR amplifications of libraries were conducted for 12 cycles each; (6) 6 amplified libraries were added to each pool; (7) the hybridization temperature was set to 55°C; and (8) a Qiagen GeneRead Size Selection Kit (#180514), following manufacturer’s protocol, was used to remove adapter dimer and fragments <150 bp in size.

Sequencing

Two independent sequencing runs were conducted. The first run included 48 libraries that were pooled along with libraries from other projects and were sequenced on a lane of an Illumina HiSeq 3000 flow cell at the Oklahoma Medical Research Foundation (Oklahoma City, OK, USA), generating 150-bp paired-end reads. The second run included 10 dual-indexed samples that were pooled with samples from other projects on a single lane of an Illumina HiSeq 3000 flow cell at the same facility and also generated 150-bp paired-end reads.

Sequence Data Filtering, Ultraconserved Element Assembly, and Alignment

Summary characteristics of the alignments analyzed are provided in Supplementary Material Table S3, along with unique dataset identifiers. These identifiers will be used throughout the text to avoid ambiguity. Given that UCE alignments that contain toepad samples are likely to contain a high proportion of missing data, and because of the failure to enrich UCEs and to assemble full-length contigs from such samples (Hosner et al. 2016), we constructed the following UCE datasets: (1) containing UCEs enriched from all taxa (hereafter tissue + toepad datasets, datasets 1a–1g, 3a–3g), (2) containing UCEs enriched from muscle/blood samples that had been stored in ethanol and/or liquid nitrogen (hereafter tissue datasets, datasets 5 and 6), and (3) containing UCEs that had been “corrected,” while implementing the “correction” workflow in PHYLUCE (hereafter corrected tissue + toepad datasets, datasets 2a–2b, 4a–4b).

Reads were de-multiplexed by the Oklahoma Medical Research Foundation using bcl2fastq2 (Illumina). We trimmed low-quality bases and adapter sequences from raw reads using Illumiprocessor v2.0.9 (Faircloth 2013), which incorporates Trimmomatic (Bolger et al. 2014). This step conducted the following: (1) any read below 40 bp in length was dropped; (2) removed leading low quality or N bases were below a quality of 5; (3) removed trailing low quality or N bases were below a quality of 15; (4) scanned the read with a 4-base window, removing windows when the average quality score drops below 15; and (5) removed Illumina adapters from sequences. We followed the PHYLUCE Tutorial I guidelines (Faircloth 2015) and used several modules from the Python package PHYLUCE v1.6.6 (Faircloth 2016) for UCE processing and analysis. Contigs were assembled using SPAdes v3.12.0 (Bankevich et al. 2012). We used LASTZ v1.04 (Harris 2007) to extract contigs that matched UCE loci. After extracting UCE loci, we aligned them with MAFFT v7.407 (Katoh and Standley 2013). We did not trim nucleotides from the alignment ends during this step. For the tissue dataset, we trimmed alignments using Gblocks v0.91b (Castresana 2000), using default parameters, with the exception of the minimum number of sequences required for a flank position (b2), which was set at 65% of taxa. To explore if varying the b2 assembly parameters ameliorated phylogenetic artifacts typically associated with toepad samples (Moyle et al. 2016, Oliveros et al. 2019, Andersen et al. 2019), we constructed datasets from the tissue + toepad dataset while varying the following b2 values: 0.45, 0.55, 0.65, 0.75, 0.85, 0.95, and 1. This was done for both a complete matrix, that included UCE loci present in all taxa (datasets 3a–3g), and an incomplete matrix, that included UCE loci present in at least 75% of taxa (datasets 1a–1g).

To examine whether low-quality sites, not just missing data, were contributing to these phylogenetic artifacts, we also estimated trees from “corrected” alignments using the “mapping” and “correction” workflows as implemented in PHYLUCE v.1.7.1 (Faircloth et al. 2016). The mapping workflow maps reads to contigs, followed by the marking of duplicates and calculation of coverage. During the correction workflow, contig positions with quality scores <20, a depth of coverage <5, and the total number of alleles in called genotypes >2 were hard masked. We implemented this workflow and then continued with the standard analysis as described previously. To evaluate if phylogenetic artifacts were found in alignments with high-quality sites but varying amounts of missing data, we estimated 4 trees from concatenated alignments. These included (1) UCEs found in all samples and trimmed with a Gblocks b2 parameter = 0.85 (dataset 4a); (2) UCEs found in all samples and trimmed with a Gblocks b2 parameter = 1 (dataset 4b), (3) UCEs found in at least 75% of samples and trimmed with a Gblocks b2 parameter = 0.85 (dataset 2a), and (4) UCEs found in at least 75% of samples and trimmed with a Gblocks b2 parameter = 1 (dataset 2b).

Ultraconserved Element Trees

For each concatenated alignment we used RAxML v8.2.10 (Stamatakis 2014) to estimate a species tree, while assuming a general time-reversible (GTR) + Gamma model of molecular evolution. We carried out 20 maximum likelihood (ML) searches for the best-fit tree. We then generated nonparametric bootstrap replicates with the autoMRE function. Following the best-tree search and bootstrapping, we printed bootstrap values on the best-fitting tree. In an attempt to account for sources of gene-tree discordance driven by incomplete lineage sorting, we also inferred species trees with the coalescent-based approaches SVDQuartets (Chifman and Kubatko 2014) and ASTRAL III v5.6.3 (Zhang et al. 2018). SVDQuartets analyzes quartets of species using singular value decomposition of the matrix of site pattern frequencies and assembles a species tree from the resulting quartets and was implemented in PAUP* v4a159 (Swofford 2003). Quartet trees were combined into a full species tree by Quartet FM (Reaz et al. 2014). We estimated trees with SVDQuartets for the tissue (datasets 5 and 6a) and tissue + toepad (datasets 1c and 3c) datasets. In each instance, support was assessed with 100 nonparametric bootstrap replicates.

We estimated the best ML gene tree for each UCE locus in the incomplete tissue dataset (dataset 5) in RAxML v8.2.10 (Stamatakis 2014), assuming a GTR model of rate substitution and gamma-distributed rates among sites. Additionally, we generated 500 bootstrap replicates in RAxML v8.2.10 (Stamatakis 2014). We used these gene trees and bootstrap replicates as input for ASTRAL III v5.6.3 (Zhang et al. 2018). We assessed support for this phylogeny using 100 multi-locus bootstraps, with gene and site resampling. Multi-locus bootstrapping resamples sites within a locus and loci within a dataset (Seo 2008). Although ASTRAL is not strictly considered a coalescent method, it is statistically consistent with the multispecies coalescent model (Mirarab et al. 2014).

Gene Tree/Species Tree Discordance Analyses

High bootstrap support values are a hallmark of trees estimated from large concatenated datasets (e.g., Roycroft et al. 2020). They may provide inflated confidence and may occur even when the topology is incorrect (Kubatko and Degnan 2007, Liu and Edwards 2009). To provide a more detailed picture of the underlying variance in topologic support for nodes across the heron tree, we also evaluated the conflict between gene trees and different phylogenetic hypotheses.

Firstly, we evaluated the conflict between gene trees and a species tree estimated from the incomplete tissue dataset (dataset 5). Discordance was evaluated using PhyParts (Smith et al. 2015) and visualized with the script phypartspiecharts.py (https://github.com/mossmatters/MJPythonNotebooks). We used the previously estimated RAxML species tree and gene trees as input. Prior to running PhyParts, we rooted all trees (including the species tree) with the pxrr function from phyx (Brown et al. 2017). Gene trees that did not include any of the outgroup taxa were not included in the conflict analysis. As a result, we used a total of 4,748 gene trees as input for PhyParts (Smith et al. 2015).

Secondly, we evaluated gene tree/species tree discordance with IQTree v2.1.3 (Minh et al. 2020). Here, we quantified gene-concordance factors (gCF) and site-concordance factors (sCF) across nodes of trees resulting from incomplete and complete tissue datasets (datasets 5 and 6a). sCF and gCF factors are observed measures of variance in support, whereas bootstrap values are measures of the sampling variance of support (Minh et al. 2020). As a result, bootstrap values can be inflated in large datasets, even when sCF and/or gCF support is low. We first estimated the species and gene trees for both datasets with IQTree v2.1.3 (Minh et al. 2020). Subsequently, we performed a concordance factor analysis. We evaluated support for the species trees with 1,000 rapid bootstraps. The topologies of the resulting species trees were consistent with those found initially with RAxML.

Mitochondrial Genome Assembly and Alignment

We extracted and assembled off-target mitochondrial reads from genomic sequences of each individual that matched the mitochondrial genome of the Malayan Night Heron (Gorsachius melanolophus) (GenBank # KT364531.1) with MITObim v1.9.1 (Hahn et al. 2013) using the quick strategy and up to 50 iterations. MITObim is a Perl wrapper that employs the assembler MIRA v4.0.2 (Chevreux et al. 1999, Chevreux et al. 2004) to reconstruct mitochondrial genomes from raw data. In addition, we downloaded a mitochondrial genome from a White-eared Night Heron (G. magnificus) (GenBank # KT364529) and included it in this dataset. Sequences that had 65% or greater similarity to the G. magnificus mitochondrial genome for 13 mitochondrial coding regions (ND1, ND2, COX1, COX2, ATP8, ATP6, COX3, ND3, ND4L, ND4, ND5, Cytb, and ND6) across 49 taxa were extracted in Geneious Prime v2019.0.4. Each gene was aligned separately in Geneious, while implementing the default parameters of MAFFT v7.388 (Katoh and Standley 2013). Alignments were subsequently concatenated in Geneious. We examined alignments in Geneious and ensured that the sequences were contained within open reading frames. The resulting alignment included 49 taxa, 2 of which were outgroups (dataset 7). Mitochondrial gene sampling details are included in Supplementary Material Table S1.

Mitochondrial Genome Trees

We partitioned the mitochondrial alignment by codon position and performed a ML analysis using RAxML v8.2.11 (Stamatakis 2014), as implemented in Geneious v2019.0.4. For each codon, we used a GTR model of rate substitution and gamma-distributed rates among sites (GTR+G) and estimated 100 rapid bootstrap replicates for nodal support. We also performed a Bayesian inference (BI) analysis in BEAST v2.5.2 (Bouckaert et al. 2019) using the following models: GTR+G, relaxed log normal clock, and birth–death tree. Here, we used the default mutation rate of 1.0 substitutions per site. In addition, we included a prior forcing monophyly of the ingroup. We carried out two independent Markov chain Monte Carlo runs of 30 million generations each, sampling every 5,000 generations. We used LOGCOMBINER v2.5.2 (Bouckaert et al. 2019) to combine log and tree files, while discarding the first 25% of Markov chain Monte Carlo generations from each run as burn-in and resampling states every 10,000 generations. We assessed convergence and confirmed that the post-burn-in effective sample sizes of the tree likelihoods and parameters of the birth–death model were >200 by evaluating the combined log file in TRACER v1.7.1 (Rambaut et al. 2018). Finally, we used TREEANNOTATOR v2.5.2 (Bouckaert et al. 2019) to summarize the remaining trees as a maximum clade credibility tree. To examine whether the mtDNA topology was robust to genotype quality, we removed taxa with ≤10x mean coverage (dataset 8), and performed the same phylogenetic analyses as described above.

Tests of Molecular Evolution

Given previous phylogenetic studies demonstrating variation in rates of evolution among heron lineages (Sheldon 1987b, Sheldon et al. 2000), we compared rates across branches in the UCE tree by implementing a two-cluster (Takezaki) test with the program LINTRE (Takezaki et al. 1995). The two-cluster test examines whether the average substitution rate of two clades separated by a given node is equivalent. Using the complete tissue (dataset 6a) UCE RAxML tree as the input topology, we computed pairwise distances via the Tamura-Nei gamma option in LINTRE. Differences in delta (δ) and deviations in sister branches from zero (Z) were calculated for each node in the tree by LINTRE.

In addition, we estimated branch-wise substitution rates using the uncorrelated relaxed clock model, as implemented in BEAST v2.5.2 (Bouckaert et al. 2019). Because of computational constraints, we only included a total of 500 UCEs in this analysis. These UCEs were drawn from the complete tissue dataset (dataset 6a). Here, we created 5 unique sets of 100 UCEs each (datasets 6b–6f), which were randomly selected without replacement and concatenated with the AMAS alignment tool (Borowiec 2016). For each set, we conducted at least 2 independent Markov chain Monte Carlo runs of 40 million generations each, sampling every 5,000 generations. For the first set, we carried out 3 independent runs, and for the rest, we carried out 2. We unlinked site models across UCE loci while linking the clock and tree models across loci. We used a HKY (Hasegawa et al. 1985) model of substitution for each locus, and an uncorrelated relaxed clock model and birth–death tree model. Default priors for these models were used except for the following modifications. We set the mean of the ucld.mean parameter to follow a normal distribution prior with a mean of 0.0005 substitutions per site per million years and a standard deviation of 0.0001, which was estimated from a phylogeny of birds (Prum et al. 2015, Berv and Field 2018). We also included the RAxML topology of the concatenated complete tissue dataset (dataset 6a) as a starting tree, along with a prior that forced the monophyly of the outgroup taxa. We used LOGCOMBINER v2.5.2 (Bouckaert et al. 2019) to combine log and tree files, while discarding the first 10% of Markov chain Monte Carlo generations from each run as burn-in and resampling every 10,000 states. We assessed the convergence of parameter estimates across runs and confirmed that the post-burn-in effective sample sizes for most parameters across runs were >200 by evaluating the resulting combined log file in TRACER v1.7.1 (Rambaut et al. 2018). The only exception was the mrca.age parameter of the outgroup, which had effective sample size values below 200 for sets 1 and 2, at 181 and 179, respectively. We used TREEANNOTATOR v2.5.2 (Bouckaert et al. 2019) to summarize the remaining trees as a maximum clade credibility tree and visualized branch-wise substitution rates in FigTree v1.4.4 (Rambaut 2009).

RESULTS

Ultraconserved Element Recovery

Following the trimming of low-quality reads and removal of adapters, we kept an average of 2.9 million paired reads from the tissue samples, ranging from 502,049 to 10,730,448 reads, with an average read length of 142.6 bp (Table 1). From the toepad samples, we obtained an average of 1.89 million paired reads, ranging from 588,998 to 4,518,881, and an average of 107.5 bp in length (Table 1).

We recovered an average of 57,007 (range: 8,651–265,132) contigs assembled from reads from tissue samples (Supplementary Material Figure S9), with an average length of 436 bp and an average depth of 12.5x. We obtained an average of 77,145 (range: 13,611–213,423) contigs from reads from toepad samples (Supplementary Material Figure S9), with an average length of 263 bp and an average depth of 5.60x.

From the tissue sample contigs, we enriched an average of 4,787 UCEs, ranging from 4,692 to 4,835, which were an average length of 1,139 bp and an average depth of 68x (Table 1, Supplementary Material Figure S10). From the toepad sample contigs, we enriched an average of 3,765 UCEs, ranging from 1,656 to 4,822, with an average length of 371 bp and an average depth of 28x (Table 1, Supplementary Material Figure S10).

Ultraconserved Elements and Mitochondrial Alignments

For the dataset deriving only from tissue samples, the complete alignment consisted of 3,681 UCE loci with a length of 4,032,900 characters (dataset 6a). The incomplete alignment (loci found in at least 75% of all taxa) of the same dataset consisted of 4,773 UCE loci with a length of 5,113,334 characters (dataset 5). For the tissue + toepad dataset, the complete alignments consisted of 466 UCE loci with lengths of 99,432– 558,073 characters (datasets 3a–3g). The incomplete alignments consisted of 4,695–4,756 UCE loci with lengths of 1,056,772–5,203,332 characters (datasets 1a–1g). For the corrected dataset, the complete alignments consisted of 77–79 UCE loci with lengths of 7,727–48,900 characters (datasets 4a–4b). The incomplete alignments (loci found in at least 75% of all taxa) of the corrected dataset consisted of 4,695–4,701 UCE loci with lengths of 689,798–3,054,227 characters (datasets 2a–2b).

The average depth of coverage of the mitochondrial contigs was 38x, ranging from 3 to 319 (Table 1). After gene extraction and alignment, we used 2 data matrices that included alignments from 13 protein-coding mtDNA genes: One including 49 taxa with a length of 11,579 characters (dataset 7) and the other including 37 taxa (≥10x coverage) also 11,579 characters in length (dataset 8).

Ultraconserved Element Trees

ML and coalescent analyses of both the incomplete and complete tissue datasets (datasets 5 and 6a) produced a mostly congruent and well-supported phylogenetic hypothesis (Figure 2; RAxML tree of dataset 5). The only taxa whose positions were equivocal were Egretta thula and Gorsachius melanolophus. Egretta thula was resolved by SVDQuartets in both analyses (Supplementary Material Figures S2 and S3) as sister to E. sacra + E. ardesiaca with modest bootstrap support (76% and 82% for the complete and incomplete datasets, respectively). Alternatively, both RaxML trees (Figure 2 and Supplementary Material Figure S1) and the ASTRAL tree (Supplementary Material Figure S4) placed E. thula as sister to E. gularis + E. garzetta, with 100% bootstrap support. Gorsachius melanolophus appeared as a sister of the rest of Ardeinae in all trees, except for the RAxML tree of the complete tissue dataset (dataset 6a; Supplementary Material Figure S1). In that case, G. melanolophus was a sister of a clade consisting of Ardea, Bubulcus, Butorides, Ardeola, Nycticorax, and Nyctanassa, albeit with low bootstrap support (56%). All other nodes were congruent across trees and well-supported (>85% bootstrap support). These analyses found each of the 5 subfamilies to be monophyletic, with the Tigriornithinae sister to the rest of the herons. Cochleariinae is sister to Agamiinae + Botaurinae + Ardeinae, and Botaurinae + Ardeinae are sister to each other. This result provides stability in a region of the heron tree that changed frequently in previous classifications.

FIGURE 2.

Maximum likelihood tree of the incomplete tissue dataset (dataset 5), estimated in RAxML. All nodes have 100% bootstrap support unless otherwise labeled. Illustrations by Bennu Birdy. The classification of taxa in this figure coincides with the one recommended in the text.

ML and SVDQuartets analyses of the incomplete and complete tissue + toepad datasets (datasets 5 and 6a) produced trees that were inconsistent with each other and had comparatively poor nodal support (Supplementary Material Figures S5–S8). Support was particularly poor for deeper nodes in trees estimated from the complete dataset (Supplementary Material Figures S5 and S7). Toepad samples were prone to long-terminal branch lengths, which were pulled towards the root of the tree and most closely allied to other toepad samples, often with high bootstrap support (Supplementary Material Figures S5 and S6). These phylogenetic artifacts were not discernibly improved by how alignment ends were trimmed by Gblocks. Given that ASTRAL is a gene-tree reconciliation method and is expected to be negatively affected by samples with short contig lengths (Hosner et al. 2016, Moyle et al. 2016), we decided to forgo phylogeny estimation in ASTRAL using matrices from this dataset.

ML analyses of the corrected incomplete and complete tissue + toepad datasets (datasets 2a–2b, 4a–4b) resulted in trees that resolved some of the phylogenetic artifacts, but no one tree resolved these issues entirely (Supplementary Material Figures S16). For example, all of these trees were still plagued by toepad samples with long-terminal branch lengths, suggesting the masking approach employed was perhaps too conservative.

Mitochondrial DNA Trees

ML and BI trees estimated from 13 coding regions extracted from 49 mtDNA genome assemblies (Figure 3) were generally well resolved and congruent with the UCE trees derived from the tissue dataset. The mtDNA ML tree, however, did have some nodes that were poorly supported (<70% bootstrap). When evaluating nodes that were either modestly or strongly supported in the mtDNA trees, only the placement of the clade consisting of Gorsachius melanolophus + G. goisagi and that consisting of Nyctanassa violacea + Nycticorax caledonicus were equivocal. In the ML tree, Nyctanassa and Nycticorax formed the sister group of the rest of the Ardeinae, whereas, in the BI tree, they were the sister to a clade consisting of G. goisagi + G. melanolophus. Both ML and BI trees of the reduced (37 taxa; Supplementary Material Figures S14 and S15) mtDNA dataset largely recapitulated the topology found in the expanded mtDNA dataset (49 taxa; Figure 3). The only exception was the relationship between Egretta rufescens and E. caerulea in ML trees. These were found to belong in a clade with Egretta tricolor in the expanded mtDNA tree. However, the reduced mtDNA tree did not match this topology, but rather placed E. caerulea as sister to the clade containing E. thula, E. garzetta, and E. gularis, although this grouping had modest bootstrap support (71%; Supplementary Material Figure S14).

FIGURE 3.

Maximum likelihood (A) and Bayesian inference (B) trees of the mitochondrial dataset (dataset 7). All nodes have 100% bootstrap or 1.0 posterior probability support unless otherwise specified. Nodes with bootstrap/posterior probability support <70% have collapsed. The Tigrisomatinae (brown), Cochleariinae (green), Botaurinae (yellow), and Ardeinae (blue) are highlighted. Figure 3 is continued on the next page.

Gene Tree/Species Tree Discordance

We also experienced a substantial amount of gen-tree discordance in nodes of the heron tree, including those previously reconstructed with high bootstrap support (Figure 4). For example, in the PhyParts analysis of the incomplete dataset, the average concordance in the tree (percentage of gene trees that supported the shown topology) was 51.6% (range: 9.03–95.1, Supplementary Material Table S5). Four branches within Ardeinae were particularly plagued by gene-tree discordance and corresponded with the 4 lowest gCF and sCF scores (Supplementary Material Figures S12 and S13). These included branches placing: (1) Egrettini as sister to the Ardeini + Nycticoracini, (2) Gorsachius melanolophus as either sister to the rest of the Ardeinae (Supplementary Material Figure S13) or to Egrettini + Nycticoracini (Supplementary Material Figure S12), (3) Nycticoracini as sister to the Ardeini, and (4) Egretta thula as sister to E. gularis + E. garzetta.

FIGURE 4.

Patterns of gene tree discordance in the heron phylogeny. Pie charts indicate degree of gene discordance at each node. Percentage of concordant gene trees are indicated in blue, discordant gene trees for the top alternative bipartition in green, all other genes supporting conflicting bipartitions in red, and uninformative gene trees in gray. Each branch is labeled with the number of genes in support (top) and in conflict (bottom) with the given clade.

Tests of Molecular Evolution

The Takezaki rate test produced the asynchronous tree shown in Figure 5. The delta (δ) values (absolute value of the difference between distances of clades to the outgroup), deviations from zero (Z) at each node, and other data produced by LINTRE are provided in Supplementary Material Table S2. The 4 largest differences, all with δ values >0.02, are: (1) a slowdown in the rate of Tigrisoma tiger herons relative to other species; and (2) a speedup in bitterns (Ixobrychus and Botaurus) relative to the rest of the herons; a (3) slowdown in Cochlearius relative to the clade consisting of Agamia, Ardeinae, and Botaurinae; and (4) a speedup in Gorsachius melanolophus relative to the clade consisting of Ardea, Butorides, Ardeola, Nycticorax, and Nyctanassa. Slow rates were also observed for Zebrilus and Agamia.

FIGURE 5.

Tree showing branch length differences among heron lineages based on a LINTRE analysis (Takezaki et al. 1995). Black circles at nodes indicate the degree of rate change in adjacent branches (larger = greater) and are proportional to the delta values. Blue branches indicate a rate slowdown, red a speedup. Data pertaining to individual branches and nodes are provided in Supplementary Material Table S2. To aid with interpretation select delta values and associated circle sizes are displayed.

Branch-wise substitution rates estimated from the relaxed clock model analysis in BEAST produced results that were broadly consistent with the Takezaki rate test results (Figure 6; Supplementary Material Figures S17–S20). In particular, higher substitution rates, in all 5 sets, were noted for internal branches leading to the bitterns, with the branch leading to Ixobrychus + Botaurus often exhibiting the fastest rate across the tree. In addition, we noted comparatively slower rates on the branches leading to Tigrisoma, Agamia, and Cochlearius.

FIGURE 6.

Branch-wise substitution rates in the heron tree. Units are in substitutions per site per million years. Depicted here is a set of 100 unique and randomly chosen UCEs from the complete tissue dataset (dataset 6b; Set 1). Branch widths and color are proportional to the median of the molecular rate. Red colors indicate faster rates; blue colors indicate slower rates. Branch-wise substitution rates of Sets 2–5 are depicted in Supplementary Material Figures S17–S20.

DISCUSSION

We present a well-resolved estimate of heron phylogeny based on a dataset including >90% of heron species and two independent sources of molecular data: UCEs and mtDNA. Although molecular phylogenetic inquiries have vastly improved our understanding of heron phylogeny over the last 35 years, they have failed to (1) estimate consistent relationships among the major clades of herons; (2) place several enigmatic, iconic heron genera (e.g., Agamia, Gorsachius, and Zebrilus); and (3) undertake any comparisons of Zonerodius. These shortcomings cast substantial uncertainty on the evolutionary history of herons and have led to fluctuation in their classification (Kushlan and Hancock 2005, Clements et al. 2019, Gill et al. 2021). Taxonomists are now positioned to develop a stable, accurate classification based on phylogeny and to employ that phylogeny for studies of heron ecology, behavior, biogeography, and molecular evolution.

Relationships Among Major Groups

Our estimate of heron phylogeny (Figure 2) offers a clear picture of the composition and position of the main heron groups. The tiger herons sensu stricto (Tigriornis and Tigrisoma) form the sister group of the rest of the herons, and the Cochleariinae (Cochlearius) and Agamiinae (Agamia) are successive sister groups of the Botaurinae (Ixobrychus and Botaurus) + Ardeinae (Gorsachius, Syrigma, Pilherodius, Nycticorax, Nyctanassa, Ardeola, Butorides, Egretta, Bubulcus, and Ardea). The relationships among major groups are consistent between and well supported in both UCE and mtDNA trees. The branching pattern between major groups largely recapitulates the mtDNA trees of Sheldon et al. (2000) and Huang et al. (2016), but notably disagrees with Huang et al. (2016) in the placement of Zebrilus. Our results place Zebrilus in the Botaurinae, concordant with Sheldon et al. (1995) and Päckert et al. (2014). To some degree, the composition of the major groups differs from the classification of Kushlan and Hancock (2005). For example, the mtDNA trees and corrected UCE trees provide evidence that Zonerodius is a member of the Ardeinae, not the Tigrisomatinae. Also, the Nycticoracini sensuKushlan and Hancock (2005) include Nycticorax and Gorsachius, but not Nyctanassa, which is placed in a separate tribe (Egrettini). We found Nycticorax and Nyctanassa to be sisters and Gorsachius to be polyphyletic, with G. goisagi and G. melanolophus in a different clade than G. leuconotus and G. magnificus. In a previous mtDNA study, G. goisagi and G. melanolophus formed the sister group of Nycticorax, although this relationship was only strongly supported in one of the analyses (Zhou et al. 2016).

Composition and Placement of Enigmatic Taxa

Although some taxa, including Agamia, Bubulcus, Gorsachius, Ixobrychus involucris, I. exilis, and Zebrilus, were included in previous molecular tree-building efforts, their relationships have remained ambiguous and their classification, therefore, not solidly founded. Zonerodius was wholly lacking in previous molecular comparisons; this study is the first to include this odd Papua New Guinean species and place it objectively in a phylogenetic tree.

Agamia agami (Agami Heron).

All analyses of the UCE data placed Agamia agami as a sister to the Ardeinae + Botaurinae with high support. This position disagrees with the hypothesis that Agamia is most closely allied with the Ardeinae (Bock 1956, Payne and Risley 1976). Huang et al. (2016) found Agamia to be the sister of Ardeola, although this relationship was not strongly supported. Kushlan and Hancock (2005) placed Agamia agami in its own subfamily, Agamiinae. Our results agree that Agamia lies outside of the Ardeinae and is best considered a monotypic subfamily.

Bubulcus ibis (Cattle Egret).

The relationships of Bubulcus have long been debated. Bock (1956) thought it was close to Ardeola, but all molecular studies indicate it belongs in or as a sister to Ardea (Sheldon 1987a, Chang et al. 2003, Zhou et al. 2014). Our mitochondrial analyses show Bubulcus to be the sister of Ardea (Figure 3), whereas UCE analyses all indicate that Bubulcus is embedded within Ardea (Figure 2, Supplementary Material Figures S1–S4). DNA-DNA hybridization comparisons (with a limited sampling of Ardea species) suggested that Bubulcus was the sister of a clade consisting of Ardea herodias, A. cocoi, A. sumatrana, A. melanocephala, Casmerodius albus (= A. alba), and Egretta intermedia (= A. intermedia) (Sheldon 1987a). The mtDNA data of Zhou et al. (2014) placed Bubulcus as sister to a clade consisting of Ardea cinerea, A. purpurea, A. intermedia, and A. modesta [alba], and the 12S rRNA analyses of Chang et al. (2003) placed Bubulcus as sister to a group that included Ardea cinerea, A. purpurea, Egretta alba (= A. alba), and E. intermedia (= A. intermedia). However, a more-broadly sampled COI tree by Huang et al. (2016) disagreed with this placement and suggested that Bubulcus is the sister of Ardea alba + A. intermedia. While the exact placement of Bubulcus with respect to Ardea remains in doubt, we advocate the inclusion of Bubulcus in Ardea because of the clear close relationship of the two taxa. Interestingly, although we have known for 35 years that Bubulcus lies next to or within Ardea, its behavioral and morphological distinctiveness have prevented a generic name change in the major bird checklists (Kushlan and Hancock 2005, Clements et al. 2019, Gill et al. 2021).

Gorsachius (Old World Forest Night herons).

We found Gorsachius sensu Gill et al. (2021) to be polyphyletic. MtDNA data placed G. melanolophus + G. goisagi and G. leuconotus + G. magnificus as distinct clades (Figure 3) in different parts of the heron tree. Gorsachius melanolophus was equivocally placed by the UCE data as either sister to the rest of the Ardeinae (Figure 2, Supplementary Material Figures S2–S4) or as sister of a clade consisting of Ardea + night herons + (Ardeola + Butorides), with low bootstrap support (Supplementary Material Figure S1). This ambiguity is reflected in mtDNA trees, which place G. melanolophus + G. goisagi as sister of the night herons (Figure 3B), or as sister of Ardea + night herons + pond herons (Figure 3A). Although our mtDNA BI analysis suggested the clade of G. melanolophus + G. goisagi is sister to the night herons, this result was not found in any of our other analyses. Diminished gene/site discordance of G. melanolophus as sister to the rest of the Ardeinae provides greater confidence in this placement (Supplementary Material Figures 12 and 13). Gorsachius leuconotus occurred as sister to the Egrettini in all UCE trees (Figure 2; Supplementary Material Figures S1–S4). In mtDNA trees, the clade comprising G. leuconotus + G. magnificus likewise was sister to the Egrettini. In general, these findings support the results obtained by Zhou et al. (2016), who found G. goisagi and G. melanolophus are sister taxa and G. magnificus is the sister of Egretta. Our study is the first to use molecular data to examine and discover the sister relationship between G. leuconotus and G. magnificus, which were both previously placed in the genus Calherodius (Peters 1931). For this reason, we recommend resurrecting Calherodius for these taxa.

Ixobrychus exilis (Least Bittern).

In most of our trees, Ixobrychus exilis appears as the sister of Botaurus lentiginosus and B. poiciloptilus with high support. The one exception is the ML analysis of mtDNA data, which failed to resolve the position of I. exilis relative to the Botaurinae. The placement of I. exilis as sister of B. lentiginosus + B. poiciloptilus corroborates previous studies that have suggested I. exilis is more closely allied to Botaurus than to Ixobrychus (Sheldon 1987a, Päckert et al. 2014).

Ixobrychus involucris (Stripe-backed Bittern).

All UCE analyses indicated that Ixobrychus involucris is the sister of a clade of bitterns that includes I. exilis + Botaurus with high support. BI analysis of mtDNA data corroborated these UCE results, but ML analysis of mtDNA did not resolve the position of I. involucris relative to Botaurus. BI analyses of two mtDNA fragments by Päckert et al. (2014) provided equivocal support for the placement of I. involucris; one analysis corroborated our UCE and BI mtDNA results, whereas the other placed I. involucris as sister of other members of Ixobrychus. Both relationships found by Päckert et al. (2014) received poor posterior probability support. Huang et al. (2016) found that I. involucris is more closely related to Botaurus pinnatus + B. stellaris than to 5 members of Ixobrychus (sinensis, minutus, flavicollis, cinnamomeus, and eurhythmus). In general, most analyses and data agree with the hypothesis that I. involucris is more closely allied to members of Botaurus (including I. exilis) than to all other members of Ixobrychus.

Zebrilus undulatus (Zigzag Heron).

All trees consistently place Zebrilus as sister of the Botaurinae (bitterns) with high support. Zebrilus has been considered close either to tiger herons or bitterns because it shares characteristics with both groups, especially a barred, cryptic plumage. Bock (1956) placed Zebrilus within the tiger heron tribe Tigriornithini, based on shared ecological and plumage traits. Payne and Risley (1976), however, considered it to be a monotypic tribe, Zebrilini, within the Botaurinae. This arrangement was supported by DNA hybridization and mitochondrial COI data (Sheldon et al. 1995, Päckert et al. 2014). It disagreed, however, with the COI tree presented by Huang et al. (2016), which placed Zebrilus as sister of rest of the herons.

Zonerodius heliosylus (Forest Bittern).

Our findings vis-à-vis Zonerodius are, unfortunately, equivocal because of bias imposed by the poor quality of the UCE data, which were generated from a toepad sample. A relatively small number of UCE sequences were recovered, and they were generally short in length (Table 1, Supplementary Material Figures S5 and S6). RAxML and SVDQuartets trees from the tissue + toepad dataset all placed Zonerodius at the base of the heron tree with high support. We found in these trees that Zonerodius commonly clustered with other toepad taxa and was invariably found as sister to Ardea humbloti. In contrast, the corrected trees provided alternative topologies, and did not recover Zonerodius at the base of the heron tree (Supplementary Material Figure S16). Instead, 3 of the 4 trees here placed Zonerodius as embedded within Ardeola. These results support the mtDNA analyses, which placed it as sister to members of Ardeola (Figure 3).

Zonerodius has long been considered a member of Tigriornithinae, based on its forest habitat and skull morphology (Bock 1956, Payne and Risley 1976). However, Zonerodius never appeared close to members of the Tigriornithinae in any of our trees, suggesting that it is convergent to tiger herons in these traits. Further comparisons either including higher-quality genomic data or the sampling of more than one individual will be required before Zonerodius is to be definitively located within the heron phylogeny. However, we recommend it be considered as closely allied to Ardeola based on our mtDNA comparisons.

Bootstrap Values Obscure Gene Tree/Species Tree Discordance

Several nodes of the heron tree, even those with high bootstrap support, were plagued by gene tree discordance. Not surprisingly, the greatest discordance was found at nodes separated by short branches and characterized by instability. In particular, the sister relationships of Gorsachius melanolophus, Egretta thula, and Nycticoracini were inconsistent in UCE and mtDNA trees. Despite a large amount of discordance (Figure 4, Supplementary Material Figures S12 and S13), the percentage of gene trees, as recovered by the PhyParts analysis of the incomplete tissue dataset, that supported the most common conflicting bipartition at a node was generally low (mean: 0.83, range: 0.23–3.07). A large percentage of gene trees supported discordant bipartitions, excluding the most common bipartition (mean: 46.7, range: 4.61–87.8). What these results demonstrate is that while discordance was common, few nodes had a strong signal for any one particular alternative topology, and that discordance consisted largely of many low-frequency conflicting topologies. The two nodes with the highest proportion of gene trees that supported the most common conflicting bipartition were the relationships between Egrettini and the Ardeini + Nycticoracini (3.07%) and of Gorsachius melanolophus as sister to the rest of the Ardeinae (3.05%). The high proportion may explain why G. melanolophus was ambiguously placed between the incomplete and complete datasets, and it indicates that there is some phylogenetic signal for the two topologies in the nuclear genome. Of these two topologies, either sister to the rest of the Ardeinae or to a clade consisting of Ardea, Bubulcus, Butorides, Ardeola, Nycticorax, and Nyctanassa, the relationship of G. melanolophus as sister to the rest of the Ardeinae had proportionally less discordance across genes and sites (Supplementary Material Figures S12 and S13) in the incomplete dataset (gCF: 11.3, sCF: 36.3; Supplementary Material Figure S13) vs. the complete dataset (gCF: 8.72, sCF: 32.4; Supplementary Material Figure S12), providing a greater degree of confidence in the relationship of G. melanolophus as sister to the rest of Ardeinae.

Toepad Phylogenetic Artifacts Are Influenced by Missing Data and Site Quality

Among trees, we noticed three anomalous phylogenetic patterns that were consistently associated with toepad samples: (1) unusually long-terminal branch lengths, (2) long-branch attraction, often between individuals from distinct clades, and (3) a tendency to branch near the root of the tree. This last pattern manifested itself idiosyncratically; some toepad individuals appeared as sister to small clades (e.g., a genus) and others occurred as sister to larger clades (e.g., most other herons). These anomalous patterns are commonly associated with historical samples (Moyle et al. 2016, Oliveros et al. 2019, Andersen et al. 2019, Salter et al. 2022), suggesting they are caused by data biases to which toepads are prone. They clearly do not reflect phylogenetic signals.

Phylogenetic artifacts such as long-terminal branches, long-branch attraction, and erroneous topologies can be driven by poor data quality and missing data, or a combination of the two. Artificially long-terminal branches, which often lead to long-branch attraction and incorrect topologies, can occur in trees constructed from poorly aligned data (Hossain et al. 2015) or with systematically missing data (Lemmon et al. 2009, Simmons 2012, Darriba et al. 2016). Missing data can also cause taxa to be pulled toward the root (Hosner et al. 2016, Moyle et al. 2016). This interaction between poor quality and missing data can make it difficult to diagnose the driver of a particular anomalous phylogenetic pattern a posteriori. We found it challenging to ascribe particular patterns in our trees to one or a combination of these phenomena, but note that some of our results suggest that data quality was an important driver of long-branch attraction. While long-branch attraction occurred in all of our tissue and toepad RAxML trees, we found that only one of the four trees estimated after the correction-workflow demonstrated long-branch attraction between toepad samples (Supplementary Material Figure S16). The other 3 trees were the first to “break up” the attraction between Ardea humbloti and Zonerodius heliosylus. Instead, these 3 trees indicated that Zonerodius heliosylus is allied with Ardeola, supporting our mtDNA results.

Prior to implementing the correction workflow, the sister relationship of Ardea humbloti and Zonerodius heliosylus proved recalcitrant. It occurred in all UCE trees, irrespective of the amount of missing data. This suggests that the clustering of these two taxa is not driven primarily by missing data. We hypothesize, rather, that this pattern is mostly influenced by shared single nucleotide polymorphisms introduced during the de novo assembly of contigs implemented in the standard PHYLUCE workflow. This workflow ignores some variants at heterozygous positions, with the most numerous variants being called and the less numerous alternative being discarded (Iqbal et al. 2012). Downstream, the effect is that low-coverage sites are more likely to introduce erroneous genotype calls into the alignment. This issue is exacerbated when using historical samples, which often contain fewer reads. Indeed, 2 of the 3 lowest-ranking toepads, in terms of cleaned reads, were Ardea humbloti and Zonerodius heliosylus (Table 1). In particular, we hypothesize that the single nucleotide polymorphisms shared between these two taxa might have resulted from PCR errors caused by toepad-specific patterns of DNA damage, such as cytosine deamination (Sefc et al. 2007). In light of these findings, we recommend researchers who compare historical samples to estimate phylogeny consider not only the effects of missing data, but also of site quality (Smith et al. 2020).

Ultraconserved Elements Corroborate Differences in Molecular Rates Among Lineages

The early use of genomic data to reconstruct heron phylogeny played a key role in the discovery of variable rates of genomic evolution, not only in birds but in all organisms. The first molecular reconstruction of heron phylogeny (Sheldon 1987a) employed DNA–DNA hybridization, a method that compared entire single-copy genomes among species. At that time, it was known that rates varied among some genes (e.g., protein genes; Brownell 1983), and between nuclear and mitochondrial DNA (Brown et al. 1979). Moreover, nuclear genomes of major groups, such as hominoids and rodents, were suspected to evolve at different rates (Britten 1986), but there was no direct evidence (i.e., without reliance on vague fossil or biogeographic dates) that rates of genomic evolution varied among different groups of organisms. Indeed, resistance to such an idea was strong (Wilson et al. 1977, Sarich and Cronin 1980, and Sibley and Ahlquist 1984). The heron rate-study (Sheldon 1987b) was the first to demonstrate directly (using a genetic distance-based ratio test) not only that genomic rates varied among groups of organisms but that they varied among lineages as closely related as species within a single family of birds. Subsequently, Sheldon et al. (2000) demonstrated that variability in mtDNA lineage-based rate patterns matched those of nuclear DNA in herons, indicating a previously unknown rate relationship between the two genomes.

Our UCE comparisons validate these early single-copy nuclear and mitochondrial DNA rate discoveries. All 3 datasets show that bitterns (Ixobrychus and Botaurus) had faster rates of sequence evolution than other herons and that tiger herons (Tigrisoma), and Boat-billed Heron had slower rates of sequence evolution. In addition, some of the newly added taxa have been found to contain either faster (Gorsachius melanolophus) or slower (Zebrilus and Agamia) rates of sequence evolution than their sister groups. The cause, or combination of causes, for such variation in rates of sequence evolution, however, remains uncertain. The search for explanations of molecular rate variation has inspired many suggestions, including differences between taxa in metabolic rate, generation time, rate of germ cell division, body temperature, DNA repair efficiency, population size, and clade size (Martin and Palumbi 1993, Sheldon and Bledsoe 1993, Rand 1994, Omland 1997, Bleiweiss 1998, and Gillooly et al. 2005).

Conclusion

Using thousands of UCEs and 13 coding regions of mtDNA, we produced a relatively well-resolved phylogenetic tree for the herons. Unfortunately, UCE trees that included toepad samples were incongruent with trees generated from samples that had been preserved (stored in liquid nitrogen and/or ethanol). In addition, trees with toepad samples contained phylogenetic artifacts, which we suspected were driven by missing data and low sequence quality. However, we consistently reconstructed, using ML, BI, and coalescent phylogenetic approaches, well-supported trees when toepad samples were excluded. In addition, toepad-free UCE trees were largely concordant with our mtDNA trees, which placed several toepad-based taxa (e.g., Ardeola idae, Ardea humbloti, Zonerodius heliosylus) within expected clades. These results support our hypothesis that the anomalous topologies recovered when the toepad data were included in the UCE tree building were artifacts, not true biological groupings.

In addition to providing a well-supported picture of the relationships among subfamilies and tribes, we found that the genera Ardea, Gorsachius, and Ixobrychus are not monophyletic, contra to current classification (e.g., Gill et al. 2021). We also suggest that the enigmatic Zonerodius heliosylus, for which we provide the first molecular data, is not a tiger heron (Tigrisomatinae) as previously thought (Kushlan and Hancock 2005), but is best considered as either a member of or closely allied to the genus Ardeola. Future work, in addition to confirming our hypothesis regarding Zonerodius with improved sampling, should focus on clarifying taxonomic issues at the species level, particularly in species with high subspecific diversity. A thorough sampling of the Ardea intermedia, Butorides virescens/striata, and Egretta thula/gularis/garzetta complexes, for example, would help to clarify outstanding taxonomic questions within these groups (Kushlan and Hancock 2005).

Based on our phylogenetic analyses, we recommend the following classification of herons. Taxonomic changes we recommend are denoted with a *.

Family Ardeidae Leach 1820

Subfamily Tigriornithinae Bock 1956
Tigrisoma Swainson 1827
Tigrisoma lineatum
Tigrisoma mexicanum
Tigrisoma fasciatum
Tigriornis Sharpe 1895
Tigriornis leucolopha
Subfamily Cochleariinae Chenu and Des Murs 1854
Cochlearius Brisson 1760
Cochlearius cochlearius
Subfamily AgamiinaeKushlan and Hancock 2005
Agamia Reichenbach 1853
Agamia agami
Subfamily Botaurinae Reichenbach 1849-50
Zebrilus Bonaparte 1855
Zebrilus undulatus
Botaurus Stephens 1819
Botaurus lentiginosus
Botaurus poiciloptilus
Botaurus stellaris
Botaurus pinnatus
Botaurus exilis*
Botaurus involucris*
Ixobrychus Billberg 1828
Ixobrychus sturmii
Ixobrychus sinensis
Ixobrychus minutus
Ixobrychus flavicollis
Ixobrychus cinnamomeus
Ixobrychus eurhythmus
Ixobrychus dubius
Subfamily Ardeinae Leach 1820
Gorsachius Bonaparte 1855
Gorsachius goisagi
Gorsachius melanolophus
Calherodius Bonaparte 1855
Calherodius leuconotus*
Calherodius magnificus *
Syrigma Ridgway 1878
Syrigma sibilatrix
Pilherodius Reichenbach 1853
Pilherodius pileatus
Egretta Forster 1817
Egretta picata
Egretta novaehollandiae
Egretta rufescens
Egretta tricolor
Egretta caerulea
Egretta ardesiaca
Egretta sacra
Egretta thula
Egretta garzetta
Egretta gularis
Egretta vinaceigula
Egretta dimorpha
Egretta eulophotes
Nycticorax Forster 1817
Nycticorax nycticorax
Nycticorax caledonicus
Nyctanassa Stejneger 1887
Nyctanassa violacea
Zonerodius Salvadori 1882
Zonerodius heliosylus
Ardeola Boie 1822
Ardeola ralloides
Ardeola grayii
Ardeola speciosa
Ardeola bacchus
Ardeola idae
Ardeola rufiventris
Butorides Blyth 1852
Butorides virescens
Butorides sundevalli
Butorides striata
Ardea Linnaeus 1758
Ardea cinerea
Ardea herodias
Ardea cocoi
Ardea pacifica
Ardea melanocephala
Ardea purpurea
Ardea goliath
Ardea alba
Ardea ibis*
Ardea intermedia
Ardea insignis
Ardea sumatrana
Ardea humbloti

Acknowledgments

We would like to thank the following collection managers, staff, and curators for their assistance with the processing of requests and shipping of samples: Paul Sweet, Peter Capaniolo, Tom Trombone, American Museum of Natural History; Ben Marks, Field Museum of Natural History; Donna Dittmann, Louisiana State University Museum of Natural Science; Mark Robbins, Rob Moyle, Town Peterson, University of Kansas Natural History Museum; Garth Spellman, Denver Museum of Nature and Science; Phil Unitt, Kevin Burns, San Diego Natural History Museum; Chris Milensky, Smithsonian Natural Museum of Natural History; Molly Hagemann, Bernice Pauahi Bishop Museum; Kristof Zyskowski, Yale Peabody Museum; Kevin Winker, John Withrow, University of Alaska Museum. We also like to thank the following institutions for providing samples to Kevin McCracken, Fred Sheldon, and Philip Lavretsky: Academy of Natural Sciences of Philadelphia, University of Washington Burke Museum, Museum of Southwestern Biology, Peabody Museum of Yale University, University of Michigan Museum of Zoology, University of Florida Natural History Museum. We thank Paul Hime (University of Kansas Biodiversity Institute) and Brant Faircloth (Louisiana State University) for assistance with several bioinformatic roadblocks. We’d like to thank members of the Moyle Lab: Luke Campillo, Luke Klicka, and Lucas DeCicco, for helpful discussions and James Kushlan for comments on the manuscript. We are indebted to Sarah Maclean (https://www.bennubirdy.com), who produced the remarkable illustrations. Lastly, we thank reviewers and editors for providing thoughtful and constructive feedback.

Funding statement

This work was funded in part by NSF grant DEB-1557053 to R.G.M, a Frank M. Chapman Memorial Fund Grant from the American Museum of Natural History to J.P.H, and the James A. Kushlan Endowment for Waterbird Biology & Conservation at the University of Miami to K.G.M. Portions of the analyses were conducted on the high-performance cluster at the University of Kansas and at Texas Tech University.

Ethics statement

Samples included in this study were derived from museum specimens and were not subject to Institutional Animal Care and Use Committee (IACUC) approval.

Author contributions

(1) J.P.H and R.G.M conceived of the project; (2) J.P.H and J.H. generated the data; (3) J.P.H, J.H., S.S., and C.O. analyzed the data; (4) J.P.H., R.G.M., F.H.S., and K.G.M. wrote the manuscript with feedback from all the remaining authors; (5) K.G.M, P.L., F.H.S., J.P.H., and R.G.M. contributed substantial funding and/or materials.

Data availability

All raw sequencing reads have been uploaded to NCBI’s Sequence Read Archive (SRA), accessioned under BioProject PRJNA658323. PHYLUCE code used in this study can be accessed at https://github.com/faircloth-lab/phyluce. Contigs, alignments, gene trees, and xml files have been deposited on Dryad and are available from Hruska et al. (2023). Data from: Ultraconserved elements resolve the phylogeny and corroborate patterns of molecular rate variation in herons (Aves: Ardeidae). Ornithology 140:ukac000. https://doi.org/10.5061/dryad.0gb5mkm5c. Custom code used in this study can be accessed at https://github.com/jphruska/Ardeidae.