SDL homologous sequence analysis

2.4.1 The homologous sequences of RIFA and outgroup species

In this study, the sequences of the Taiwan population and the population of the original habitat, South America, were analyzed. There are 25 sequences from the South American population, all from northern Argentina (A. Cohanim and E. Privman, personal communication). Since the coverage of these sequences was different, the analysis of different regions would not include all of them but the sequences which cover the analyzed homologous region. There were 10 sequences from the Taiwan population representing 10 different CSD allele sequences and all of them were include in every analysis (Huang et al., unpublished).

Monomorium pharaonis was used as the out-group species in this study. This is the species which is categorized in the same tribe, Solenopsidini, which includes TFA and RIFA. Additionally, the published genome of this species is the closest available one to Solenopsis genus (Mikheyev and Linksvayer, 2015) (GenBank assembly accession:

GCA_000980195.3). The choice of this out-group species was to minimize possible long branch attraction which may result in random rooting (Kinene et al., 2016). For genomic regions lacking data (sub-fragments, CoP, HVR and UFO) unrooted trees were constructed.

2.4.2 DNA extraction

In order to get better sequence information, higher quality DNA was required. For this reason I modified the Puregene® DNA purification Kit (Gentra systems). In this part, all samples, include alates and workers, were disrupted and lysed as above, then the DNA was extracted with the kit. The DNA pellet was dissolved in 100 μl of DNA hydration solution and then stored at 4⁰C.

2.4.3 Primer design

The master trigger gene for sex determination may be located in the hypervariable region. For convenience, we have divided it into three regions (CoP, HVR and UFO, see Fig 3) based on 2 conserved sequences found separating these fragments in RIFA (Huang et al., unpublished). Based on an alignment of 18 RIFA and 2 TFA hypervariable region sequences currently available, multiple primer pairs were designed to these conserved sequences to amplify each of the CoP, HVR, and UFO sub-regions. EGF-like is a gene adjacent to the hypervariable region. RIFA data showed the DNA sequence of the gene is slightly polymorphic associated with the different haplotype of the hypervariable region.

Despite EGF-like being a conserved gene, since it is close to the hypervariable region, primer pairs for this gene were also designed.

2.4.4 Selection of candidate neutral loci

First, 7 contigs from the TFA PacBio genome assembly were selected randomly.

Then I checked whether any >5kb region within each contig lacked any predicted open reading frame (ORF) >300 bp using the orffinder tool (https://www.ncbi.nlm.nih.gov/orffinder/). The candidate loci were also examined for any potential long coding RNAs by mapping to RIFA RNAseq data (Huang et al., unpublished) from each larval instar and sex. I chose a 1kb region satisfying these criteria on five of

these contigs as the neutral control fragments for this study.

2.4.5. Cloning and sequencing

2.4.5.1 Specific fragment PCR

The EGF-like gene was amplified using the Kapa HiFi HotStart PCR Kit (Kapa Biosystems). The mix was 30 μl, consisting of 1 μl 0.19~159 ng/μl DNA, 6 μl 5x κ HiFi Buffer (Fidelity), 0.9 μl dNTPs (10 mM each dNTP), 0.5 μl κ HiFi HotStart DNA Polymerase (2.5 U/μl) and 0.9μl forward and reverse primers (10 mM/μl), and water for the remaining volume. PCR reactions consisted of an initial denaturation temperature of 95⁰C for 3 min, followed by 35 normal cycles and a final extension at 72⁰C for 30 min.

The normal cycles were: denaturation at 98⁰C for 20 sec, annealing at 58⁰C for 15 sec and extension at 72⁰C for 1.5 min. PCR reactions were carried out in an ABI 9700 thermal cycler (Applied Biosystems).

The conditions for the neutral control PCRs was the same as above, except from a premix containing both forward and reverse primers (10 mM /μl) was used. The PCR reactions consisted of an initial denaturation temperature of 95⁰C for 4 min, followed by 35 normal cycles and a final extension at 72⁰C for 10 min. The normal cycles were:

denaturation at 98⁰C for 30 sec, annealing at 60⁰C for 30 sec and extension at 72⁰C for 5min. PCR reactions were carried out in a SuperCycler Trinity (Kyratec.).

The target DNA was amplified using the Kapa Long Range HotStart PCR Kit (Kapa Biosystems). The mix was 30 μl, consisting of 1 μl 0.19~159 ng/μl DNA, 6 μl 5x κ LongRang buffer (without Mg²⁺), 0.9 μl dNTPs (10mM each dNTP), 2.1 μl MgCL2 (25 mM), 0.5 μl κ LongRang HotStart DNA Polymerase (2.5 U/μl) and 0.5 μl forward and reverse primers (10 mM /μl), and water for the remaining volume. PCR reactions consisted of an initial denaturation temperature of 94⁰C for 4 min, followed by 35 normal

cycles and a final extension at 72⁰C 7 min. The normal cycles were: denaturation at 94⁰C for 20 sec, annealing at 50⁰C for 15 sec and extension at 72⁰C for 7 min. PCR reactions were carried out as above PCR machine.

2.4.5.2 Transformation, culturing, and sequencing

The EGF-like gene and neutral control fragment PCR products were separated by 1% TAE agarose gel electrophoresis at 60 V, 220 mA for 2 hr. Then, the gels with DNA were cut and DNA was purified using the Viogene® Gel / PCR DNA isolation systems.

The DNA fragment was inserted into the pCRTM TOPO II blunt vector (Zero Blunt®

TOPO® PCR Cloning Kit) and then plasmid was transformed into the E. coli competent cell DH5α (EGF-like gene was transformed into Fast-TransTM competent cells and neutral control fragments were transformed into RBC HIT competent cells). After 37⁰C plate culture overnight the transformed colonies were confirmed by colony PCR. Positive colonies were used to inoculate a liquid culture overnight. Plasmid DNA was extracted using the QIAprep mini kit (Qiagen system) and then sent to the company Genomics BioSci&Tech for Sanger sequencing.

Cloning of the target (CoP, HVR, and UFO) DNA PCR products was similar as above with differences as follows. First, all transformations were with RBC HIT E. coli competent cell DH5α. Second overnight plating was at 30⁰C, and then culture plates were incubated in 37⁰C 4~6 hr prior to inoculating colonies into a liquid culture overnight.

Third, all plasmids were checked for successful inserts by EcoRI restriction enzyme digestion for 2 hr in 37⁰C followed by gel electrophoresis. Only clones with inserts were sent to Genomics BioSci&Tech for Sanger sequencing.

2.4.6 Analysis of nucleotide diversity and molecular evolution

All TFA sequences (this study) and RIFA sequences were combined and aligned by

Multiple Alignment using Fast Fourier Transform program (MAFFT) (Katoh and Standley, 2013). Then I used the software Jmodeltest (Guindon and Gascuel, 2003) to determine the best substitution model under the Akaike information criterion. After, all the sequences were analyzed to make the phylogenic tree with 10,000 bootstrap replications in best model through the software MEGA 7 (Kumar et al., 2016). For the nucleotide substitution models of the hypervariable region, gamma distributed General Time Reversible model (CoP, HVR) and gamma distributed with invariant site General Time Reversible model (UFO, EGF-like) were used. For the neutral regions, gamma distributed General Time Reversible model (LG1N), Hasegawa-Kishino-Yano model (LG3N, LG5N1 and LG10N) and uniformed General Time Reversible model (LG5N4) were used. The partial deletion option and the BioNJ initial tree were used on all trees.

All trees were created with the online tool Interactive Tree Of Life (iTOL version 3.5.4) (Letunic and Bork, 2016).

Finally, the nucleotide diversity (π) was calculated by averaging the nucleotide difference per site between two randomly picked sequences extracted respectively from these two species or from one species. The Tajima’s D statistic was calculated and show the p-value according to the coalescent model. These statistics were calculated using the DnaSP (Rozas et al., 2003) software.