Performing clustering after "seeing" that there are 3 clusters on the PCA scatter plot I have five (normalized) variables related to family history of insurance applicants. n = 59,831.
I performed a Principal Components Reduction, reduced it to 3 dimensions, put together explains 90.1% variance. When i plot the principal component vectors in 3-D, i get what i can see as 3 clear clusters, as below, which are split along the PC2 dimension (1 & 3 cross-section does not see as clear split) (plotted using scatterplot3d library in R):
The question is how do i now assign the labels to these clusters? I kind of 'know' that there are 3 clusters, so this seems to add a supervised flavour to the problem. I tried k-means along just the PC2 dimension, and also along all 3 dimensions but even with a lot of nstarts (=1000) the clusters dont seem to come as expected: below is clustered along PC2 dimension using K-means = 3, I am guessing that K-means is trying to divide the clusters 'equally' (same number of n in each cluster?)

And below is clustered using kmeans=3 , along PC1 and PC2 dimensions (similar results for all 3 dimensions). The new split along PC1 suggests that adding more dimensions to the clustering exercise may not be a good idea.

So, using just PC2 as the clustering variable, how would i get the correct clusters? Thanks in advance...
UPDATE & EDIT - Thanks to good suggestions, I tried applying the GMM with k = 3. I am getting results like the below. These are much better than using k-means. However, the 'ends' of the clusters dont seem to be assimilated within their main distributions..I'm getting the feeling that all the 'outliers' are getting their own cluster here. It seems to me that none of the surrounding distributions wants to take 'responsibility' for these red points in between the green and black clusters.

I used the mvnormalmixEM and the normalmixEM (for just PC2 dimension) to get the posterior probabilities after convergence, and then assigned the label where the posterior probability was the highest among the three GMM components. I didn't set any initial values of the prior responsibilities/weights [lambda parameter]. Would anyone have ideas on how to fine-tune? 
It feels to me that the red points arent falling within the distribution parameters that the EM algorithm has discovered, and that I might need to suggest to the algorithm to search for larger variance parameters, but I am unsure how to do this
Thanks in advance
 A: K means is very sensitive to feature scaling. So it's hard to get mathematically meaningful results when your axes have different units!
Furthermore, the PCA plot is likely misleading. It is probably showing some artifact of your data. There probably is some attribute with very high variance and just three levels (e.g. income 1000, 10000, 100000) that dominates this visualization in an undue way. The other attributes then just add a "gaussian blur" to this. So all you are doing is reverse engineer that attribute that already is in your data. My guess is that you can simply identify it by looking at the highest variance attribute.
A: Finally got the below using GaussianMixture in python's scikit-learn v0.19. This allowed me to specify that the covariance matrix should be shared across clusters (option = 'tied' in the covariance_type param) - something which doesn't seem available as an option in R.
I also perhaps, 'cheated' a little bit by roughly specifying the initial distribution weights and the initial rough centroid locations by visual inspection of the PCA scatterplot. 

A: It will be quite useful to look also at the loadings. Which is the variable that create the scatter plot you see ? I always do not trust only scatter plots. Changing axes you can find cluster almost wherever you would like to and for sure you can find a methodology that can give you 3 clusters. 
