Overexpression of a novel candidate oncogene KIF14 correlates with tumor progression and poor prognosis in prostate cancer

Prostate cancer (PCa) is the second leading cause of death from cancer in men. The mechanism underlying tumorigenesis and development of PCa is largely unknown. Here, we identified Kinesin family member 14 (KIF14) as a novel candidate oncogene in PCa. We found that KIF14 was overexpressed in multiple PCa cell lines and primary PCa tissues. Knockdown of KIF14 in DU145 and PC3 prostate cancer cells suppressed cell proliferation, induced cell cycle arrest and apoptosis. Transcriptome analysis by RNA-sequencing demonstrated that KIF4 suppression led to transcriptional changes of genes involved in p53 and TGF-beta signaling pathway. In addition, upregulated expression of GADD45A, GADD45B, p21, PIDD and Shisa5, which contribute to growth arrest and apoptosis induction, and downregulated CCNB1 that promotes cell cycle progression were confirmed by quantitative real-time PCR after KIF4 knockdown. We further found that KIF14 protein level was positively correlated with T stage and Gleason Score. Patients with higher KIF14 expression had shorter overall survival time than those with lower KIF14 expression. Thus, our data indicate that KIF14 could act as a potential oncogene that contributes to tumor progression and poor prognosis in PCa, which may represent a novel and useful prognostic biomarker for PCa.

INTRODUCTION
Prostate cancer (PCa) is one of the most commonly diagnosed malignant tumor and the second leading cause of cancer death among men in the United States. In 2016, about 180,890 new PCa cases were diagnosed and 26,120 PCa deaths were occurred in USA [1]. About one-third of patients with organ-confined PCa fail surgical treatment and progress to aggressive or metastatic lethal disease in 10 years [2]. Thus, understanding mechanisms underlying the tumorigenesis and development of PCa is urgently needed for early diagnosis, risk stratification and selection Research Paper www.impactjournals.com/oncotarget of therapy strategies against PCa. However, despite the impressive progress in discovery of genomic alterations in PCa, novel genomic bio-markers and proteins still remain to be identified.
The kinesin family proteins (KIFs) are ATP-and microtubule-associated motor proteins involved in cell division, microtubule polymer dynamics, intracellular transportation and signal transduction [3]. Kinesin family member 14 (KIF14) contains a C-terminal citron kinase binding region and an N-terminal protein regulating cytokinesis binding region [4,5], and was found to bind to chromatin and microtubules during the formation of the bipolar spindle [4,6]. Elevated KIF14 expression can induce rapid and error-prone mitosis, leading to aneuploidy during tumorigenesis [6]. Although the exact molecular functions of KIF14 are still under investigation, increasing evidences has implicated KIF14 as an oncogene whose overexpression has been found in multiple types of tumors, including hepatocellular carcinoma, lung cancer, breast cancer, glioma, laryngeal carcinoma and ovarian cancer, etc [6][7][8][9][10][11][12]. Moreover, KIF14 levels are prognostic for the outcome of hepatocellular carcinoma and ovarian cancer, where KIF14 overexpression enhances tumor growth, while its knockdown decreases tumorigenicity in vitro and in xenografts [13]. However, the expression of KIF14 and its role in prostate carcinoma has not been reported thus far.
In this study, we found that KIF14 was upregulated in human PCa cell lines and primary PCa tissues. Knockdown of KIF14 in DU145 and PC3 prostate cancer cells reduced cell proliferation, induced cell cycle arrest and apoptosis. Transcriptome analysis by RNAsequencing demonstrated that KIF4 suppression led to transcriptional changes of genes involved in p53 and TGF-beta signaling pathway. In addition, gene expression of a few of cell cycle regulated proteins (GADD45A, GADD45B, p21 and CCNB1) and apoptotic proteins (PIDD and SHISA5) were confirmed by quantitative real-time PCR. We further found that KIF14 expression was significantly upregulated in primary human PCa specimens. KIF14 protein level was positively correlated with clinical stage and Gleason Score. Patients with higher KIF14 expression had shorter 5-year overall survival time than those with lower KIF14 expression. Thus, our data indicate that KIF14 overexpression contributes to prostate tumor progression and is associated with poor prognosis for PCa patients, which could be served as a novel and useful prognostic indicator for PCa.

KIF14 was overexpressed in PCa
Previous studies suggested that KIF14 was increased in multiple cancer including hepatocellular carcinoma (HCC), lung cancer, breast cancer, glioma, laryngeal carcinoma and ovarian cancer [6][7][8][9][10][11][12]. However, the role of KIF14 in prostate cancer remains unknown yet. We first checked the gene expression of KIF14 in a serial PCa cell lines and normal prostate tissue and found that KIF14 was significantly upregulated in PCa cells as compared to that in normal prostate tissue ( Figure 1A). We further performed a summary for KIF14 gene expression extracted from to the Oncomine database, a cancer microarray database and Integrated Data-Mining Platform [14]. A fold-change of 2 for KIF14 gene expression compared to control was required for criterion. Based on analysis on two datasets derived from samples of patients with PCa [15,16], we determined that KIF14 is highly overexpressed in prostate carcinoma as compared with normal prostate tissues ( Figure 1B). Thus, the overexpression of KIF14 in multiple PCa cell lines and primary PCa suggested that KIF14 might be a potential oncogene in PCa.

Inhibition of KIF14 leads to G2 arrest and reduced proliferation in PCa cells
To explore the possible mechanism of KIF14 as an oncogene, we adopted a small interfering RNA (siRNA) approach to knockdown KIF14 and examined the effect of KIF14 on cell cycle distribution and proliferation. In the experimental set-up, cell lines of DU145 and PC3, where KIF14 was markedly upregulated as compared to normal control, were chosen for the functional study. The knockdown efficiency was first validated after 48 h transfection with siRNA. Quantitative real-time PCR analysis showed that the mRNA levels of KIF14 were reduced by 89 ± 1% (P < 0.01) and 91 ± 5% (P < 0.01) in DU145 and PC3, respectively ( Figure 2A). We examined whether reduced KIF14 expression would affect prostate cancer cell growth by counting the cell numbers for three constitutive days after siRNA transfection. DU145 and PC3 cells exhibited reduced proliferation after KIF14 expression was knocked down ( Figure 2B and Figure 2C). Flow cytometric analysis of propidium iodide (PI)-stained cells revealed that cells accumulated and arrested in G2 as evidenced by a significant increase in the fraction of cells in G2 in siKIF14-trasnfected DU145 and PC3 cells as compared to that in scramble control transfected cells ( Figure 2D and Figure 2E), indicating that loss of KIF14 induces G2 arrest in PCa.

Knockdown of KIF14 leads to apoptosis of PCa cells
We further examined the effect of KIF14 on apoptosis. Flow cytometric analysis of Annexin V-PE/7-AAD-stained cells after 48 hours after siRNA transfection discovered that loss of KIF14 produced more apoptotic cells of DU145 and PC3 ( Figure 3A-3C). Taken together, these data indicate that KIF14 overexpression www.impactjournals.com/oncotarget  is necessary for G2/M cell cycle progression, promotion of cell proliferation, and prevention of cell apoptosis, which ultimately leads to indefinite growth of prostate cancer cells.

KIF14 regulated genes involved in cell cycle and apoptosis
Our functional study suggested that KIF14 is a potential oncogene in PCa. To explore the underlying mechanism by which KIF14 exerts oncogenic function in PCa, we performed RNA-Sequencing analysis to identify genes that were differentially expressed in KIF14 knockdown DU145 cells and control DU145 cells. Genes with 2 fold changes were considered as significant ( Figure  4A). Firstly, we analyzed the candidate genes by KEGG (Kyoto Encyclopedia of Genes and Genomes) Pathway enrichment analysis. We found that knockdown of KIF14 led to transcriptional changes of genes involved in p53, TGF-beta, Hippo and PI3K-AKT signaling pathways ( Figure 4B). The involvement of these cancer related pathway indicated that KIF14 has functional role in carcinogenesis of PCa. We further confirmed the expression of potential KIF14 target gene by quantitative real-time PCR. In both DU145 and PC3 cells, knockdown of KIF14 upregulated the gene expression of GADD45A, GADD45B, p21, PIDD and Shisa5, which contribute to growth arrest and apoptosis induction, while downregulated CCNB1 that promotes cell cycle progression ( Figure 4C and Figure  4D). Our data suggested that KIF14 may function as a candidate oncogene through regulating genes involved in cell cycle modulation and apoptosis, thus leading to tumor progression in PCa.

KIF14 was overexpressed in primary PCa and positively correlated with clinical stage in patients with PCa
Our study suggested that KIF14 was upregulated in PCa cell lines. We further investigated whether KIF14 was overexpressed in primary human PCa specimens by immunohistochemical staining. KIF14 was highly expressed in 57% (37/65) primary PCa tissues, while it was either absent or weakly expressed in another 20 BPH patients (20/25) (Table 1). Thus, KIF14 was significantly overexpressed in primary human PCa (P = 0.0017, Table 1). To further clarify the relationship between KIF14 expression and PCa progression, we investigated its expression and clinicopathological variables in 65 patients ( Table 2). KIF14 expression was positively correlated with histological grade ( Figure 5A). Moreover, KIF14 expression was significantly correlated with clinical stage (P = 0.0002) and Gleason score (P = 0.0041). However, no correlation was found between KIF14 expression and age and metastasis.

Prognostic significance of KIF4 expression
Clinical Gleason Score 8-10 is a negative independent prognostic factor on the progression free survival. In the present study, whereas 36% (5/14) of patients with Gleason Score≤6 showed low KIF14 expression, 75% (21/28) of   In 20/25 BPH specimens, KIF14 was either absent or weakly expressed. However, the PCa tissue samples (37/65) were scored strong positive or moderate positive expression. The difference was significant. (P < 0.05).
men with Gleason Score of 8-10 displayed high KIF14 expression (P = 0.0041, Table 2), indicating a positive correlation of KIF14 overexpression with poor prognosis. Moreover, Kaplan-Meier analysis demonstrated that patients with higher KIF14 expression had a significantly shorter 5-year overall survival time ( Figure 5B). Thus, these data emphasize the close association between KIF14 expression and prognosis in patients with PCa.

DISCUSSION
Tumorigenesis is a complex process which contains lots of abnormal disruptions of tumor suppressor genes, oncogenes and their related signaling pathways [17][18][19][20][21][22][23][24][25][26]. PCa is one of the most common non-cutaneous malignancies among men worldwide with a peak incidence in men of ~ 70 years [27]. Though some of the prostate cancer patients are no need to exposure to treatment, about one-third of patients with organ-confined PCa fail surgical treatment and progress to aggressive or metastatic lethal disease in 10 years [2]. Therefore, understanding mechanisms underlying the tumorigenesis and development of PCa is needed for early diagnosis, risk assessment and selection of therapeutic strategies against PCa. In this study, we identified KIF14 as a potential candidate oncogene involved in tumor progression and poor prognosis in prostate cancer.
KIF14 is a member of kinesin family proteins involved in cell division, microtubule polymer dynamics, intracellular transportation and signal transduction [3]. There is increasing evidence that KIF14 is overexpressed in multiple types of tumors, including hepatocellular carcinoma, lung cancer, breast cancer, etc [6][7][8][9][10][11][12], implicating its role as an potential oncogene. The upregulation of KIF14 might be caused by genomic amplification [6,9,12]. It was reported that in hepatocellular carcinoma, loss of KIF14 downregulates the expression of Skp2 and Cks1, leading to accumulation of p27Kip1. Downregulation of Skp2 and Cks1 also resulted in cytokinesis failure, which may inhibit tumor growth [6]. Inhibition of KIF14 inhibited the migration and induces apoptosis via inactivation of Akt kinase in HCC cells [7]. The effect of KIF14 on AKT phosphorylation was also reported in triple-negative breast cancer [9]. However, the expression of KIF14 and its role in prostate carcinoma remain unclear thus far. The present study demonstrated for the first time, to our knowledge that KIF14 is a candidate oncogene in PCa, which is upregulated in multiple cell lines of PCa and primary human PCa tissues. Remarkly, KIF14 depletion induced cell cycle arrest, reduced proliferation and apoptosis in PCa cells at least partially through upregulating GADD45A, GADD45B, p21, PIDD and Shisa5, which contribute to growth arrest GADD45 mediated G2/M checkpoint in human and murine cells [28], this effect may be due to GADD45 ability to dissociate complexes of Cyclin B1 and Cdc2 [29]. The dephosphorylation of the CDC2/ CCNB complex mediated by Cdc25 allows cells entry into M phase [30,31]. The main function of p21 in cell cycle regulation is inhibiting the activity of cyclin B/ Cdk1 complexes. Thus p21 is responsible for the nuclear retention of the cyclin B/Cdk1 complex and cyclin B degradation in response to DNA damage [32]. In our study, we found CCNB was repressed, while GADD45A, GADD45B and p21 were upregulated after KIF14 inhibition. The changed expression profile of these genes could explain G2 arrest in KIF14 knockdown cells. G2 arrest will in turn lead to lower proliferation rate in cells and it could explain the inhibition of proliferation in KIF14 knockdown cells. We also found that SKP2 was downregulated in our microarray study (data not shown), which is consistent with other group's study [6]. We also found that PIDD, which interact with other death domain proteins and function as an adaptor protein in cell death-related signaling processes [33], was upregulated in KIF14 downregulated cells. PIDD is an effector of p53-induced apoptosis and it could act as switch between cell survival and cell death in response to DNA damage [33]. Shisa5 is another apoptotic protein upregulated in KIF14 knockdown cells. Shisa5 is localized to the endoplasmic reticulum and together with p53 induces apoptosis in a caspase-dependent manner [34,35]. The upregulation of these two proteins will induce apoptosis in KIF14 knockdown cells. Our future study will be focused on how KIF14 regulates those genes' expression.
Except for KIF14, many other kinesin proteins are reported as oncoproteins [36][37][38][39][40][41][42][43]. For example, KIFC1 induces resistance to docetaxel and is associated with survival of patients with prostate cancer [36]; KIF11 is a driver of invasion, proliferation, and self-renewal in glioblastoma [37]; KIF3C promotes tumor growth and metastasis in breast cancer by regulating TGF-β signaling [38]; It is also reported that inhibition of KIF22 suppresses cancer cell proliferation by delaying mitotic exit through upregulating CDC25C expression [39]; KIF3a can promote proliferation and invasion via Wnt signaling in advanced prostate cancer [40]. All these further support our finding that KIF14 is a functional oncoprotein in prostate cancer.
It has been reported that KIF14 is a marker of poor prognosis in patients with hepatocellular carcinoma and ovarian cancer, where KIF14 overexpression enhances tumor growth, while its knockdown decreases tumorigenicity in vitro and in xenografts [13]. However, the role of KIF14 as an indicator in PCa is largely unknown. Here, we found that KIF14 was upregulated in most clinical PCa tissues and cell lines. Subsequently, by immunohistochemical methods, we observed a positive correlation between KIF4 expression and histological grade, clinical stage and Gleason Score. Kaplan-Meier analysis demonstrated that KIF14 was a prognostic factor for patient survival. Our results thus indicate that KIF14 could be served as a novel and useful prognostic biomarker for PCa.
In summary, our study found that KIF14 is upregulated in PCa cell lines and clinical PCa tissues. Inhibition of KIF14 suppresses cell proliferation, induces G2 arrest and apoptosis. This effect was mainly caused by regulating GADD45A, GADD45B, CCNB, p21, PIDD and SHISA5 expression. Moreover, the protein level of KIF14 is positively correlated with histological grade, clinical tumor stage and poor prognosis in patients with PCa. Thus, our results suggest that KIF14 could act as a candidate oncogene that contributes to tumor progression and poor prognosis in PCa. KIF14 may represent a novel molecular target for the treatment of PCa.

RNA extraction, reverse transcription and quantitative realtime PCR
RNA was extracted with Trizol reagent according to manufacturer's protocol (Invitrogen, Eugene, OR, USA). One normal prostate RNA was purchased from Clontech (Clontech, Palo Alto, CA). Reverse transcription was performed with Random hexamers and SuperScript-III (Invitrogen, Eugene, OR, USA). Quantitative realtime PCR was carried out with the Applied Biosystems 7300 real-time systems using real-time PCR Master Mix (SYBR Green). The standard curve was used to establish amplification efficiency [45]. All primers used are listed in Table 3. Each experiment was performed in triplicate in three independent experiments.

Growth curve
DU145 or PC3 cells were seeded in six-well plates at an appropriate density with penicillin/streptomycin free medium. Cells were transfected with control siRNA or siKIF14, respectively. After transfection, cell numbers were counted every 24 hours including the start point. Each experiment was conducted in triplicate in three independent experiments.

Flow cytometry analysis of cell cycle
About 48 hours after transfection, cells were harvested and fixed with ice-cold 70% ethanol for 24 hours. After that, cells were washed with PBS once and stained with 100 µl of 50 mg/L propidium iodide (PI) for 15 min at 4°C in dark. The cell cycle profiles were assayed by flow cytometry (excitation at 536 nm and emission at 617 nm) with a flow cytometer (FACS Calibur TM , BD Biosciences, CA, USA) and data were analyzed with the ModFit LT software. Each experiment was conducted in triplicate in three independent experiments.

Apoptosis assay
Apoptosis was assessed by flow cytometry analysis. Apoptotic cells were detected with PE Annexin V Apoptosis Detection Kit I (BD Bioscience). About 48 hours after transfection, cells were harvested and centrifuged at 1000 rpm for 10 min, rinsed with cold PBS once. About 1 × 10 5 cells were stained with 5 µl PE Annexin V and 5 µl 7-AAD working solution for 15 min at room temperature in dark. The ratio of apoptotic cells were analyzed by flow cytometry. The total number of early apoptotic cells and late apoptotic cells were considered as total apoptotic cells. Each experiment was conducted in triplicate in three independent experiments.

RNA sequencing
RNA was extracted with an RNeasy Mini kit (QIAGEN). Firstly, mRNA was purified and RNA-Sequencing library preparation was performed according to the manual of manufacturers (KAPA biosystems). Sequencing reactions were performed with the Illumina HiSeq platform. RNA-seq reads were mapped to the human genome (hg19) using Burrows-Wheeler Aligner (bwa). The duplicate reads were then marked by picard. The HTseq tool was used to calculate the reads count for each gene. Finally, we used the Reads Per Kilobase per Million mapped reads (rpkm) command in edgeR package to calculate the rpkm of each gene. Specifically, we added one to the reads count for each gene and take log values. These log values were inputted into the heatmap.2 in R package gplots to generate the heatmap figure [46]. The genes with 2 fold changes were loaded to David bioinformatics database (https:// david.ncifcrf.gov/) for pathway analysis.

Oncomine analysis
The Oncomine cancer microarray database (http:// www.oncomine.com) [14] was used to analyze expression profiles of KIF14 in a variety of human prostate cancer and normal tissues. Two representative prostate cancer datasets [15,16] were presented in this study. A foldchange of 2 for KIF14 gene expression compared to control was required for criterion.

Patient samples and IHC
All patient samples were collected in Shenzhen People's Hospital with the permission of patients. All PCa and benign prostatic hyperplasia (BPH) specimens were obtained immediately after operation (radical prostatectomy, transrectal biopsy or transurethral resection of the prostate). Tissues were then fixed in 10 % buffered formalin and embedded in paraffin. Tissue sections were stained with KLF14 antibody purchased from abcam. The protocol of IHC staining was described before [47]. The clinical characteristics of the PCa cohort are listed in Table 2. The tumor staging and differentiation of PCa was defined according to the 2010 AJCC Cancer Staging Manual seventh edition. (Well differentiation: Gleason Score ≤ 6) Moderate differentiation: Gleason Score = 7; Poor differentiation: Gleason Score ≥ 8)

Follow-up
All patients with PCa received check-ups every 2-3 months during the first 2 years and every 3-6 months there after until the follow-up period ended in January 2017. Surgeons who were blinded to the study performed the follow-up study. The overall survival was defined as the length of time between surgery and either the death of the patient or the last follow-up visit.

Statistical analysis
Data are presented as mean ± standard deviation. Statistical assessments were carried out using Student's t test. P < 0.05 was considered statistically significant.


extfieldsofstudy:
Biology
Computer Science
Medicine

text:
Semi-automated quantitative Drosophila wings measurements

Background Drosophila melanogaster is an important organism used in many fields of biological research such as genetics and developmental biology. Drosophila wings have been widely used to study the genetics of development, morphometrics and evolution. Therefore there is much interest in quantifying wing structures of Drosophila. Advancement in technology has increased the ease in which images of Drosophila can be acquired. However such studies have been limited by the slow and tedious process of acquiring phenotypic data. Results We have developed a system that automatically detects and measures key points and vein segments on a Drosophila wing. Key points are detected by performing image transformations and template matching on Drosophila wing images while vein segments are detected using an Active Contour algorithm. The accuracy of our key point detection was compared against key point annotations of users. We also performed key point detection using different training data sets of Drosophila wing images. We compared our software with an existing automated image analysis system for Drosophila wings and showed that our system performs better than the state of the art. Vein segments were manually measured and compared against the measurements obtained from our system. Conclusion Our system was able to detect specific key points and vein segments from Drosophila wing images with high accuracy. Electronic supplementary material The online version of this article (doi:10.1186/s12859-017-1720-y) contains supplementary material, which is available to authorized users.

Background
Recent advances in high-throughput sequencing technology have enabled us to obtain genome information from any kind of organism [1]. For humans, more than a thousand of genome sequences and single nucleotide polymorphisms among them have been characterized [2]. This further motivated genome-wide association studies (GWAS) in the literature with the hopes of discovering the variations responsible for genetic diseases and traits in the genome sequences. The strategy of investigating the relationship between genotype and phenotype should succeed in discovering the genetic basis for any inherited traits in any kinds of organisms. However, in many higher animals *Correspondence: leehk@bii.a-star.edu.sg † Equal contributors 1 Imaging Informatics Division, Bioinformatics Institute, 30 Biopolis Street, #07-01, Matrix, Singapore 138671 Singapore, Singapore 3 Research Center for Genomics and Bioinformatics, Tokyo Metropolitan University, Hachioji, Tokyo 192-0397, Japan Full list of author information is available at the end of the article other than human beings, the description of phenotype needs manpower with expertise in morphology of the animals under consideration; such a dependency of this approach on available manpower may be a critical bottleneck of GWAS. Hence, development of an automated system for describing phenotype may be one of the best solutions to this problem. Fruit flies, in particular Drosophila melanogaster is an important organism used for biological research fields particularly in genetics and developmental biology. Moreover, the complete genome sequences of this species and its related species [3] opened up an ideal opportunity for comparative genomics analyses for the systematic understanding of phenotype and genotype relationships. Species belonging to the genus Drosophila are particularly useful for this purpose because there are more than 2000 described species in the genus Drosophila [4], each of which has distinct phenotypic characters. On the other hand, this situation makes it difficult for ordinary Drosophila researchers to distinguish morphological differences among different species and only a limited number of researchers who are experts in Drosophila taxonomy can recognize key morphological characters to identify the species.
The use of wing morphometrics for species identification have been demonstrated on a variety of insect species [5][6][7][8][9], because of the simple two dimensional structures with clear-cut patterns of veins and their crossings, which are frequently utilized as landmarks to measure distances and angles among them. There is much interest in quantifying wing structures of Drosophila [10][11][12][13][14] and Drosophila wings have been widely used to study the genetics of development, morphometrics and evolution [15][16][17][18][19][20]. Advances in technology have increased the ease in which images of insect wings can be acquired. However the extraction of phenotypic information for further research and analysis is a manual and time consuming process. This has motivated the creation of several software which aims to automate this task. One common method involves using software for easy annotation of landmarks on the wing (http://life.bio.sunysb.edu/ ee/rohlf/software.html, http://www.hockerley.plus.com/). However annotating landmarks is time consuming and prone to errors [21]. More advanced software attempt to find the location of landmarks and vein lengths using a variety of image processing and optimization methods. These methods tend to involve fitting a set of bezier splines onto the wing veins. WINGMACHINE [22] is a program that automatically detects and measures the positions of veins and edges of wing blades from live flies. MorphoJ [23] performs geometric morphometrics analysis on morphological landmark inputs. Although MorphoJ does not detect and extract landmarks, it is a useful tool for studying data on shape combined with molecular genetics or ecological information. Automated systems to detect and extract information from wings from a variety of other insect species exist [24][25][26][27][28]. However only a few of these programs are able to extract specific unique key points that are used for the prediction of the fly species. Crnojevic et al. [28] trained a svm classifier based on Histogram of Oriented Gradient (HOG) and Complete Local Binary Pattern (CLBP) features for vein junction detection of hoverflies. In addition the junctions were used to construct convex hulls which were used to discriminate 4 different hoverfly species. Thus a more specific and customized software is needed to automate the detection and prediction process of Drosophila wings. As mentioned previously, the wing morphology of Drosophila species is highly favorable for automated image analysis [22]. Therefore, in this paper we introduce an automated system that locates specific keypoints on a Drosophila wing image and performs morphometrics measurement of several important vein segments. The extraction of these phenotypic data serves as a base for future studies and research of Drosophila wings such as specie prediction.
The software requires the user to annotate three specific key points on the fly wing before key points and vein segments are calculated automatically. The program provides an intuitive user interface and users can manually edit key points and vein segment lengths if the prediction is inaccurate. In this contribution, we report the performance of our program and perform a comparison against WING-MACHINE [22] an existing program that automates the measurement of Drosophila wings.

Data set
A total of 959 flies from 16 (Table 1). A one-sided wing was taken from each individual and placed in 8 μl mounting solution (ethanol : glycerol = 2 : 3) dropped on a glass slide. After a coverslip was placed, the slide was placed under a stereomicroscope (Nikon SMZ-10A) and the wing photograph was taken by a digital camera (Canon EOS 40D), where the original image resolutions was 1936×1288, 2816×1880, or 3888×2592 pixels. A total of 600 images were used for image processing. Figure 1 shows the arc-lengths to be measured as well as the mean and standard deviation of the pixel intensities of all 600 images. To conserve memory usage and improve processing speed, small image patches around the manually annotated key points were cropped and stored in small image patch files. In practise, only image patches near the key point locations were needed for key point detection. All images are also converted to gray scale images. Let annotated images be denoted by D m , where m = 1 · · · N and N = 600 is the total number of images.

Key point detection
To identify a species from the wing image 13 key points are needed. The key points are labelled 'a' to 'm' and can be found in Fig. 1. Key point 'k' is the boundary of the costal fringe and the others are intersections of veins. Given a new, unannotated fly wing image I, key points detection is carried out in two steps. The first step aligns the annotated images D m and the input image I. In the next step, key points are detected using template matching.
To align the input image I with D m , the user is required to click the location of the three key points 'a' , 'h' and 'k' on I (see Fig. 1). A screen shot of the software can be found in Fig. 2. An affine transformation matrix Q m is then calculated by matching the key points 'a' , 'h' and 'k' in I and D m . Hence, each image in the data set D m is associated with one affine transformation Q m , m = 1, · · · N. This transformation matrix is then used to transform all images in the data set from D m to D m . Figure 3 shows an example of transforming D m to D m . In our implementation, only pixels within a small image patch near the annotated key points of D m are transformed. This process greatly increase the processing speed. As an affine transformation is used to align the images, the software is capable of handling fly wings of any orientation.

Template matching
Image patches between I and D m are matched and the best match locations are predicted to be the key point locations. Define w I (x, y, θ) to be an image patch centred at (x, y), oriented at an angle θ extracted from the image I. Similarly, let w D m (x, y, θ) be an image patch extracted from the image D m . The objective is to find the best match between w I and w D m by adjusting the image patch parameters m, x, y, θ.
We shall describe the template matching procedure using key point j as illustration. The procedure for matching other key points is similar. The location of key point j in D m , x m j , y m j is known since D m has been manually annotated. The transformed location x m j , y m j can be calculated using Q m . Five templates w D m x m j , y m j , θ are constructed with θ = −10, −5, 0, 5, 10 degrees. Fig. 4c,  examples of these five templates. Next, we want to match these templates to a target location w I x m j , y m j , 0 in image I. One last step before the matching process is to normalize the pixel intensities within the patches by performing a linear transformation to make the range of pixel intensities within the image patch between 0 and 255.
Matching is done by sliding the windoww D m x m j , y m j , θ with respect tow I x m j , y m j , 0 .w D m x m j , y m j , θ and w I x m j , y m j , 0 are the normalized image patches. Figure 5 illustrates the matching process. The matching score is given by, c x , c y is the shift of the center locations between the image patches, with −15 ≤ c x , c y ≤ 15. · is the Euclidean norm. Finally, the predicted key point location corresponds to the best match among all shifts, orientations and template images D m .
and the key point j's coordinate is, The above process is repeated for the remaining key points on I. The details of the above algorithm can be found in Algorithm 1 KEYPOINT-DETECTION in Additional file 1.

Arc length calculations
Tracing the arcs as depicted in Fig. 1 is done using active contours while the Euclidean distances between key points are calculated trivially. In the active contour method, an estimate of the desired curve is first generated using a template. The curve is then evolved using a variational principle via a gradient descend method. Two essential ingredients for applying the active contour method are estimation of the initial curve and the existent a b c

Initialization of active contours
The initialization of active contour curves is critical for finding the fly wing vein accurately. If the initial curve lies far away from the vein, active contour might evolve this arc to another vein.
For example, setting the initial curves to straight lines works well for arcs fg, gh, jl and kl as these arcs are relatively straight. This is not the case for arc ml as there is another vein just below it. Active contour will tend to evolve this curve towards the vein below ml.
A good initialization of curves can be obtained by mapping veins of a template fly wing image to the input image I. This template wing should provide a good representation of wings for all different fly wing species. This mapping is done by aligning the template fly wing image with the input image. The alignment process has been described in the key point detection section.  Figure 6 illustrates this concept. First a fly wing image is chosen to be the template wing image (Fig. 6a). We have selected our template wing such that key point k is equidistant to both key point j and l. Mapped arcs will be warped badly if key point k is not equidistant to both key points j and l. Different arcs are drawn using different colours in Fig. 6 for clarity. Arc ml is drawn outside of the template wing due to two reasons. The first is that the arc ml of the template wing that we have chosen is quite flat and does not provide a suitable representation of other fly wings' ml arcs, which have more curvature. Secondly, it can be noticed that there is another vein just below arc ml that lies close to arc ml. Active Contour might evolve the mapped arc to the wrong vein if the mapped arc happens to lie closer to that vein and not arc ml. To prevent this, the arc ml is drawn outside the template wing image so that the mapped arc also lies outside the fly wing. Some results of this mapping are shown in Fig. 6b -d. One can now notice how the mapped arc ml of Fig. 6c lies outside the wing, rather than in between the two veins. Figure 7 shows how the input image is processed to produce a gradient for active contour optimization. The main idea is to make the veins' pixels very bright and for them to be surrounded by pixels that become darker as they become further away from the veins. A step-by-step work flow of how the preprocessed image is obtained is shown in Fig. 7a. The original image (Fig. 7b) is blurred and inverted to obtain an image where the veins' pixels are very bright. The result is shown in Fig. 7c. A second image is needed to help produce a smooth downwards gradient away from the vein. Edges (Fig. 7d) are found on the blurred image before applying a threshold to obtain a binary image (Fig. 7e). The binary image is dilated to make the fly wing veins thicker before inverting the image (Fig. 7f). A distance transform (Fig. 7g) is applied before inverting the image one last time (Fig. 7h). The two images ( Fig. 7c and h) are then added pixel-wise to produce the preprocessed image P (Fig. 7i). This process is integrated in the software and is done automatically before performing Active Contour.

Active contour formalism
Given an image P, the arc length between 2 key points can be found by finding the path that follows high intensity pixels between the two points. Curves are also represented by connecting short straight line segments. Given 2 key points p a = (x a , y a ) and An objective function for the active contour can then be defined as, that f is small along paths of high image intensity. To make the active contours invariant with respect to fluctuations of the overall image intensity, we make f to be invariant to image intensity multiplication, i.e. f (x 0 , . . .y n , P) = f (x 0 , . . .y n , aP). Hence f can be defined as: P is the normalization on the image intensity such that f becomes invariant under image intensity multiplication. 0 < 1 is a small regularizer to prevent numerical overflow. The desired curve can be found by minimizing the objective function: using gradient descend method, n−1 (8) Finally, the arc length can then be calculated by: The details of the algorithm can be found in Algorithm 2 ARC-LENGTH in Additional file 1.

Species identification from fly wing images
The Drosophila species can be identified from the wing image using 13 key points (Fig. 1). The centroid of wing was estimated as the centroid of the octagon whose vertices were the key points 'a' , 'h' , 'i' , 'j' , 'k' , 'l' , 'm' and 'b' . The centroid of the octagon was computed from the areaweighted average of the centroids of the five triangles split from the octagon. Then, Euclidean distances from the wing centroid to the 13 key points were measured and normalized by dividing by the distance between 'a' and 'j' (wing length).
In addition to the normalized distances from 13 key points to the centroid, we measured three traditional wing indices, costal-index, c3 fringe and 5x-index [29][30][31] to identify the species from a wing image. The costal-index is defined as the ratio of the length of the vein l-m to that of the vein j-l. The C3 fringe is the ratio of the length of the vein k-l to that of the vein j-l. The 5X-index is the ratio of the length of the vein g-h to that of the vein f-g. Because arc lengths of wing veins were required to compute these indices, we used a fitting of veins to Bezier paths. As a result, for each wing, we obtained 16 variables to construct the training dataset and measured their values for 589 flies from 15 species excluding Microdrosophila sp., which was used as a sample for a negative control as will be explained later.
With this size of the training dataset, however, utilization of all of the 16 variables did not show the best performance in the species identification in our preliminary test, where the data from 588 files were used to construct the training dataset and the remaining data from one fly was used as the test data. Therefore, we examined the conditions that may give the best result among all the possible numbers (2-16) and possible combinations of these numbers of variables (65519 patterns in total) and found that the best result was obtained when six variables, the distances from the key points e, f, g and m to the centroid, costal-index and c3-fringe, were employed. Accordingly, we used only these six variables for the following discriminant analysis.
The mean vector and the variance and covariance matrix of the six variables were computed for the 15 species and for male and female separately excluding Drosophila hydei female and Scaptodrosophila coracina female due to a small sample size. Therefore, we obtained the vectors for 28 groups in total. Then, the squared Mahalanobis' distance, D 2 , [32] from each wing to each group was computed. Following De Maesschalck et al.
denotes the vector of j-th wing,μ i andŜ i denote the mean vector and the variance and covariance matrix for i-th group, with the aid of the R package ( [34]). For each wing image from unknown sample, the D 2 values were computed withμ i andŜ i for the 28 groups and the group that gave the smallest D 2 value was identified to be the species of the sample. The statistical significance of the goodness of fit to the inferred group was obtained by the chi-square test with n degrees of freedom. When the smallest D 2 value was larger than the expected value for the best-fit species at the 0.1% level, the sample's species was identified as 'unknown' .
To examine the efficiency of species identification by the discriminant analysis, we used 370 additional wing images not included in the training data from 15 species including Microdrosophila sp., which were absent in the training data and thus expected to be identified as 'unknown' species if the method worked correctly.

Evaluation metric of key point detection
The key points of each fly wing in the data set is found using the methods described above. To test the accuracy of the key point detection algorithm, 15 fly wing images, one from each species, is manually annotated by 10 users. This give us 10 X,Y coordinates for each key point for each of the 15 fly wing images. The covariance matrix, which represents the spread of the 10 annotated points for a single key point, can be obtained from the manually annotated data. These covariances matrices are then used for benchmarking the accuracy of the automated methods.
A total of 15 covariance matrices, 1 for each annotated image, are obtained for a single key point. The 15 covariance matrices are summed and averaged to obtain the final covariance matrix for a single key point. This process is done for all key points. Table 2 shows the covariance matrix obtained for each key point. One might notice that the covariance matrix of key points 'a' , 'b' , 'c' , 'k' and 'm' contains larger values. This indicates that the manual annotations of these key points have a large spread and it is not easy even for humans to locate these key points consistently.

Key point detection results
Some results of Key Point Detection are shown in Fig. 8. To determine if a predicted key point is accurate, we check if the euclidean distance (in pixels) between the predicted key point and the ground truth is less than some threshold. The threshold used in our experiment is 2 standard deviations away from the corresponding key point's covariance matrix. The results of this experiment are summarized in Table 3. The mean deviation in pixels for each point is also reported. Accuracy for key points 'a' , 'h' and 'k' are not reported as these points were marked by the user. The mean pixel deviation is also slightly higher for key points 'b' , 'c' , 'l' and m' . These key points are located on slightly thicker fly wing veins as compared to other key points. A thicker vein would increase the candidate area for a key point and this explains the larger values in pixel deviations.
Different sized data sets were used to predict key points. For the data set of size 100, 100 images had key points manually annotated and used for template matching. The remaining 500 images were used to test the accuracy of the algorithm. Cross validations was also done for this data set. Similarly, for the data set of size 599, 599 images had key points manually annotated and the remaining image was used for testing the algorithm. We expect that the data set of size 599 will yield better results. The accuracy results are shown in Table 3. We also downsized our images further and repeated the above experiments. The results of the downsized dataset can be found in Table 4.
We also compared our method against an existing fly wing fitting software, WINGMACHINE [22]. WINGMA-CHINE first detects curves in a fly wing image that fit a priori template before fitting spline curves to the image. The location of two landmarks -the humeral break and alula notch, need to be known before curve fitting can take place. We ran the Wings software on ten images in our dataset. The results are shown in Fig. 8g -i and Table 5. WINGMACHINE tries to find the fly wing veins using edge detection methods and results were not so ideal in cases where there were background artefacts. Moreover, the edge detection algorithm sometimes fail to pick up some veins that are too faint. This can lead to a wrong prediction on the joint locations which results in a bad fit for the spline curves.

Arc length results
Some examples of good results of Active Contour are shown in Fig. 8. Note that even though there are dark patches (Fig. 8d), background artefacts (Fig. 8e) or spots on the wing (Fig. 8f), Active Contour is still able to find the fly wing vein accurately most of the time.
We found a total of 4 bad results out of 600 images in our dataset. Users are allowed to edit bad results of arc lengths in the fly wing software. Some examples of bad results of Active Contour are also displayed in Fig. 8. As can been seen in Fig. 8j and k, air bubbles are present in the background while dust is present in Fig. 8l. It can be observed that only the location of arc lm is incorrect and that artefacts in these images lie above arc lm. It has been mentioned previously that the arc lm is initialized above the fly wing. These artefacts in the background produces large peaks in the active contour gradient image, resulting in Active Contour getting stuck in the local maxima of the artefact. Other arcs do not suffer from this problem as these arcs are initialized onto the foreground which has little artefacts. Thus only artefacts that lie just above the fly wing affects the accuracy of Active Contour.
The accuracy of our Active Contour method was also tested. Arc lengths ground truth of fifteen different fly wings, one for each species, are annotated manually and compared with the arc lengths found using Active Contour. The error is calculated as the area bounded by the annotated arc and the predicted arc. The error is then divided by the length of the annotated arc. The results are tabulated in Table 6. This is normal as a longer arc spans a larger area of the image and this increases the amount of noise and background artefact that might affect the Active Contour algorithm.

Species identification results
The efficiency of the wing image analysis for species identification is shown in Table 7. Out of 370 flies examined, 346 flies were identified as the correct species (success rate was 94%), whereas the capability of identifying correct sex was poor (success rate was 56%). This suggests that the wing morphology is not significantly different between female and male, whereas the between-species difference is significantly larger than the within-species difference among all of the 15 species examined. In addition, it is noteworthy that the present system perfectly recognized Microdrosophila sp., which was not included in the training data, as 'unknown' species. This suggests that the present system has a discrimination power to give an appropriate answer even when the species of sample is not included in the training data. Therefore, the present system is expected to be useful to analyze wing morphology for species identification among species belonging to Drosophila and related genera, even after more data for more species are added. shows arc lengths found using our method. The colours red, blue, orange, green and magenta represents the arc segments fg, gh, jl, kl, lm respectively. g -i shows key point and arc length predicted for the same images using the WINGMACHINE software. Red arrows are used to indicate the errors for clarity. The errors in g and i are due to incorrect prediction of joint locations while the error in h is due to a bad fitting results of an a priori template wing. Lastly, j -l shows three examples of bad arc length estimation results found using Active Contour due to artefacts in the background such as air bubbles and dust particles. The red arrows show the artefacts that affect the Active contour algorithm. There are a total of four images with bad arc length estimation out of our dataset of 600 images

Discussion
Our software has managed to successfully predict landmarks on up to 15 different species of fly wing images. Incorrect landmarks detected can be corrected using the software user interface with ease. In addition, our software is able to accurately calculate and measure specific vein segments which are required for fly specie detection. Incorrect vein segments found can also be corrected using the user interface of the software. In contrast to WINGMACHINE software, which makes use of high level features to optimize fit of an a priori model of a fly wing, our software utilizes low level features such as edges and pixel values to find land marks and vein segments. As mentioned by WINGMACHINE, developing detection algorithms using low levels features tend to encounter problems such as varying thickness of veins, uneven lighting and background artefacts which cause inaccuracies during prediction. However, our software manages to overcome the majority of these difficulties and performs well for images in our dataset.
One advantage of using WINGMACHINE over our software is that only two land marks are needed for detection instead of three. However, it should be noted that WINGMACHINE software may require pretreatment of image data and this may require more manpower as compared to our software. The advantage of using MorphoJ is that phenotypic data can be easily read and analysed within the software. However, MorphoJ does not have the capabilities to detect and extract these information from images automatically.
There are a few drawbacks for our algorithm. The first is the need for a large dataset. As demonstrated in the key point detection results section, a large dataset containing fly wing images of numerous fly species would help increase the accuracy of our algorithm. This results in a slow detection process as each image in the dataset has to be loaded. To reduce the time taken for the detection process, image patches centered on the fly wing's key points are cropped from the original fly wing image and saved in a separate image. This separate image is loaded

Conclusion
In this paper, we have demonstrated how to locate Drosophila wing key points and obtain arc lengths of the wing veins using image processing methods. Our results show that our method is accurate and we compared our The results of using different a dataset of size 599 and 100 are shown as well. As shown in the table, using a larger dataset yields a more accurate result. Also note that key points a, h, k are excluded as these key points are annotated by the users and are thus not being predicted   Table 7 Efficiency of species identification using the wing image analysis


extfieldsofstudy: Biology Medicine

text:
The impact of host sex on the outcome of co-infection

Males and females vary in many characteristics that typically underlie how well a host is able to fight infection, such as body-size, immune capacity, or energy availability. Although well studied in the context of sexual signalling, there is now growing recognition that these differences can influence aspects of pathogen evolution as well. Here we consider how co-infection between multiple pathogen strains is shaped by male-female differences. In natural populations, infections by more than one pathogen strain or species are believed to be a widespread occurrence. Using the water flea, Daphnia magna, we exposed genetically identical males and females to replicated bacterial co-infections. We found that pathogen transmission and virulence were much higher in females. However, males did not simply lower average pathogen fitness, but rather the influence of co-infection was more varied and less defined than in females. We discuss how pathogens may have more fitness benefits to gain, and consequently to lose, when infecting one sex over the other.

Males and females of the same species can be strikingly different from each other. They may differ in many overt traits, such as size, morphology or colouration 1 ; but these differences can also extend to the immune response, with one sex often being more severely affected by infections 2,3 . In birds and mammals, for example, males are often more susceptible to infection and develop higher parasite loads than females 4,5 , whereas the reverse pattern is often observed in invertebrates [6][7][8] . When combined with the diversity of sex ratios, sex biases in dispersal, and general behavioural difference between males and females 1,9,10 , there is ample opportunity for dimorphism between the sexes to influence many aspects of disease epidemiology and evolution 3,8,11 .
Here we consider how differences between males and females can influence a widespread feature of host-pathogen interactions -the influence of co-infections on the outcome of infection. Multiple infections are believed to be more common in nature than infections composed of only a single parasite species or genotype 12,13 . Both theory and empirical studies have shown that the overall virulence of a co-infection can increase, decrease, or remain intermediate to the intrinsic properties of each infecting pathogen 14 . Rarely, however, have these concepts been explored within the context of sexual dimorphism in host resistance and exploitative potential. Co-infection studies typically either explore infection outcomes in only one sex 15,16 , or have not explicitly accounted for males-female differences 17 .
We explored the interplay between host sex and co-infection using the well-studied crustacean Daphnia magna and its obligate bacterial pathogen Pasteuria ramosa 15,18,19 . In Daphnia up to 50 percent of all individuals in a population can be infected with two or more pathogen species 20 , or eight or more genotypes of the one pathogen species 21 . Male Daphnia, however, are significantly smaller, more resistant, and die younger than females 22 ; potentially offering a parasite a more difficult pool of resources to exploit than females. To study how the contrasting environments provided by males and females impacted on host and pathogen fitness, we manipulated pathogen competition in genetically identical male and female Daphnia (genotype HU-HO-2), using three P. ramosa genotypes (C1, C14, and C19) with known infection characteristics. In line with theory 14 , the reduction in host lifespan (relative to uninfected controls) and the production of pathogen transmission spores were used as estimates of co-infection outcomes. We hypothesised that any differences between single and multiple infections seen in females would be dampened in males as a result of the sexual dimorphism in physiological characteristics.

Results
Our results reveal that the production of transmission spores under single and co-infections depended on the sex of the host (Fig. 1). While all pathogen genotypes performed significantly better in females (overall sex effect P < 0.01, Table 1), the influence of the pathogen combinations was sex-specific (overall infection term, P < 0.05, Table 1). This significant interaction was driven by the co-infection combinations where the two pathogens differ   Table 1. Influence of host sex and co-infection on the production of transmission spores, and the relative reduction in host lifespan. In bold are the overall tests using the entire dataset for the influence of host sex on the outcome of all infection treatments (all single and co-infections together). Below these are the individual analyses for each co-infection combination.
We compared virulence under coinfection with that expected based on interference competition (intermediate values), competition for resources (higher values than most virulent) or spite (lower values than the least virulent). Post-hoc comparisons, as indicted in Figs 1 and 2, revealed that co-infection scenarios similar to resource competition or spite are least likely, as co-infection never significantly exceeded the most virulent pathogen, nor was lower than the least virulent. Instead, co-infection largely matched the average of the innate characteristics of each infecting genotype in isolation; albeit with some sex-specific differences. When the two pathogens differed in their virulence or spore production, co-infection in females was indistinguishable from the average of the two single infections, and different to that predicted by either the least or most virulent strain (see Table 2). In contrast, coinfection in males was more varied, with one co-infection (C1 + C19) not significantly different to the predicted intermediate values, and another (C14 + C19) more similar to the most virulent and least transmissive strain (C19).

Discussion
Males and females vary in many characteristics that combine to offer pathogens two environments that differ in their exploitative potential 3 . Here, males of Daphnia magna are smaller and more resistant to infection than their genetically identical female counterparts 22 , a pattern common to many other arthropods 6,7 . Our results show that these differences change how well a pathogen can exploit a host, with pathogen transmission (mature spore loads) and virulence (reduction in expected lifespan) being much higher in females. However, the presence of the significant interaction term for pathogen transmission indicates that males do not simply lower average pathogen fitness, but rather the transition from the worst performing pathogen, to the mixed co-infection, to the best performing pathogen was different in each sex.
Emerging from our results is that the best performing pathogens have the most to gain from single infections in females, but also the most to lose from co-infections in the same sex. High-growth pathogen strains (C1 or C14) outperformed their lesser performing competitor (C19) in single infections by 2.1 million spores on average in females, versus only 0.55 million spores in males. As the evolutionary advantage of a pathogen depends its success relative to other competitors, not its absolute transmission, this represents a nearly four-fold improvement in relative fitness. Conversely, the same pathogens also experienced the greatest reduction in transmission success when co-infections occurred in females. Total spore loads for a co-infection were always less than that of the high-growth pathogen strain in isolation; a substantial loss of potential fitness even if the other pathogen was outcompeted completely. We also cannot rule out that the lesser performing pathogen in isolation (C19) is a superior competitor in a co-infection, a fact that would further decrease the fitness of the pathogen ranked best in isolation. Overall, this loss of fitness for high growth pathogens (and potential gain for low growth pathogens) was three-fold higher in females (1.2 million spores) than males (0.41 million spores).
The distinction between males and females was more subtle when considering the virulence of co-infections. For females, virulence in co-infections closely matched those expected by interference competition or immune-mediated competition between competing strains 14 . Overall virulence of a co-infection was firmly intermediate to that of the two competing pathogens, ruling out alternate explanations that predict either an increase (competition for a finite pool of resources, or facilitation via immune depletion) or a decrease (generation of spite) in virulence 14 . Evidence for interference competition-like mechanisms in males was less well-supported, as certain co-infection combinations tracked the overall characteristics of the least virulent pathogen. Given the lack of an overall significant interaction term for virulence, however, it appears that co-infection in males does not fundamentally change the driver of overall virulence, but rather that the patterns of co-infection in males may be less defined or extreme than in females.
As with any form of sexual dimorphism, underlying the observed differences in pathogen virulence and transmission may be a range of morphological, physiological, or immune based processes. In this system, males are more resistant (in terms of infection success) and so presumably also better equipped to limit the growth of a pathogen. For male and female Daphnia, however, an obvious point of difference is the divergence in body size 22 , whereby males are approximately 25 percent smaller than females. This could also explain the results we observe; irrespective of any immune response. With the larger body size and greater capacity for gigantism (i.e. pathogen induced growth) in females, there is naturally more resources available for a pathogen to exploit (energy or space for replication, for example). Infections in males may just be more constrained by the host, limiting both the potential fitness gains, and even loses, that a pathogen can derive from infection in this sex.
Our findings divert somewhat from previous studies in this system which have typically assumed that virulence tracks the most virulent co-infecting strain 15 , but increasingly it is becoming clear that co-infections depend as much on the specific genotype of the pathogens involved, as well as that of the host, and their interaction 17,23 . We suggest that heterogeneity between the sexes is an additional factor that is likely to influence co-infection outcomes by reducing the intensity of co-infection scenarios when males are infected. Extending this work to include genetic difference in the degree of sexual dimorphism, which we expect given the strong clonal differences observed in Daphnia, would further integrate the study of host-pathogen co-evolution with the evolution of male-female differences (see ref. 3 ).

Female co-infections
Male co-infections t-value df P-value t-value df P-value

Methods
The host and pathogen genotypes used in this study were chosen as completely compatible and displayed strong genetic interactions for the severity of infectious disease in earlier studies 19,24 . The host genotype (HU-HO-2) originates from Hungary, while the three pathogen genotypes were C1 originating from Moscow, Russia; C19, derived from North Germany; and, C14 from Finland, Tvärminne. Various combinations of these clones, in some form, have been used previously in studies of co-infection 23,25 .

Production of male and female Daphnia.
To reliably induce the production of males and females from a single Daphnia genotype, we used a modified protocol based on the work of Olmstead and Leblanc 26 . Briefly, we exposed the parental generation to a juvenoid hormone, methyl farnesoate (Echelon Biosciences, product number S-0153), after their first and second clutches, and then collected the third and fourth clutches for use in this experiment. Individuals were placed in 20-mL of the standard Daphnia media (ADaM 27 , modified after 28 ), but supplemented with methyl farnesoate at a concentration of 300 μg/L. This media was changed either every two days, or when an animal produced a clutch. Importantly, we found that the hormone had no detectable effect on the subsequent fitness of the offspring. A short pulse of hormone has been used in previous studies to study male-female differences in genetically identical Daphnia 29 . To confirm that the hormone had negligible impact on the life-history of exposed animals, we conducted a series of experimental trials (see Table 3) using female Daphnia raised from mothers that were either exposed to the hormone or raised normally in standard Daphnia media. We found that there was no significant difference between the animals collected from treated and untreated mothers in the trait means and variances of both general life-history traits (lifespan and clutch production), and the characteristics of infection (infection rates; spore loads, lifespan, and clutch production of infected animals). Using this approach removes the need to stress mothers (low-light, high-density, or starvation) to induce male production; a process which can also influence infection characteristics 19,30 . Experimental infection trials. Animals were first raised individually in 60-mL jars containing 50 mL of artificial Daphnia medium following standard protocols 18,19 , and fed daily with algae (Chlorella vulgaris, up to 5 million cells animal −1 day −1 ). At four days old, individual males and females were exposed to 20,000 spores from a single P. ramosa genotype (C1, C14, and C19), or 20,000 spores containing an equal mix of two genotypes (C1 + C14, C1 + C19, or C14 + C19, 10,000 per genotype). As per standard approaches, this process was repeated the following day when the animals were five days old (40,000 spores in total) to ensure that an internal immune or physiological response, rather than moulting, contributed to differences in infection (see ref. 31 ). In total, there were 12 treatment groups (two host sex groups × [3 single infections + 3 multiple infections]), with 36 animals initially assigned to each, plus 20 unexposed controls per sex (472 initial animals).
Survival was monitored daily and animals frozen on the day of death for the assessment of infection status and spore production using a Neubauer haemocytometer. As expected, overall infection rates were lower in males (mean ± SE: 0.58 ± 0.06; range: 0.35 to 0.70) than females (mean ± SE: 0.83 ± 0.04; range: 0.69 to 0.92). Analysis of variance revealed that infection rates indeed varied significantly with host sex (χ 2 = 33.18, df = 1, P < 0.01) and the co-infection treatments (χ 2 = 18.91, df = 5, P < 0.01), but were otherwise largely consistent (interaction term: χ 2 = 6.13, df = 5, P = 0.29). Due to difference in the average lifespan between males and females (males: 33 days ± 1.9; females: 67 days ± 2.0), we subsequently focused on two traits of common currency: the production of transmission spores at host death; and, the reduction in lifespan as compared to the average of the unexposed controls.
Statistical analysis. All statistical analyses were performed in R (ver. 3.2.4; R Development Core Team, available at www.r-project.org). Initially there were 472 animals exposed to the co-infection treatments, arising from approximately 36 animals each assigned to 12 treatment groups (two host sex groups × [3 single  Table 3. The influence of the hormone treatment on the fitness characteristics of the following offspring generation. Equality of means was assessed using a least-squares linear model in the case of lifespan and spore loads, and a generalised linear model for the number of clutches produced (Poisson distribution, log link function) and infection rates (Binomial distribution, logit link function). Where appropriate, Levene's test was used to assess homogeneity of variances between the offspring of untreated and treated mothers. All measures were estimated using a single genotype of both D. manga (HU-HO-2) and P. ramosa (C14), following the same basic infection process outlined in this study.
infections + 3 multiple infections]) plus 20 male and 20 female unexposed controls. However, this sample size was reduced by the end of the experiment as: i) not all animals became infected; ii) they died before day 15 post-infection and so could not reliably have their infection status assessed; or iii) were removed from the study due to sampling errors (see Supplementary Information for more details). Before analyses, we checked for the assumptions of linear models (normality, homogeneity of variances). In the case of the spore loads at host death and the reduction in average lifespan, we used a white-adjusted analysis of variance to correct for unequal variances as implemented using the car package of R 32 . Traits were first analysed using a full-factorial analysis of variance (Type III) with the infection treatment (various single or multiple infections), host sex (male or female), and their interaction as fixed effects. For each trait of interest, we began first with an overall test, using the entire dataset (both single and multiple infection together). If a significant effect was detected (host sex, infection treatment, or their interaction) we then explored each subset of pathogen combinations (e.g. C1 versus C14 versus C1 + C14) to see if specific single versus multiple infections combinations were contributing to the overall pattern. In this way, we avoid issues of multiple testing that would arise if we explored every related combination of single and multiple infections in isolation. Finally, one sample t-tests were used to compare the values for co-infection treatment groups against an expected value based on interference competition (intermediate values), competition for resources (higher values than most virulent) or spite (lower values than the least virulent).
Availability of materials and data. Datasets supporting the conclusions of this article are available via a figshare repository at https://dx.doi.org/10.4225/03/585732493b9ad.


Identifying differential isoform abundance with RATs: a universal tool and a warning

Motivation The biological importance of changes in gene and transcript expression is well recognised and is reflected by the wide variety of tools available to characterise these changes. Regulation via Differential Transcript Usage (DTU) is emerging as an important phenomenon. Several tools exist for the detection of DTU from read alignment or assembly data, but options for detection of DTU from alignment-free quantifications are limited. Results We present an R package named RATs – (Relative Abundance of Transcripts) – that identifies DTU transcriptome-wide directly from transcript abundance estimations. RATs is agnostic to quantification methods and exploits bootstrapped quantifications, if available, to inform the significance of detected DTU events. RATs contextualises the DTU results and shows good False Discovery performance (median FDR ≤0.05) at all replication levels. We applied RATs to a human RNA-seq dataset associated with idiopathic pulmonary fibrosis with three DTU events validated by qRT-PCR. RATs found all three genes exhibited statistically significant changes in isoform proportions based on Ensembl v60 annotations, but the DTU for two were not reliably reproduced across bootstrapped quantifications. RATs also identified 500 novel DTU events that are enriched for eleven GO terms related to regulation of the response to stimulus, regulation of immune system processes, and symbiosis/parasitism. Repeating this analysis with the Ensembl v87 annotation showed the isoform abundance profiles of two of the three validated DTU genes changed radically. RATs identified 414 novel DTU events that are enriched for five GO terms, none of which are in common with those previously identified. Only 141 of the DTU evens are common between the two analyses, and only 8 are among the 248 reported by the original study. Furthermore, the original qRT-PCR probes no longer match uniquely to their original transcripts, calling into question the interpretation of these data. We suggest parallel full-length isoform sequencing, annotation pre-filtering and sequencing of the transcripts captured by qRT-PCR primers as possible ways to improve the validation of RNA-seq results in future experiments. Availability The package is available through Github at https://github.com/bartongroup/Rats.

A b s t r a c t

Motivation:
The biological importance of changes in gene and transcript expression is well recognised and is reflected by the wide variety of tools available to characterise these changes. Regulation via Differential Transcript Usage (DTU) is emerging as an important phenomenon. Several tools exist for the detection of DTU from read alignment or assembly data, but options for detection of DTU from alignment-free quantifications are limited. Results: We present an R package named RATs -(Relative Abundance of Transcripts) -that identifies DTU transcriptome-wide directly from transcript abundance estimations. RATs is agnostic to quantification methods and exploits bootstrapped quantifications, if available, to inform the significance of detected DTU events. RATs contextualises the DTU results and shows good False Discovery performance (median FDR ≤ 0.05) at all replication levels. We applied RATs to a human RNA-seq dataset associated with idiopathic pulmonary fibrosis with three DTU events validated by qRT-PCR. RATs found all three genes exhibited statistically significant changes in isoform proportions based on Ensembl v60 annotations, but the DTU for two were not reliably reproduced across bootstrapped quantifications. RATs also identified 500 novel DTU events that are enriched for eleven GO terms related to regulation of the response to stimulus, regulation of immune system processes, and symbiosis/parasitism. Repeating this analysis with the Ensembl v87 annotation showed the isoform abundance profiles of two of the three validated DTU genes changed radically. RATs identified 414 novel DTU events that are enriched for five GO terms, none of which are in common with those previously identified. Only 141 of the DTU evens are common between the two analyses, and only 8 are among the 248 reported by the original study. Furthermore, the original qRT-PCR probes no longer match uniquely to their original transcripts, calling into question the interpretation of these data. We suggest parallel full-length isoform sequencing, annotation pre-filtering and sequencing of the transcripts captured by qRT-PCR primers as possible ways to improve the validation of RNA-seq results in future experiments. Availability: The package is available through Github at https://github.com/bartongroup/Rats .
High-throughput transcriptomics experiments have until recently focused on quantifying gene expression and calculating Differential Gene Expression (DGE) between samples in different groups, conditions, treatments, or time-points. In higher eukaryotes, alternative splicing of multi-exon genes and alternative transcript start and end sites can lead to multiple transcript isoforms originating from each gene. Since isoforms can have different and distinct functions [1][2][3], analysis of Differential Transcript Expression (DTE) is preferable to DGE. Unfortunately, isoform-level transcriptome analysis is more complex and expensive since, in order to achieve similar statistical power in a DTE study, higher sequencing depth is required to compensate for the expression of each gene being split among its component isoforms. In addition, isoforms of a gene typically have a high degree of sequence similarity and this complicates the attribution of reads among them. Despite these challenges, several studies have shown that shifts in individual isoform expression represent a real level of gene regulation with phenotypic consequences [4][5][6][7], suggesting there is little justification for choosing DGE over DTE in the study of complex transcriptomes.
It is possible to find significant DTE among the isoforms of a gene, even when the gene shows no significant DGE. This introduces the concept of Differential Transcript Usage (DTU), where the individual isoform abundances of a gene can change relative to one another, with the most pronounced examples resulting in a change of the dominant isoform (isoform switching). The definitions of DGE, DTE and DTU are illustrated in Figure 1.
To quantify the isoforms and assess changes in their abundance, most existing tools for DTU analysis (e.g. Cufflinks [8], DEXSeq [9], LeafCuttter [10]) rely on reads that either span splice-junctions or align to unique exons. However, with the newest generation of transcript quantification tools (Kallisto [11], Sailfish [12], Salmon [13]), reads are not aligned to the transcriptome or the genome. Instead, these tools combine a pseudo-mapping of the k-mers present within each read to the k-mer distributions from the transcript annotation with an expectation maximization algorithm, to infer the expression of each transcript model directly. Alignment-free methods are much faster than traditional alignment-based methods (RSEM [14], TopHat2 [15], STAR [16]) or assembly-based methods (Cufflinks [8], Trinity [17]), but the lack of alignments prevents these new methods from being compatible with differential expression methods such as Cufflinks, DEXSeq and Leafcutter. Instead, Sleuth [18] is a tool that handles DTE analysis from alignment-free transcript quantifications, but DTU analysis is currently less straight-forward. SwitchSeq [19] and iso-kTSP [6] focus on a particular subset of DTU analysis from alignment-free data, namely isoform switching. SUPPA [20,21], on the other hand, identifies differential splicing events.
In this paper, we present RATs, a tool for identifying DTU directly from isoform quantifications. It is designed primarily for use with alignment-free abundance data and can take advantage of bootstrapped quantifications, but remains agnostic to the quantification method and can also be used in alignment-based workflows. We compare the results from RATs to published and validated instances of DTU and demonstrate that the results of both RNA-seq based and qRT-PCR based analyses are sensitive to the annotation used for transcript quantification and primer design, respectively.  RATs applies filters to the data prior to any statistical testing in order to reduce both the number of low quality calls and the number of tests carried out. These filters are (i) isoform ratio changes can only be defined for genes that are expressed in both conditions and (ii), transcript abundances need to exceed a minimum count threshold. The workflow of RATs is illustrated in Figure 2.
Significant relative transcript abundance changes are detected using two separate approaches: one detects DTU at the gene level and the other detects DTU at the transcript level. The Gtest of independence [22], without continuity corrections, is implemented in RATs and used by both approaches. At the gene level, RATs compares the set of each gene's isoform abundances between two conditions to identify if the abundance ratios have changed. At the transcript level, RATs compares the abundance of each individual transcript against the pooled abundance of its sibling isoforms to identify changes in the proportion of the gene's expression attributable to that specific transcript. Both methods include the Benjamini-Hochberg false discovery rate correction for multiple testing [23], and achieve median false positive rates ≤ 0.05 (FPR = false positives / input negatives) even with only three replicates per condition, with notable improvements at higher replication levels, as shown in Figure 3 (panel A). These two approaches to identifying DTU have different strengths and weaknesses when applied to an RNA-seq dataset. Simultaneously utilizing the expression information across all the isoforms in a gene makes the gene-level test sensitive to smaller changes in relative expression, compared to testing transcripts individually, but also more prone to false positives. Figure 3, shows that the gene-level test has a higher mean and median FPR than the transcript-level test irrespective of replication or effect size, although the two methods converge for highly replicated experiments and/or large effect sizes. Furthermore, it is only applicable to genes with at least two isoforms that both pass the pre-filtering criteria imposed by RATs, potentially limiting its utility for genes with low read coverage. Finally, the genelevel test only identifies the presence of a shift in the ratios of the isoforms belonging to the gene, without identifying which specific isoforms are affected. The transcript-level test, by contrast, directly identifies the specific isoforms whose proportions are changing and has a lower False Positive Rate than the gene-level test. However, considering each isoform independently requires a larger number of tests to be performed, thus resulting in a greater multiple testing penalty and a lower statistical power. RATs currently does not combine the results of the two methods into a single verdict, but rather presents both as complementary views of the data.  [24].

[A] FPR of each bootstrap iteration, for default values of all RATs parameters (v0.4.4), across a range of replication levels. [B] Mean FPR by replication level, as a function of effect size threshold (effect size = difference between an isoform's proportions). For genes, the effect size is defined as the largest proportion difference observed among each gene's isoforms. In every iteration, the FPR was calculated against the number of genes or transcripts that were eligible for testing each time, but that number remained very stable across bootstraps and replication levels (Supplementary file 1). In this case the FPR calculated equates to the False Discovery Rate, because the p-values are appropriately corrected for multiple testing.
Alignment-free transcript quantification is hundreds of times faster than traditional read alignment methods [11,13], making it feasible to repeat the process many times on iterative subsets of the read data. This bootstrapping approach quantifies the technical variance in the transcript abundance estimates. RATs provides the option to use these bootstrapped abundance estimates to apply a reproducibility constraint on the DTU calls, by replacing the mean abundances with the values from random individual iterations and measuring the fraction of these iterations that result in a positive DTU classification. This bootstrapping process reveals the extent to which the variability in the quantification of the expression of each transcript affects the DTU identification of the gene and a threshold can be set on the reproducibility of this classifications. Similarly, RATs also optionally measures the reproducibility of the DTU results relative to the inter-replicate variation by sub-setting the samples. Finally, two further user-defined thresholds are applied: one filters out lowexpression transcripts, while the other excludes transcripts according to the change in their DTU proportion (effect size).
RATs is implemented in R [25]. It is distributed through Github as an R source package (https://github.com/bartongroup/RATS) and has been freely available since August 2016. A detailed manual is included in the package.
RATs accepts as input either a Sleuth object [18] or, alternatively, a set of tables of fragment count estimates (with or without bootstrap information). An annotation mapping the correspondence between transcript and gene identifiers is also required. This can be given directly as a table or inferred from a GTF file.
RATs results are returned in the form of R data.table objects [26]. Along with the DTU calls, the results record the full provenance of the calculations, allowing the user to trace back the decision steps behind each call and gain better understanding of their data. Additionally, RATs employs the ggplot2 package [27] to provide visualisations of the results. In addition to a "volcano plot" view of the DTU effect size against the significance, RATs also visualises the DTU classifications for a given gene in the context of the respective isoform abundances and their consistency across the replicates. An example, drawn from the case study in the following section, is shown in Figure 4. Panel A shows the absolute (upper panels) and relative (lower panels) abundances for each isoform within each condition. Panel B shows the same information, now reorganised to facilitate the comparison of transcript abundance between conditions. From these plots it is straightforward to identify the dominant isoform in each condition and the cases of isoform up-or down-regulation. In this example ENST00000457634 is the dominant transcript from this gene in the control condition but is down-regulated in the IPF condition, while ENST00000490573 is strongly up-regulated in the IPF condition becoming the dominant isoform. The lines connecting transcript abundance points in panel A show the transcript expression profile for this gene for each of the three replicates, making it easy to identify outlier behaviour that may warrant further investigation. In this example, the transcript expression profile from replicate 1 does not agree well with the other two replicates in the control condition. Finally, the DTU classification of the isoforms is encoded as colour or shape, highlighting the changes that are considered significant for the defined criteria. Once created, all plots produced by RATs remain customisable via standard ggplot2 operations.
To test the ability of RATs to identify correctly instances of DTU, we compared it against published and validated instances of DTU from publicly available RNA-seq data. We took RNA-seq data from Deng et al. [28], who identified non-DGE changes in the isoform levels of genes. predictions. The majority of the rejected DTU instances show changes in proportion smaller than our threshold of 0.2, despite high statistical significances, and a also fail to meet our bootstrap reproducibility criterion.
Repeating the quantification and DTU detection using the current version of the human genome assembly and annotation, Ensembl v87, reveals a different picture of DTU in these data. Despite a 25% increase in the number of annotated transcripts in Ensembl v87 compared to v60 (  In light of the expanded annotation, the DTU results of the three genes of interest based on the same RNA-seq data changed considerably (Table 3). TOM1L1 no longer demonstrates a switch in primary isoform from ENST00000445275 in the controls to The updates in the genome assembly and annotation are also likely to impact the result of the qRT-PCR validation of the DTU identified in these genes, since PCR primers are designed to target specific transcripts all in a defined annotation. To examine this, the reported primer sequences, designed according to Ensembl v60, were searched for in the respective Ensembl v87 transcript sequences ( Figure 5). Two pairs of primers were originally used for each gene [28], each pair consisting of a forward and a reverse primer: One pair was common between the two annotated isoforms of each gene and one was unique to one of them. The unique pair used for TOM1L1 targeted isoform ENST00000445275 in the Ensembl v60 annotation. In the Ensembl v87 annotation, however, these primers match two other isoforms (ENST00000570371 and ENST00000575882). The originally targeted ENST00000445275 has been re-annotated with a truncated 5' end in Ensembl v87 and now matches only the reverse primer. Thus, it would not be expected to be detected by the given qRT-PCR design. The reverse primer matches an additional isoform in Ensembl v87 (ENST00000570965) that is also annotated with a short 5' end. Given that isoforms are revised and their ends can change between annotation versions, it is uncertain how many more isoforms may match the primer pair in full or in part. Additionally, the common primers used by Deng et al. [28] to detect all the isoforms of the gene match fewer than half of the annotated isoforms in Ensembl v87, and the isoforms captured by the unique primers are not a subset of those captured by the common primers. Interpreting the given qRT-PCR results in the context of the most recent annotation is thus extremely difficult, if not impossible, and any conclusions reached would certainly differ substantially from those published in the original study. For CMTM4, only one pair of primers is listed by Deng et al. and it matches two isoforms in Ensembl v87, one of which is the originally targeted ENST00000330687. Finally, in PEX11B, the unique primer pair still uniquely matches isoform ENST00000369306 and the common primer pair matches the new isoform as well as the old ones, lending confidence to the original interpretation of the qRT-PCR quantifications for this gene. G e n e T r a n s c r  RATs fills a gap in the line-up of differential expression tools, enabling the identification of DTU from isoform expression quantification data. It is based on established statistical methods for ratio comparisons and provides results in formats suitable for both downstream computational analysis and visual inspection. The package includes plotting and summarization routines designed to encode several layers of information, allowing users to examine quickly and easily both the high level DTU picture of their data and to drill down into the individual details and provenance of each result.
The capability of RATs to reliably identify DTU depends critically on the upstream isoform expression quantification tools and the quality of the input data they use. This limitation is common to all tools that use an annotation, and any downstream analysis can be affected by the choice of annotation [38]. We showed this by applying RATs to publicly available RNAseq data with validated instances of DTU, using two annotation versions separated by six years. All three validated genes in the study contained additional isoforms in the newer annotation and only one of them displayed similar isoform abundance shifts with both annotations. A second gene showed DTU attributable to different isoforms, depending on the annotation version, whereas the third gene showed no significant DTU with the newer annotation. We also demonstrated that qRT-PCR, often considered the gold standard for differential expression validation, is subject to the same limitation, as evidenced by the multiple hits of the previously assumed unique primer sequences in the newer annotation. qRT-PCR primers are designed to target unique sequences in the transcriptome. However, the transcriptome remains a work in progress even for the most intensively studied model organisms. The natural extensive sequence overlap between isoforms together with the ongoing discovery of additional isoforms may mean that, for many genes, qRT-PCR is not a suitable method for the validation of transcript abundance changes (DTE and DTU alike).
For hybridisation-based methods like qRT-PCR to serve as a reliable validation method of RNA quantification, the suitability of the primers should first be validated by sequencing the captured transcripts. Additionally, it has been reported that pre-filtered annotations improve quantification performance [39], and are likely to be helpful in primer design as well. Such annotations could be obtained by including a parallel set of full-length isoform RNA-seq data in the experimental design, such as via PacBio sequencing (http://www.pacb.com) or Oxford Nanopore Direct RNA-seq (https://nanoporetech.com). An additional advantage of this approach is that it effectively defines the relevant transcriptome for the specific experiment [40][41][42][43]. This may be of importance for experiments focussing on specific tissues or developmental stages of an organism where the transcriptome of the sample is likely to be only a subset of the global reference transcriptome of the organism.
Despite these concerns over the annotation, DTU results by RATs were shown to be enriched for GO terms that may be relevant to the symptoms and risk factors of IPF or lung cancer. As the controls in the study were not "healthy" individuals, it is likely that the DTU results from this comparison include genes related to the cancer pathology in addition to those associated with the IPF pathology. It is worth keeping in mind that DTU analysis only considers shifts in the ratios of isoforms and that there may be additional genes and transcripts differentially regulated that are not reported here. In fact, no single analysis type among DGE, DTE and DTU can give a complete picture of differential expression; at least two of these analyses must be carried out.
In conclusion, we offer a computational method for identification and visualisation of differential transcript usage and recommend that caution and scrutiny must be exercised in the interpretation of quantifications, whether they be from RNA-seq or qRT-PCR.