Dear user, Finding divergent sites between two species is what we did for D. mojavensis and D. arizonae in our paper. Like you, we had two paired-end samples for each species. 

Below is an example of the KisSplice output we obtain: >bcc_9|Cycle_0|Type_0a|upper_path_Length_83|C1_0|C2_0|C3_0|C4_0|C5_86|C6_187|C7_102|C8_161|Q1_0|Q2_0|Q3_0|Q4_0|Q5_60|Q6_59|Q7_58|Q8_59|rank_1.00000 CCTCAAGCAATGTCCATGCGACTTATTGAAAAATTTACACTGGAATTAGTCCAAGCGTCCATTGAGAACGACAAAAGACTGCA >bcc_9|Cycle_0|Type_0a|lower_path_Length_83|C1_80|C2_165|C3_100|C4_168|C5_0|C6_0|C7_0|C8_0|Q1_59|Q2_59|Q3_61|Q4_58|Q5_0|Q6_0|Q7_0|Q8_0|rank_1.00000 CCTCAAGCAATGTCCATGCGACTTATTGAAAAATTTACACTAGAATTAGTCCAAGCGTCCATTGAGAACGACAAAAGACTGCA

The "SNP" is G/A. The upper path corresponds to G. The lower path to A. C1 to C4 correspond to read counts in mojavensis. C5 to C8 correspond to read counts in arizonae. G is supported by 0 reads in mojavensis, and by (86+187, 102+161) reads in arizonae. A is supported by ~(80+165, 100+168) reads in mojavensis and 0 reads in arizonae. This "SNP" is clearly species specific. One allele is present in one species. The other allele is present in the other species. If you test it in KissDE, you will find a p-value close to 0, and a Deltaf (difference of allele frequency) of 1