Waterlily (Nymphaea spp.), a perennial herbaceous aquatic plant, is divided into two ecological groups: hardy waterlily and tropical waterlily. Although the hardy waterlily has no attractive blue flower cultivar, its adaptability is stronger than tropical waterlily because it can survive a cold winter. Thus, breeding hardy waterlily with real blue flowers has become an important target for breeders. Molecular breeding may be a useful way. However, molecular studies on waterlily are limited due to the lack of sequence data. In this study, six cDNA libraries generated from the petals of two different coloring stages of blue tropical waterlily cultivar Nymphaea ‘King of Siam’ were sequenced using the Illumina HiSeq™ 2500 platform. Each library produced no less than 5.65 Gb clean reads. Subsequently, de novo assembly generated 112,485 unigenes, including 26,206 unigenes annotated to seven public protein databases. Then, 127 unigenes could be identified as putative homologues of color-related genes in other species, including 28 up-regulated and 5 down-regulated unigenes. In petals, 16 flavonoids (4 anthocyanins and 12 flavonols) were detected in different contents during the color development due to the different expression levels of color-related genes, and four flavonols were detected in waterlily for the first time. Furthermore, UA3GTs were selected as the most important candidates involved in the flavonoid metabolic pathway, UA3GTs induced blue petal color formation in Nymphaea ‘King of Siam’. This study will improve our understanding of the molecular mechanism of blue flowers in waterlily and provide the basis for molecular breeding of blue hardy waterlily cultivars.