To evaluate the milk thistle transcriptome under drought stress in field conditions, irrigation was applied using a weighted method at three levels: 100% F.C every 2 days, 70% F.C every 4 days, and 40% F.C every 8 days. Sampling was performed after 8 days at the flowering stage. Plant leaves were collected for RNA-seq analysis, seeds for oily and methanolic extracts, and downstream analyses were performed. Since there was no annotated reference genome for this plant, the De novo Assembly method was implemented to assemble the transcriptome. Contigs were blasted against five databases: NT, NR, Uniprot, and protein databases of Arabidopsis thaliana and Helianthus annuus. A total of 9,517 genes (~ 73% of Uniprot genes) were common across all databases and selected for further analysis due to their comprehensive annotation. Then, DEGs were identified and functionally annotated using Gene Ontology (GO) analysis with the ShinyGO platform, biological pathway analysis through KEGG, and transcription factor identification via PlantTFDB. Next, silybinin content was measured using HPLC. Generally, the most repeated pathways in all treatments include the Biosynthesis of secondary metabolites and the MAPK signaling pathway. Also, most biological processes are related to the oxidation-reduction process, and response to stress, and most molecular functions are protein and mRNA binding. Our results indicate the active role of transcription factors ERF, C3H, and bHLH in drought stress tolerance. Silybin a and b showed that severe drought stress enhanced the accumulation of silybinin compared with seeds from the control. Eight differentially expressed genes (CYP86A1, CYP710A1, FATA2, LACS3, LOX2, PAL, PLA2-ALPHA, and PXG3) were used to validate the RNA-Seq data. qRT-PCR results confirmed strong consistency with the RNA-Seq findings. Finally, the genes involved in the silymarin pathway were identified, and their expression was determined through RNA-Seq data and compared with the silymarin contents.