Artifex, Kallal & Hormiga, 2018

† Caerostris View in CoL morphology information arbitrarily appended to this taxon for total evidence analysis. †† Morphological coding supplemented by Ambicodamus sp. following Álvarez-Padilla et al. (2009).

methods of resampling were used to evaluate nodal support. Morphological data for Caerostris from Kuntner et al. (2008) was arbitrarily assigned to sequence data associated with Caerostris sp. 1243; the taxon ‘Nicodamidae’ is a composite following Álvarez-Padilla et al. (2009). One of the conspecific sequenced Phonognatha , Deliochus and Artifex gen. nov., was arbitrarily assigned morphological and behavioural data. Morphological and behavioural character state changes were be mapped on to the total evidence tree using WINCLADA ( Nixon, 1999) and ambiguous optimizations were examined using both ACCTRAN and DELTRAN.

MODEL-BASED ANALYSES

Specimens used in this study are presented in the Supporting Information (Appendix S2), including locality and voucher information. Specimens preserved in 95% ethanol were used for DNA extraction using the Qiagen DNEasy kit. Left legs were used for extractions and the remainders are preserved as vouchers. Six markers were amplified for analyses. These are mitochondrial ribosomal markers 12S rRNA (~400 bp) and 16S rRNA (~550 bp), cytoplasmic ribosomal markers 18S rRNA (~1800 bp) and 28S rRNA (~2700 bp), nuclear protein-coding gene histone H3 (~327 bp), and mitochondrial protein-coding gene cytochrome c oxidase subunit I or COI (~800 bp). PCR was completed using the Promega GoTaq kit, using the primers in Table 2. The generalized thermocycle was initial denaturation at 94 °C for 2 minutes, followed by a cycle of denaturation at 94 °C for 30 s, annealing at 40–56 °C for 35 s (see details following) and elongation at 65 °C for 30 s repeated 34 times, with a final elongation step at 72 °C for 3 min, followed by cooldown for 10 °C for 30 min, then held at 4 °C. Annealing temperatures were started at the following temperatures, then modified up or down by 2 °C to achieve clean, thin bands on a 1.5% agarose gel: 12S (46–48 °C), 16S (46– 48 °C), 18S (48–52 °), 28S (48–52 °C), COI (40–48 °C) and H3 (53–56 °C). Amplified products were sent to Macrogen USA in Rockville, MD for sequencing. Contigs were formed using GENEIOUS 6.0.6 (http:// www.geneious.com; Kearse et al., 2012), then queried against NCBI BLAST nucleotide database to check for contamination.

Multiple sequence alignments were completed using MAFFT v.7.017 ( Katoh & Standley, 2013). Alignments of 12S and 16S were completed using L-INS-i, ideal for samples with a single conserved domain. For 18S and 28s, the E-INS-i method was used, which is used to account for missing fragments. The protein-coding genes COI and H3 were aligned using MACSE ( Ranwez et al., 2011), which can detect frameshifts and stop codons in alignments. To account for missing data and poor alignment, TRIMAL v.1.2 was used with the gappyout setting ( Capella-Gutierrez, Silla-Martinez & Gabaldon, 2009). The concatenated molecular sequence matrix consisted of 6028 bp.

The maximum likelihood analysis was conducted using RAxML 8.2.8 ( Stamatakis, 2014) on CIPRES ( Miller, Pfeiffer & Schwartz, 2011). A total of 93 taxa with nucleotide sequence data were included. Partition schemes were tested using PARTITIONFINDER 1.1.1 ( Lanfear et al., 2012), with ten possible partitions: four non-protein coding markers (12S, 16S, 18S, 28S) and two protein-coding markers (COI, H3) further partitioned by codon position. The rapid bootstrapping algorithm option to find a single best-scoring tree with the GTRGAMMA model was used. Bootstrap iterations were set to 1000. The root was set as Uloborus glomosus ( Walckenaer, 1841) ( Uloboridae ).

Relaxed clock Bayesian analyses were conducted using MRBAYES 3.2.6 ( Ronquist & Huelsenbeck, 2003) on the high-performance cluster Colonial One at The George Washington University. A molecular matrix of 93 taxa and a total evidence matrix, including 95 taxa (two additional Phonognatha that could not be sequenced), were analysed. As in the maximum likelihood analysis, PARTITIONFINDER 1.1.1 ( Lanfear et al., 2012) was used on the molecular partitions; the morphological partition was analysed using the Mk model ( Lewis, 2001). We used the fossilized birthdeath prior for under-sampled lineages and the independent gamma rates model with broad speciation priors (exp[10]), extinction (beta[1,1]) and fossilization (beta[1,1]) following Zhang et al. (2015) and Pyron (2017). The tree age prior is set to 175–200 Mya based on analyses by Dimitrov et al. (2017). The clock rate prior is set by log normalizing the estimated substitution rate of cytochrome c oxidase I ( Bidegaray-Batista & Arnedo, 2011). A total of 16 chains (4 cold, 12 heated) were run for 100 million generations, with the first 25% discarded as burn-in. Convergence was considered achieved when estimated sample sizes (ESS) were above 100 and traces from log files examined in TRACER v.1.6 ( Rambaut et al. 2014) were plateaued. The early fossil record of many orb-weaving spider clades is sparse, and interpreting the fossils described so far is difficult due to the quality of preservation. There are two fossils currently described as early araneids: Mesozygiella dunlopi Penney & Ortuño, 2006 and Olindarachne martinsnetoi ( Downen, 2011), both from the late Aptian Age (115–121 Myr). Nevertheless, the placement of these taxa as araneids seems somewhat tenuous, at least in part due to the synapomorphies allowing attribution to Araneidae being very limited (e.g. Coddington, 1986; Dimitrov et al., 2017). Therefore, despite the strengths of additional calibrations, we have opted to calibrate solely based on the substitution rate, rather than introduce uncertainty from incorrect fossil placement. Trees were rooted as in maximum likelihood analyses.

BIOGEOGRAPHY

The R package BIOGEOBEARS was used for ancestral range estimation as it includes a framework for comparison of various models of biogeography ( Matzke, 2013). This method allows for modelling using the LAGRANGE DEC model ( Ree & Smith, 2008) and maximum-likelihood approximations of DIVA ( Ronquist, 1997) and BAYAREALIKE ( Landis et al., 2013). These biogeographic models have various strengths and weaknesses: DEC has parameters for dispersal and extinction that may change at cladogenetic events ( Ree & Smith, 2008; Matzke, 2014); DIVALIKE approximates the parsimony effects of DIVA ( Ronquist, 1997), namely dispersal and vicariance; BAYAREALIKE ( Landis et al., 2013), which is suited to broad sympatry. Furthermore, BIOGEOBEARS implements another free parameter, j, which is designed for ‘jump’ speciation, or a founder effect from colonization from a main population. The j parameter can be applied to each of the three biogeographic models, allowing for six possible configurations for comparison that were applied to the time-calibrated total evidence results from MRBAYES.

Seven areas were allocated for the 20 zygielline taxa included in this study. These are Palearctic (P), Nearctic (N), Indomalayan (I), New Caledonia (C), Queensland north of the St. Lawrence Gap (Q), eastern Australia south of the St. Lawrence Gap (S) and Western Australia (W). The first three areas represent biogeographic realms, and the latter four are more specific to Australia and New Caledonia, and have been informed by prior studies of biogeography (including spiders) in the area ( Crisp et al., 2004; Rix & Harvey, 2012). Queensland north of the St. Lawrence Gap (Q) represents tropical wet forests that sometimes have a distinct fauna from more temperate forests to the south (S) and is thought to be biologically relevant since the early Tertiary Period. The appearance of the Nullarbor Plain is thought to occur during the early Miocene, and could be a viciarance event separating Western Australia (W) from the south-east (S). The 20 zygielline taxa were coded as present or absent for each of these areas based on material examined for the taxonomic revision. Additionally, a dispersal matrix approximating relative probability of movement between two areas was also applied, with intercontinental dispersal coded as less likely (0.1) than intracontinental dispersal (1.0). The nested models (with and without the founding parameter) were compared using likelihood ratio tests. AIC scores were calculated for comparison of the six models.

RETREAT EVOLUTION

We coded the 89 unique species in the total evidence data matrix based on the retreat form used by that species. The terminal taxa were coded in one of the following five character states based on our own observations or literature records: no snare web (e.g. arkyids, mimetids), snare web and no retreat (e.g. tetragnathids), web with adjacent retreat (typically composed of leaves modified with silk, e.g. Araneus marmoreus , Fig. 3E), web with integrated leaf retreat [e.g. Phonognatha melania ( Fig. 2A), Phonognatha graeffei ( Fig. 3A), Acusilas sp. ( Fig. 3B)], and web with integrated detritus retreat (e.g. Metepeira labyrinthea and Spilasma duodecimguttata , Fig. 3C, D). This coding is a modification of the analysis of Gregorič et al. (2015), who broadly classified all relatively complex retreats as silk tubes. We divided the ‘silk tube’ class further, and decided to include relatively simple leaf-curling as a retreat adjacent to the web, to more finely examine the variation of retreats within Araneidae . It must be emphasized that retreat type can vary within genera, and the coding is for the specifically sampled taxon, where possible.

The ultrametric tree from the total evidence analysis was analysed in R (R Core Team, 2015) using the packages APE ( Paradis, Claude & Strimmer, 2004), PHYTOOLS ( Revell, 2012) and GEIGER ( Harmon et al., 2007). Three discrete character likelihood models were applied: equal transition rates between all states (ER), equal symmetric transition rates between two states that differ from a different relationship between two states (SYM) and different rates for each state, where each rate is parameterized separately (ARD). This was carried out using the rerootingMethod and fitMk functions in phytools, and results were scored using AIC. The best scoring model was used to perform an ancestral state reconstruction on web retreats using stochastic character mapping using the phytools’ function make.simmap. The time-calibrated total evidence tree from the Bayesian analysis was overlaid with 1000 stochastic character maps from the posterior distribution (nsim = 1000, q = ‘mcmc’), allowing a probability at each node for the likelihood of the five retreat codings to be represented by a pie chart.